diff --git a/.buildkite/Manifest.toml b/.buildkite/Manifest.toml
index 75e6dbcbd4..bb70456681 100644
--- a/.buildkite/Manifest.toml
+++ b/.buildkite/Manifest.toml
@@ -373,7 +373,7 @@ weakdeps = ["CUDA", "Krylov"]
     KrylovExt = "Krylov"
 
 [[deps.ClimaLand]]
-deps = ["Adapt", "ArtifactWrappers", "ClimaComms", "ClimaCore", "ClimaUtilities", "DataFrames", "Dates", "DocStringExtensions", "Insolation", "Interpolations", "IntervalSets", "LinearAlgebra", "NCDatasets", "SciMLBase", "StaticArrays", "SurfaceFluxes", "Thermodynamics"]
+deps = ["Adapt", "ArtifactWrappers", "ClimaComms", "ClimaCore", "ClimaUtilities", "DataFrames", "Dates", "DocStringExtensions", "Insolation", "Interpolations", "IntervalSets", "LinearAlgebra", "NCDatasets", "SciMLBase", "StaticArrays", "SurfaceFluxes", "Thermodynamics", "UnrolledUtilities"]
 path = ".."
 uuid = "08f4d4ce-cf43-44bb-ad95-9d2d5f413532"
 version = "0.12.1"
@@ -2669,6 +2669,11 @@ git-tree-sha1 = "6cc9d682755680e0f0be87c56392b7651efc2c7b"
 uuid = "9602ed7d-8fef-5bc8-8597-8f21381861e8"
 version = "0.1.5"
 
+[[deps.UnrolledUtilities]]
+git-tree-sha1 = "b73f7a7c25a2618c5052c80ed32b07e471cc6cb0"
+uuid = "0fe1646c-419e-43be-ac14-22321958931b"
+version = "0.1.2"
+
 [[deps.UnsafeAtomics]]
 git-tree-sha1 = "6331ac3440856ea1988316b46045303bef658278"
 uuid = "013be700-e6cd-48c3-b4a1-df204f14c38f"
diff --git a/Project.toml b/Project.toml
index b7bdede713..1fb822dc87 100644
--- a/Project.toml
+++ b/Project.toml
@@ -21,6 +21,7 @@ SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 SurfaceFluxes = "49b00bb7-8bd4-4f2b-b78c-51cd0450215f"
 Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c"
+UnrolledUtilities = "0fe1646c-419e-43be-ac14-22321958931b"
 
 [weakdeps]
 ClimaParams = "5c42b081-d73a-476f-9059-fd94b934656c"
@@ -59,5 +60,6 @@ StaticArrays = "1"
 StatsBase = "0.33, 0.34"
 SurfaceFluxes = "0.10, 0.11"
 Thermodynamics = "0.12"
+UnrolledUtilities = "0.1"
 cuDNN = "1"
 julia = "1.9"
diff --git a/docs/Manifest.toml b/docs/Manifest.toml
index b0832a53d1..ea9e697546 100644
--- a/docs/Manifest.toml
+++ b/docs/Manifest.toml
@@ -370,7 +370,7 @@ weakdeps = ["CUDA", "Krylov"]
     KrylovExt = "Krylov"
 
 [[deps.ClimaLand]]
-deps = ["Adapt", "ArtifactWrappers", "ClimaComms", "ClimaCore", "ClimaUtilities", "DataFrames", "Dates", "DocStringExtensions", "Insolation", "Interpolations", "IntervalSets", "LinearAlgebra", "NCDatasets", "SciMLBase", "StaticArrays", "SurfaceFluxes", "Thermodynamics"]
+deps = ["Adapt", "ArtifactWrappers", "ClimaComms", "ClimaCore", "ClimaUtilities", "DataFrames", "Dates", "DocStringExtensions", "Insolation", "Interpolations", "IntervalSets", "LinearAlgebra", "NCDatasets", "SciMLBase", "StaticArrays", "SurfaceFluxes", "Thermodynamics", "UnrolledUtilities"]
 path = ".."
 uuid = "08f4d4ce-cf43-44bb-ad95-9d2d5f413532"
 version = "0.12.1"
@@ -2692,6 +2692,11 @@ git-tree-sha1 = "6cc9d682755680e0f0be87c56392b7651efc2c7b"
 uuid = "9602ed7d-8fef-5bc8-8597-8f21381861e8"
 version = "0.1.5"
 
+[[deps.UnrolledUtilities]]
+git-tree-sha1 = "b73f7a7c25a2618c5052c80ed32b07e471cc6cb0"
+uuid = "0fe1646c-419e-43be-ac14-22321958931b"
+version = "0.1.2"
+
 [[deps.UnsafeAtomics]]
 git-tree-sha1 = "6331ac3440856ea1988316b46045303bef658278"
 uuid = "013be700-e6cd-48c3-b4a1-df204f14c38f"
diff --git a/docs/src/APIs/Soil.md b/docs/src/APIs/Soil.md
index 8eae45f61f..acc45e5975 100644
--- a/docs/src/APIs/Soil.md
+++ b/docs/src/APIs/Soil.md
@@ -96,5 +96,5 @@ ClimaLand.Soil.RootExtraction
 ## Soil Jacobian Structures
 
 ```@docs
-ClimaLand.Soil.RichardsTridiagonalW
-```
\ No newline at end of file
+ClimaLand.Soil.ImplicitEquationJacobian
+```
diff --git a/docs/src/APIs/shared_utilities.md b/docs/src/APIs/shared_utilities.md
index d09551d3c3..e316e296b9 100644
--- a/docs/src/APIs/shared_utilities.md
+++ b/docs/src/APIs/shared_utilities.md
@@ -61,8 +61,7 @@ ClimaLand.boundary_var_domain_names
 ClimaLand.boundary_var_types
 ClimaLand.make_tendency_jacobian
 ClimaLand.make_update_jacobian
-ClimaLand.∂tendencyBC∂Y
-ClimaLand.AbstractTridiagonalW
+ClimaLand.set_dfluxBCdY!
 ClimaLand.get_drivers
 ```
 
diff --git a/docs/tutorials/standalone/Soil/layered_soil.jl b/docs/tutorials/standalone/Soil/layered_soil.jl
index 40945f9169..8a602e2d9f 100644
--- a/docs/tutorials/standalone/Soil/layered_soil.jl
+++ b/docs/tutorials/standalone/Soil/layered_soil.jl
@@ -29,7 +29,7 @@ t0 = Float64(0)
 tf = Float64(60 * 60)
 dt = Float64(30);
 
-# Define the domain 
+# Define the domain
 zmax = FT(0)
 zmin = FT(-1.1)
 nelems = 75
@@ -109,7 +109,7 @@ soil = Soil.RichardsModel{FT}(;
 Y, p, cds = initialize(soil);
 
 # Initial conditions
-Y.soil.ϑ_l .= 0.0353; # read from Figure 4 of Huang et al. 
+Y.soil.ϑ_l .= 0.0353; # read from Figure 4 of Huang et al.
 
 # We also set the initial conditions of the auxiliary state here:
 set_initial_cache! = make_set_initial_cache(soil)
@@ -131,9 +131,10 @@ ode_algo = CTS.IMEXAlgorithm(
 )
 exp_tendency! = make_exp_tendency(soil)
 imp_tendency! = make_imp_tendency(soil)
-update_jacobian! = make_update_jacobian(soil)
+tendency_jacobian! = make_tendency_jacobian(soil)
+
 jac_kwargs =
-    (; jac_prototype = RichardsTridiagonalW(Y), Wfact = update_jacobian!)
+    (; jac_prototype = ImplicitEquationJacobian(Y), Wfact = tendency_jacobian!)
 prob = SciMLBase.ODEProblem(
     CTS.ClimaODEFunction(
         T_exp! = exp_tendency!,
diff --git a/docs/tutorials/standalone/Soil/richards_equation.jl b/docs/tutorials/standalone/Soil/richards_equation.jl
index e7663e7abd..1eecb1144b 100644
--- a/docs/tutorials/standalone/Soil/richards_equation.jl
+++ b/docs/tutorials/standalone/Soil/richards_equation.jl
@@ -140,7 +140,7 @@ soil = Soil.RichardsModel{FT}(;
 # We also create the function which is used to update our Jacobian.
 exp_tendency! = make_exp_tendency(soil);
 imp_tendency! = ClimaLand.make_imp_tendency(soil);
-update_jacobian! = ClimaLand.make_update_jacobian(soil);
+tendency_jacobian! = ClimaLand.make_tendency_jacobian(soil);
 
 # # Set up the simulation
 # We can now initialize the prognostic and auxiliary variable vectors, and take
@@ -186,7 +186,7 @@ ode_algo = CTS.IMEXAlgorithm(
 
 # Here we set up the information used for our Jacobian.
 jac_kwargs =
-    (; jac_prototype = RichardsTridiagonalW(Y), Wfact = update_jacobian!);
+    (; jac_prototype = ImplicitEquationJacobian(Y), Wfact = tendency_jacobian!);
 
 # And then we can solve the system of equations, using
 # [SciMLBase.jl](https://github.com/SciML/SciMLBase.jl) and
diff --git a/experiments/standalone/Soil/richards_comparison.jl b/experiments/standalone/Soil/richards_comparison.jl
index 76376124c7..51ce4e8330 100644
--- a/experiments/standalone/Soil/richards_comparison.jl
+++ b/experiments/standalone/Soil/richards_comparison.jl
@@ -81,7 +81,7 @@ bonan_sand_dataset = ArtifactWrapper(
         Y.soil.ϑ_l .= FT(0.24)
         exp_tendency! = make_exp_tendency(soil)
         imp_tendency! = ClimaLand.make_imp_tendency(soil)
-        update_jacobian! = ClimaLand.make_update_jacobian(soil)
+        tendency_jacobian! = ClimaLand.make_tendency_jacobian(soil)
         set_initial_cache!(p, Y, t0)
 
         stepper = CTS.ARS111()
@@ -98,8 +98,8 @@ bonan_sand_dataset = ArtifactWrapper(
 
         # set up jacobian info
         jac_kwargs = (;
-            jac_prototype = RichardsTridiagonalW(Y),
-            Wfact = update_jacobian!,
+            jac_prototype = ImplicitEquationJacobian(Y),
+            Wfact = tendency_jacobian!,
         )
 
         prob = SciMLBase.ODEProblem(
@@ -181,7 +181,7 @@ end
         Y.soil.ϑ_l .= FT(0.1)
         exp_tendency! = make_exp_tendency(soil)
         imp_tendency! = ClimaLand.make_imp_tendency(soil)
-        update_jacobian! = ClimaLand.make_update_jacobian(soil)
+        tendency_jacobian! = ClimaLand.make_tendency_jacobian(soil)
         set_initial_cache!(p, Y, t0)
 
         stepper = CTS.ARS111()
@@ -197,8 +197,8 @@ end
         )
         # set up jacobian info
         jac_kwargs = (;
-            jac_prototype = RichardsTridiagonalW(Y),
-            Wfact = update_jacobian!,
+            jac_prototype = ImplicitEquationJacobian(Y),
+            Wfact = tendency_jacobian!,
         )
 
         prob = SciMLBase.ODEProblem(
diff --git a/experiments/standalone/Soil/richards_runoff.jl b/experiments/standalone/Soil/richards_runoff.jl
index 2c7894acc6..820763a112 100644
--- a/experiments/standalone/Soil/richards_runoff.jl
+++ b/experiments/standalone/Soil/richards_runoff.jl
@@ -218,7 +218,7 @@ Y.soil.ϑ_l .= hydrostatic_profile.(lat, z, ν, θ_r, vg_α, vg_n, S_s, f_max)
 set_initial_cache! = make_set_initial_cache(model)
 exp_tendency! = make_exp_tendency(model);
 imp_tendency! = ClimaLand.make_imp_tendency(model);
-update_jacobian! = ClimaLand.make_update_jacobian(model);
+tendency_jacobian! = ClimaLand.make_tendency_jacobian(model);
 
 set_initial_cache!(p, Y, t0)
 stepper = CTS.ARS111()
@@ -235,7 +235,7 @@ ode_algo = CTS.IMEXAlgorithm(
 
 # set up jacobian info
 jac_kwargs =
-    (; jac_prototype = RichardsTridiagonalW(Y), Wfact = update_jacobian!)
+    (; jac_prototype = ImplicitEquationJacobian(Y), Wfact = tendency_jacobian!)
 
 prob = SciMLBase.ODEProblem(
     CTS.ClimaODEFunction(
diff --git a/experiments/standalone/Soil/water_conservation.jl b/experiments/standalone/Soil/water_conservation.jl
index ab862d3350..6559a4e2f8 100644
--- a/experiments/standalone/Soil/water_conservation.jl
+++ b/experiments/standalone/Soil/water_conservation.jl
@@ -102,7 +102,7 @@ for FT in (Float32, Float64)
     exp_tendency! = make_exp_tendency(soil)
     set_initial_cache! = make_set_initial_cache(soil)
     imp_tendency! = make_imp_tendency(soil)
-    update_jacobian! = make_update_jacobian(soil)
+    tendency_jacobian! = make_tendency_jacobian(soil)
 
     rmses = Array{FT}(undef, length(dts))
     mass_errors = Array{FT}(undef, length(dts))
@@ -114,8 +114,8 @@ for FT in (Float32, Float64)
         set_initial_cache!(p, Y, FT(0.0))
 
         jac_kwargs = (;
-            jac_prototype = RichardsTridiagonalW(Y),
-            Wfact = update_jacobian!,
+            jac_prototype = ImplicitEquationJacobian(Y),
+            Wfact = tendency_jacobian!,
         )
 
         prob = SciMLBase.ODEProblem(
@@ -143,6 +143,7 @@ for FT in (Float32, Float64)
             mass_change_exp = -(flux_in - flux_out) * t_sim
             mass_change_actual = mass_end - mass_start
             relerr = abs(mass_change_actual - mass_change_exp) / mass_change_exp
+            @show relerr
             @assert relerr < FT(1e-9)
             mass_errors[i] = relerr
 
@@ -212,7 +213,7 @@ for FT in (Float32, Float64)
     exp_tendency! = make_exp_tendency(soil_dirichlet)
     set_initial_cache! = make_set_initial_cache(soil_dirichlet)
     imp_tendency! = make_imp_tendency(soil_dirichlet)
-    update_jacobian! = make_update_jacobian(soil_dirichlet)
+    tendency_jacobian! = make_tendency_jacobian(soil_dirichlet)
     update_aux! = make_update_aux(soil_dirichlet)
 
     rmses_dirichlet = Array{FT}(undef, length(dts))
@@ -225,8 +226,8 @@ for FT in (Float32, Float64)
         set_initial_cache!(p, Y, FT(0.0))
 
         jac_kwargs = (;
-            jac_prototype = RichardsTridiagonalW(Y),
-            Wfact = update_jacobian!,
+            jac_prototype = ImplicitEquationJacobian(Y),
+            Wfact = tendency_jacobian!,
         )
 
         prob = SciMLBase.ODEProblem(
diff --git a/src/ClimaLand.jl b/src/ClimaLand.jl
index 4fe4c699d9..d53db2bf3b 100644
--- a/src/ClimaLand.jl
+++ b/src/ClimaLand.jl
@@ -23,7 +23,6 @@ include("shared_utilities/drivers.jl")
 include("shared_utilities/boundary_conditions.jl")
 include("shared_utilities/sources.jl")
 include("shared_utilities/implicit_tendencies.jl")
-include("shared_utilities/implicit_functions.jl")
 include("standalone/Bucket/Bucket.jl")
 
 """
diff --git a/src/shared_utilities/Domains.jl b/src/shared_utilities/Domains.jl
index 615d92e216..371eb3be59 100644
--- a/src/shared_utilities/Domains.jl
+++ b/src/shared_utilities/Domains.jl
@@ -658,29 +658,28 @@ extruded center finite difference space, this would return a 2D field
 corresponding to the surface, with values equal to the topmost level.
 """
 function top_center_to_surface(center_field::ClimaCore.Fields.Field)
-    cs = axes(center_field)
-    face_space = obtain_face_space(cs)
-    N = ClimaCore.Spaces.nlevels(face_space)
-    interp_c2f = ClimaCore.Operators.InterpolateC2F(
-        top = ClimaCore.Operators.Extrapolate(),
-        bottom = ClimaCore.Operators.Extrapolate(),
-    )
-    surface_space = obtain_surface_space(cs)
+    center_space = axes(center_field)
+    N_minus_half = ClimaCore.Spaces.nlevels(center_space)
+    surface_space = obtain_surface_space(center_space)
     return ClimaCore.Fields.Field(
         ClimaCore.Fields.field_values(
-            ClimaCore.Fields.level(
-                interp_c2f.(center_field),
-                ClimaCore.Utilities.PlusHalf(N - 1),
-            ),
+            ClimaCore.Fields.level(center_field, N_minus_half),
         ),
         surface_space,
     )
 end
 
+"""
+    top_center_to_surface(val)
+
+When `val` is a scalar (e.g. a single float or struct), returns `val`.
+"""
+top_center_to_surface(val) = val
+
 """
     linear_interpolation_to_surface!(sfc_field, center_field, z)
 
-Linearly interpolate the center field `center_field` to the surface 
+Linearly interpolate the center field `center_field` to the surface
 defined by the top face coordinate of `z`; updates the `sfc_field`
 on the surface (face) space in place.
 """
@@ -705,7 +704,7 @@ end
     get_Δz(z::ClimaCore.Fields.Field)
 
 A function to return a tuple containing the distance between the top boundary
-and its closest center, and the bottom boundary and its closest center, 
+and its closest center, and the bottom boundary and its closest center,
 both as Fields.
 """
 function get_Δz(z::ClimaCore.Fields.Field)
diff --git a/src/shared_utilities/implicit_functions.jl b/src/shared_utilities/implicit_functions.jl
deleted file mode 100644
index 16daccfbae..0000000000
--- a/src/shared_utilities/implicit_functions.jl
+++ /dev/null
@@ -1,125 +0,0 @@
-### These functions are here for now, but linsolve!, ldiv!, and
-### thomas_algorithm! will soon be moved to ClimaCore.
-
-import LinearAlgebra
-import ClimaCore: Spaces
-
-export linsolve!, thomas_algorithm!
-
-# Function required by SciMLBase.jl
-linsolve!(::Type{Val{:init}}, f, u0; kwargs...) = _linsolve!
-_linsolve!(x, A, b, update_matrix = false; kwargs...) =
-    LinearAlgebra.ldiv!(x, A, b)
-
-# Function required by Krylov.jl (x and b can be AbstractVectors)
-# See https://github.com/JuliaSmoothOptimizers/Krylov.jl/issues/605 for a
-# related issue that requires the same workaround.
-function LinearAlgebra.ldiv!(x, A::AbstractTridiagonalW, b)
-    A.temp1 .= b
-    LinearAlgebra.ldiv!(A.temp2, A, A.temp1)
-    x .= A.temp2
-end
-
-function LinearAlgebra.ldiv!(
-    x::Fields.FieldVector,
-    A::AbstractTridiagonalW,
-    b::Fields.FieldVector,
-)
-    _ldiv!(x, A, b, axes(x.soil.ϑ_l))
-end
-
-# 2D space - only one column
-function _ldiv!(
-    x::Fields.FieldVector, # Δx
-    A::AbstractTridiagonalW, # f'(x)
-    b::Fields.FieldVector, # -f(x)
-    space::Spaces.FiniteDifferenceSpace,
-)
-    (; dtγ_ref, ∂ϑₜ∂ϑ, W_column_arrays, transform) = A
-    dtγ = dtγ_ref[]
-
-    _ldiv_serial!(
-        x.soil.ϑ_l,
-        b.soil.ϑ_l,
-        dtγ,
-        transform,
-        ∂ϑₜ∂ϑ,
-        W_column_arrays[Threads.threadid()], # can / should this be colidx? TODO do we need this in single-column case
-    )
-end
-
-# 3D space - iterate over columns
-function _ldiv!(
-    x::Fields.FieldVector,
-    A::AbstractTridiagonalW,
-    b::Fields.FieldVector,
-    space::Spaces.ExtrudedFiniteDifferenceSpace,
-)
-    (; dtγ_ref, ∂ϑₜ∂ϑ, W_column_arrays, transform) = A
-    dtγ = dtγ_ref[]
-
-    Fields.bycolumn(space) do colidx
-        _ldiv_serial!(
-            x.soil.ϑ_l[colidx],
-            b.soil.ϑ_l[colidx],
-            dtγ,
-            transform,
-            ∂ϑₜ∂ϑ[colidx],
-            W_column_arrays[Threads.threadid()], # can / should this be colidx?
-        )
-    end
-end
-
-function _ldiv_serial!(
-    x_column,
-    b_column,
-    dtγ,
-    transform,
-    ∂ϑₜ∂ϑ_column,
-    W_column_array,
-)
-    x_column .= b_column
-
-    x_column_view = parent(x_column)
-
-    @views W_column_array.dl .= dtγ .* parent(∂ϑₜ∂ϑ_column.coefs.:1)[2:end]
-    W_column_array.d .= -1 .+ dtγ .* parent(∂ϑₜ∂ϑ_column.coefs.:2)
-    @views W_column_array.du .=
-        dtγ .* parent(∂ϑₜ∂ϑ_column.coefs.:3)[1:(end - 1)]
-    thomas_algorithm!(W_column_array, x_column_view)
-
-    # Apply transform (if needed)
-    if transform
-        x_column .*= dtγ
-    end
-    return nothing
-end
-
-"""
-    thomas_algorithm!(A, b)
-
-Thomas algorithm for solving a linear system A x = b,
-where A is a tri-diagonal matrix.
-A and b are overwritten, solution is written to b.
-Pass this as linsolve to ODEFunction.
-"""
-function thomas_algorithm!(A, b)
-    nrows = size(A, 1)
-    # first row
-    @inbounds A[1, 2] /= A[1, 1]
-    @inbounds b[1] /= A[1, 1]
-    # interior rows
-    for row in 2:(nrows - 1)
-        @inbounds fac = A[row, row] - (A[row, row - 1] * A[row - 1, row])
-        @inbounds A[row, row + 1] /= fac
-        @inbounds b[row] = (b[row] - A[row, row - 1] * b[row - 1]) / fac
-    end
-    # last row
-    @inbounds fac = A[nrows, nrows] - A[nrows - 1, nrows] * A[nrows, nrows - 1]
-    @inbounds b[nrows] = (b[nrows] - A[nrows, nrows - 1] * b[nrows - 1]) / fac
-    # back substitution
-    for row in (nrows - 1):-1:1
-        @inbounds b[row] -= b[row + 1] * A[row, row + 1]
-    end
-    return nothing
-end
diff --git a/src/shared_utilities/implicit_tendencies.jl b/src/shared_utilities/implicit_tendencies.jl
index 9b22348a0e..3a0c5b13e2 100644
--- a/src/shared_utilities/implicit_tendencies.jl
+++ b/src/shared_utilities/implicit_tendencies.jl
@@ -1,5 +1,4 @@
-export make_tendency_jacobian,
-    make_update_jacobian, AbstractTridiagonalW, ∂tendencyBC∂Y
+export make_tendency_jacobian, make_update_jacobian, set_dfluxBCdY!
 
 
 """
@@ -48,7 +47,7 @@ function make_update_jacobian(model::AbstractModel)
 end
 
 """
-    ∂tendencyBC∂Y(::AbstractModel,
+    set_dfluxBCdY!(::AbstractModel,
                   ::AbstractBC,
                   ::AbstractBoundary,
                   _...)::Union{ClimaCore.Fields.FieldVector, Nothing}
@@ -57,18 +56,9 @@ A function stub which returns the derivative of the implicit tendency
 term of the `model` arising from the boundary condition,
 with respect to the state Y.
 """
-function ∂tendencyBC∂Y(
+function set_dfluxBCdY!(
     ::AbstractModel,
     ::AbstractBC,
     ::AbstractBoundary,
     _...,
 )::Union{ClimaCore.Fields.FieldVector, Nothing} end
-
-"""
-    AbstractTridiagonalW
-
-An abstract type for tridiagonal Jacobian matrices.
-"""
-abstract type AbstractTridiagonalW end
-
-Base.similar(w::AbstractTridiagonalW) = w
diff --git a/src/shared_utilities/models.jl b/src/shared_utilities/models.jl
index f9ccbcf9c6..2dc95de78c 100644
--- a/src/shared_utilities/models.jl
+++ b/src/shared_utilities/models.jl
@@ -164,7 +164,7 @@ end
      make_update_cache(model::AbstractModel)
 
 A helper function which updates all cache variables of a model;
-currently only used in `set_initial_cache` since not all 
+currently only used in `set_initial_cache` since not all
 cache variables are updated at the same time.
 """
 function make_update_cache(model::AbstractModel)
@@ -404,7 +404,6 @@ function get_drivers(model::AbstractModel)
     return (nothing, nothing)
 end
 
-
 """
     initialize(model::AbstractModel)
 
diff --git a/src/shared_utilities/utils.jl b/src/shared_utilities/utils.jl
index aec5eef230..de69945e93 100644
--- a/src/shared_utilities/utils.jl
+++ b/src/shared_utilities/utils.jl
@@ -16,14 +16,6 @@ function heaviside(x::FT)::FT where {FT}
     end
 end
 
-"""
-    to_scalar_coefs(vector_coefs)
-
-Helper function used in computing tendencies of vertical diffusion terms.
-"""
-to_scalar_coefs(vector_coefs) = map(vector_coef -> vector_coef.u₃, vector_coefs)
-
-
 """
     add_dss_buffer_to_aux(p::NamedTuple, domain::Domains.AbstractDomain)
 
diff --git a/src/standalone/Soil/Soil.jl b/src/standalone/Soil/Soil.jl
index a03969868c..1f6c3b5296 100644
--- a/src/standalone/Soil/Soil.jl
+++ b/src/standalone/Soil/Soil.jl
@@ -64,7 +64,6 @@ import ClimaCore: Fields, Operators, Geometry, Spaces
 using Thermodynamics
 
 import ClimaLand.Domains: Column, HybridBox, SphericalShell
-using ClimaLand: AbstractTridiagonalW
 import ClimaLand:
     AbstractImExModel,
     make_update_aux,
@@ -91,7 +90,7 @@ import ClimaLand:
     get_drivers
 export RichardsModel,
     RichardsParameters,
-    RichardsTridiagonalW,
+    ImplicitEquationJacobian,
     EnergyHydrology,
     EnergyHydrologyParameters,
     AbstractSoilModel,
diff --git a/src/standalone/Soil/boundary_conditions.jl b/src/standalone/Soil/boundary_conditions.jl
index e2e9f9f108..b1d50e2554 100644
--- a/src/standalone/Soil/boundary_conditions.jl
+++ b/src/standalone/Soil/boundary_conditions.jl
@@ -5,6 +5,8 @@ import ClimaLand:
     boundary_vars,
     boundary_var_domain_names,
     boundary_var_types
+using ClimaLand: Domains
+using ClimaCore: Geometry
 export TemperatureStateBC,
     MoistureStateBC,
     FreeDrainage,
@@ -161,7 +163,7 @@ end
                                               <:Runoff.TOPMODELRunoff,
                                               }, ::ClimaLand.TopBoundary)
 
-An extension of the `boundary_vars` method for RichardsAtmosDrivenFluxBC with 
+An extension of the `boundary_vars` method for RichardsAtmosDrivenFluxBC with
 TOPMODELRunoff.
 
 These variables are updated in place in `boundary_flux`.
@@ -181,7 +183,7 @@ boundary_vars(
                               ::ClimaLand.TopBoundary)
 
 An extension of the `boundary_var_domain_names` method for RichardsAtmosDrivenFluxBC
-with TOPMODELRunoff. 
+with TOPMODELRunoff.
 """
 boundary_var_domain_names(
     bc::RichardsAtmosDrivenFluxBC{
@@ -191,11 +193,12 @@ boundary_var_domain_names(
     ::ClimaLand.TopBoundary,
 ) = (:surface, :surface, :surface, :surface, :surface)
 """
-    boundary_var_types(::RichardsAtmosDrivenFluxBC{<:PrescribedPrecipitation,
+    boundary_var_types(::RichardsModel{FT},
+                        ::RichardsAtmosDrivenFluxBC{<:PrescribedPrecipitation,
                                                    <:Runoff.TOPMODELRunoff{FT},
                                                   },
-                       ::ClimaLand.TopBoundary,
-                       ) where {FT}
+                        ::ClimaLand.TopBoundary,
+                        ) where {FT}
 
 An extension of the `boundary_var_types` method for RichardsAtmosDrivenFluxBC
 with TOPMODELRunoff.
@@ -214,7 +217,7 @@ boundary_var_types(
                                               <:Runoff.AbstractRunoffModel,
                                               }, ::ClimaLand.TopBoundary)
 
-An extension of the `boundary_vars` method for RichardsAtmosDrivenFluxBC with 
+An extension of the `boundary_vars` method for RichardsAtmosDrivenFluxBC with
 no runoff modeled.
 
 These variables are updated in place in `boundary_flux`.
@@ -234,7 +237,7 @@ boundary_vars(
                               ::ClimaLand.TopBoundary)
 
 An extension of the `boundary_var_domain_names` method for RichardsAtmosDrivenFluxBC
-with no runoff modeled. 
+with no runoff modeled.
 """
 boundary_var_domain_names(
     bc::RichardsAtmosDrivenFluxBC{
@@ -244,11 +247,12 @@ boundary_var_domain_names(
     ::ClimaLand.TopBoundary,
 ) = (:surface, :surface)
 """
-    boundary_var_types(::RichardsAtmosDrivenFluxBC{<:PrescribedPrecipitation,
+    boundary_var_types(::RichardsModel{FT}
+                        ::RichardsAtmosDrivenFluxBC{<:PrescribedPrecipitation,
                                                    <:Runoff.AbstractRunoffModel,
                                                   },
-                       ::ClimaLand.TopBoundary,
-                       ) where {FT}
+                        ::ClimaLand.TopBoundary,
+                        ) where {FT}
 
 An extension of the `boundary_var_types` method for RichardsAtmosDrivenFluxBC
 with no runoff modeled.
@@ -447,7 +451,7 @@ function boundary_flux(
 end
 
 """
-    ClimaLand.∂tendencyBC∂Y(
+    ClimaLand.set_dfluxBCdY!(
         model::RichardsModel,
         ::MoistureStateBC,
         boundary::ClimaLand.TopBoundary,
@@ -457,16 +461,16 @@ end
         t,
 )
 
-Computes and returns the derivative of the part of the
-implicit tendency in the top layer, due to the boundary
-condition, with respect to the state variable in the top layer.
+Computes the derivative of the flux in the top layer (due to the
+boundary condition), with respect to the state variable in the top layer.
+This value is then updated in-place in the cache.
 
-For a diffusion equation like Richards equation with a single state
-variable, this is given by
-`∂T_N∂Y_N = [-∂/∂z(∂F_bc/∂Y_N)]_N`, where `N` indicates the top
-layer cell index.
+For Richards equation (a diffusion equation with a single state variable),
+this is given by `∂F_bc/∂Y_N= -K_N (∂ψ_bc/∂ϑ_N) / Δz`,
+where `N` indicates the top layer cell index and
+`ψ_bc` is the pressure head at the boundary condition.
 """
-function ClimaLand.∂tendencyBC∂Y(
+function ClimaLand.set_dfluxBCdY!(
     model::RichardsModel,
     ::MoistureStateBC,
     boundary::ClimaLand.TopBoundary,
@@ -476,56 +480,44 @@ function ClimaLand.∂tendencyBC∂Y(
     t,
 )
     (; ν, hydrology_cm, S_s, θ_r) = model.parameters
-    fs = ClimaLand.Domains.obtain_face_space(model.domain.space.subsurface)
-    face_len = ClimaCore.Utilities.PlusHalf(ClimaCore.Spaces.nlevels(fs) - 1)
-    interpc2f_op = Operators.InterpolateC2F(
-        bottom = Operators.Extrapolate(),
-        top = Operators.Extrapolate(),
-    )
-    K = Fields.level(interpc2f_op.(p.soil.K), face_len)
-    dψ = Fields.level(
-        interpc2f_op.(dψdϑ.(hydrology_cm, Y.soil.ϑ_l, ν, θ_r, S_s)),
-        face_len,
-    )
-    return ClimaCore.Fields.FieldVector(;
-        :soil => (; :ϑ_l => @. -K / Δz * dψ / (2 * Δz)),
+
+    # Copy center variables to top face space
+    KN = Domains.top_center_to_surface(p.soil.K)
+    hydrology_cmN = Domains.top_center_to_surface(hydrology_cm)
+    ϑ_lN = Domains.top_center_to_surface(Y.soil.ϑ_l)
+    νN = Domains.top_center_to_surface(ν)
+    θ_rN = Domains.top_center_to_surface(θ_r)
+    S_sN = Domains.top_center_to_surface(S_s)
+
+
+    # Get the local geometry of the face space, then extract the top level
+    levels = ClimaCore.Spaces.nlevels(Domains.obtain_face_space(axes(p.soil.K)))
+    local_geometry_faceN = ClimaCore.Fields.level(
+        Fields.local_geometry_field(Domains.obtain_face_space(axes(p.soil.K))),
+        levels - ClimaCore.Utilities.half,
     )
+
+    # Update dfluxBCdY at the top boundary in place
+    # Calculate the value and convert it to a Covariant3Vector
+    @. p.soil.dfluxBCdY =
+        covariant3_unit_vector(local_geometry_faceN) *
+        (KN * dψdϑ(hydrology_cmN, ϑ_lN, νN, θ_rN, S_sN) / Δz)
+    return nothing
 end
 
 """
-    ClimaLand.∂tendencyBC∂Y(
-        ::RichardsModel,
-        ::AbstractWaterBC,
-        boundary::ClimaLand.TopBoundary,
-        Δz,
-        Y,
-        p,
-        t,
-)
-
-A default method which computes and returns the zero for the
-derivative of the part of the
-implicit tendency in the top layer, due to the boundary
-condition, with respect to the state variable in the top layer.
+    covariant3_unit_vector(local_geometry)
 
-For a diffusion equation like Richards equation with a single state
-variable, this is given by
-`∂T_N∂Y_N = [-∂/∂z(∂F_bc/∂Y_N)]_N`, where `N` indicates the top
-layer cell index.
+A function to compute the unit vector in the direction of the normal
+to the surface.
 
-If `F_bc` can be approximated as independent of `Y_N`, the derivative
-is zero.
+Adapted from ClimaAtmos.jl's unit_basis_vector_data function.
 """
-function ClimaLand.∂tendencyBC∂Y(
-    ::RichardsModel,
-    ::AbstractWaterBC,
-    boundary::ClimaLand.TopBoundary,
-    Δz,
-    Y,
-    p,
-    t,
-)
-    return ClimaCore.Fields.FieldVector(; :soil => (; :ϑ_l => zeros(axes(Δz))))
+function covariant3_unit_vector(local_geometry)
+    FT = Geometry.undertype(typeof(local_geometry))
+    data =
+        FT(1) / Geometry._norm(Geometry.Covariant3Vector(FT(1)), local_geometry)
+    return Geometry.Covariant3Vector(data)
 end
 
 # BC type for the soil heat equation
@@ -665,7 +657,7 @@ abstract type AbstractEnergyHydrologyBC <: ClimaLand.AbstractBC end
     WaterHeatBC{W <: AbstractWaterBC, H <: AbstractHeatBC} <:
        AbstractEnergyHydrologyBC
 
-A general struct used to store the boundary conditions for Richards 
+A general struct used to store the boundary conditions for Richards
 and the soil heat equations separately; useful when the boundary
 conditions for each component are independent of each other.
 """
@@ -731,7 +723,7 @@ end
 
 The list of domain names for additional variables added to the
 EnergyHydrology model auxiliary state, which defaults to adding storage for the
- boundary flux field. 
+ boundary flux field.
 
 Because we supply boundary conditions for water and heat, we found it convenient to
 have these stored as a NamedTuple under the names `top_bc` and `bottom_bc`.
@@ -808,6 +800,7 @@ boundary_var_domain_names(
 ) = (:surface, :surface, :surface, :surface, :surface, :surface)
 """
     boundary_var_types(
+        ::EnergyHydrology{FT},
         ::AtmosDrivenFluxBC{
             <:PrescribedAtmosphere{FT},
             <:AbstractRadiativeDrivers{FT},
@@ -890,7 +883,7 @@ end
                                     <:Runoff.TOPMODELRunoff,
                                     }, ::ClimaLand.TopBoundary)
 
-An extension of the `boundary_vars` method for AtmosDrivenFluxBC with 
+An extension of the `boundary_vars` method for AtmosDrivenFluxBC with
 TOPMODELRunoff. This
 adds the surface conditions (SHF, LHF, evaporation, and resistance) and the
 net radiation to the auxiliary variables.
@@ -949,6 +942,7 @@ boundary_var_domain_names(
 
 """
     boundary_var_types(
+        model::EnergyHydrology{FT},
         ::AtmosDrivenFluxBC{
             <:PrescribedAtmosphere{FT},
             <:AbstractRadiativeDrivers{FT},
@@ -979,3 +973,37 @@ boundary_var_types(
     FT,
     FT,
 )
+
+"""
+    boundary_vars(::MoistureStateBC, ::ClimaLand.TopBoundary)
+
+An extension of the `boundary_vars` method for MoistureStateBC at the
+top boundary.
+
+These variables are updated in place in `boundary_flux`.
+"""
+boundary_vars(bc::MoistureStateBC, ::ClimaLand.TopBoundary) =
+    (:top_bc, :dfluxBCdY)
+
+"""
+    boundary_var_domain_names(::MoistureStateBC, ::ClimaLand.TopBoundary)
+
+An extension of the `boundary_var_domain_names` method for MoistureStateBC at the
+top boundary.
+"""
+boundary_var_domain_names(bc::MoistureStateBC, ::ClimaLand.TopBoundary) =
+    (:surface, :surface)
+"""
+    boundary_var_types(::RichardsModel{FT},
+                        ::MoistureStateBC,
+                        ::ClimaLand.TopBoundary,
+                        ) where {FT}
+
+An extension of the `boundary_var_types` method for MoistureStateBC at the
+    top boundary.
+"""
+boundary_var_types(
+    model::RichardsModel{FT},
+    bc::MoistureStateBC,
+    ::ClimaLand.TopBoundary,
+) where {FT} = (FT, ClimaCore.Geometry.Covariant3Vector{FT})
diff --git a/src/standalone/Soil/rre.jl b/src/standalone/Soil/rre.jl
index 6f3d404252..ca77bef025 100644
--- a/src/standalone/Soil/rre.jl
+++ b/src/standalone/Soil/rre.jl
@@ -1,7 +1,11 @@
+using ClimaCore.MatrixFields
+import ClimaCore.MatrixFields: @name, ⋅
+import UnrolledUtilities
+
 """
     RichardsParameters{F <: Union{<: AbstractFloat, ClimaCore.Fields.Field}, C <: AbstractSoilHydrologyClosure}
 
-A struct for storing parameters of the `RichardModel`.
+A struct for storing parameters of the `RichardsModel`.
 $(DocStringExtensions.FIELDS)
 """
 struct RichardsParameters{
@@ -118,6 +122,17 @@ function make_update_boundary_fluxes(model::RichardsModel)
             p,
             t,
         )
+
+        # Update `p.soil.dfluxBCdY`
+        set_dfluxBCdY!(
+            model,
+            model.boundary_conditions.top,
+            TopBoundary(),
+            Δz_top,
+            Y,
+            p,
+            t,
+        )
     end
     return update_boundary_fluxes!
 end
@@ -330,84 +345,76 @@ end
 
 
 """
-    RichardsTridiagonalW{R, J, W, T} <: ClimaLand.AbstractTridiagonalW
+    ImplicitEquationJacobian{M, S}
 
 A struct containing the necessary information for constructing a tridiagonal
 Jacobian matrix (`W`) for solving Richards equation, treating only the vertical
 diffusion term implicitly.
 
+`matrix` is a block matrix containing the tri-diagonal matrix `∂ϑres∂ϑ`
+in the RichardsModel case.
+`solver` is a diagonal solver for the RichardsModel case because our matrix is
+block diagonal.
+
 Note that the diagonal, upper diagonal, and lower diagonal entry values
 are stored in this struct and updated in place.
-$(DocStringExtensions.FIELDS)
 """
-struct RichardsTridiagonalW{R, J, JA, T, A} <: ClimaLand.AbstractTridiagonalW
-    "Reference to dtγ, which is specified by the ODE solver"
-    dtγ_ref::R
-    "Diagonal entries of the Jacobian stored as a ClimaCore.Fields.Field"
-    ∂ϑₜ∂ϑ::J
-    "Array of tridiagonal matrices containing W for each column"
-    W_column_arrays::JA
-    "An allocated cache used to evaluate ldiv!"
-    temp1::T
-    "An allocated cache used to evaluate ldiv!"
-    temp2::T
-    "A flag indicating whether this struct is used to compute Wfact_t or Wfact"
-    transform::Bool
-    "A pre-allocated cache storing ones on the face space"
-    ones_face_space::A
+struct ImplicitEquationJacobian{M, S}
+    "Jacobian matrix stored as a MatrixFields.FieldMatrix"
+    matrix::M
+    "Solver to use for solving the tridiagonal system"
+    solver::S
 end
 
 """
-    RichardsTridiagonalW(
+    ImplicitEquationJacobian(
         Y::ClimaCore.Fields.FieldVector;
-        transform::Bool = false
 )
 
-Outer constructor for the RichardsTridiagonalW Jacobian
+Outer constructor for the ImplicitEquationJacobian Jacobian
 matrix struct.
 
 Initializes all variables to zeros.
 """
-function RichardsTridiagonalW(
-    Y::ClimaCore.Fields.FieldVector;
-    transform::Bool = false,
-)
-    FT = eltype(Y.soil.ϑ_l)
+function ImplicitEquationJacobian(Y::ClimaCore.Fields.FieldVector)
+    FT = eltype(Y)
     center_space = axes(Y.soil.ϑ_l)
-    N = Spaces.nlevels(center_space)
-
-    tridiag_type = Operators.StencilCoefs{-1, 1, NTuple{3, FT}}
-    ∂ϑₜ∂ϑ = Fields.Field(tridiag_type, center_space)
-
-    fill!(parent(∂ϑₜ∂ϑ), FT(0))
-
-    W_column_arrays = [
-        LinearAlgebra.Tridiagonal(
-            zeros(FT, N - 1) .+ FT(0),
-            zeros(FT, N) .+ FT(0),
-            zeros(FT, N - 1) .+ FT(0),
-        ) for _ in 1:Threads.nthreads()
-    ]
-    dtγ_ref = Ref(FT(0))
-    temp1 = similar(Y.soil.ϑ_l)
-    temp1 .= FT(0)
-    temp2 = similar(Y.soil.ϑ_l)
-    temp2 .= FT(0)
-
-    face_space = ClimaLand.Domains.obtain_face_space(center_space)
-    ones_face_space = ones(face_space)
-
-    return RichardsTridiagonalW(
-        dtγ_ref,
-        ∂ϑₜ∂ϑ,
-        W_column_arrays,
-        temp1,
-        temp2,
-        transform,
-        ones_face_space,
+
+    # Cosntruct a tridiagonal matrix that will be used as the Jacobian
+    tridiag_type = MatrixFields.TridiagonalMatrixRow{FT}
+    # Create a field containing a `TridiagonalMatrixRow` at each point
+    tridiag_field = Fields.Field(tridiag_type, center_space)
+    fill!(parent(tridiag_field), NaN)
+
+    # Get all prognostic vars in soil, and create a tridiagonal matrix for each
+    soil_varnames = MatrixFields.top_level_names(Y.soil)
+    varnames = map(
+        name -> MatrixFields.append_internal_name(@name(soil), name),
+        soil_varnames,
     )
+    matrix = MatrixFields.FieldMatrix(
+        UnrolledUtilities.unrolled_map(
+            x -> (x, x) => copy(tridiag_field),
+            varnames,
+        )...,
+    )
+
+    # Set up block diagonal solver for block Jacobian
+    alg = MatrixFields.BlockDiagonalSolve()
+    solver = MatrixFields.FieldMatrixSolver(alg, matrix, Y)
+
+    return ImplicitEquationJacobian(matrix, solver)
 end
 
+Base.similar(w::ImplicitEquationJacobian) = w
+
+function LinearAlgebra.ldiv!(
+    x::Fields.FieldVector,
+    A::ImplicitEquationJacobian,
+    b::Fields.FieldVector,
+)
+    MatrixFields.field_matrix_solve!(A.solver, x, A.matrix, b)
+end
 
 """
     ClimaLand.make_update_jacobian(model::RichardsModel{FT}) where {FT}
@@ -418,35 +425,69 @@ Using this Jacobian with a backwards Euler timestepper is equivalent
 to using the modified Picard scheme of Celia et al. (1990).
 """
 function ClimaLand.make_update_jacobian(model::RichardsModel{FT}) where {FT}
-    update_aux! = make_update_aux(model)
-    update_boundary_fluxes! = make_update_boundary_fluxes(model)
-    function update_jacobian!(W::RichardsTridiagonalW, Y, p, dtγ, t)
-        (; dtγ_ref, ∂ϑₜ∂ϑ) = W
-        dtγ_ref[] = dtγ
+    function update_jacobian!(jacobian::ImplicitEquationJacobian, Y, p, dtγ, t)
+        (; matrix) = jacobian
         (; ν, hydrology_cm, S_s, θ_r) = model.parameters
-        update_aux!(p, Y, t)
-        update_boundary_fluxes!(p, Y, t)
-        divf2c_op = Operators.DivergenceF2C(
-            top = Operators.SetValue(Geometry.WVector.(FT(0))),
-            bottom = Operators.SetValue(Geometry.WVector.(FT(0))),
-        )
-        divf2c_stencil = Operators.Operator2Stencil(divf2c_op)
+
+        # Create divergence operator
+        divf2c_op = Operators.DivergenceF2C()
+        divf2c_matrix = MatrixFields.operator_matrix(divf2c_op)
+        # Create gradient operator, and set gradient at boundaries to 0
         gradc2f_op = Operators.GradientC2F(
-            top = Operators.SetGradient(Geometry.WVector.(FT(0))),
-            bottom = Operators.SetGradient(Geometry.WVector.(FT(0))),
+            top = Operators.SetGradient(Geometry.WVector(FT(0))),
+            bottom = Operators.SetGradient(Geometry.WVector(FT(0))),
         )
-        gradc2f_stencil = Operators.Operator2Stencil(gradc2f_op)
+        gradc2f_matrix = MatrixFields.operator_matrix(gradc2f_op)
+        # Create interpolation operator
         interpc2f_op = Operators.InterpolateC2F(
             bottom = Operators.Extrapolate(),
             top = Operators.Extrapolate(),
         )
-        compose = Operators.ComposeStencils()
 
-        @. ∂ϑₜ∂ϑ = compose(
-            divf2c_stencil(Geometry.Covariant3Vector(W.ones_face_space)),
-            (
-                interpc2f_op(p.soil.K) * ClimaLand.to_scalar_coefs(
-                    gradc2f_stencil(
+        # The derivative of the residual with respect to the prognostic variable
+        ∂ϑres∂ϑ = matrix[@name(soil.ϑ_l), @name(soil.ϑ_l)]
+
+        # If the top BC is a `MoistureStateBC`, add the term from the top BC
+        #  flux before applying divergence
+        if haskey(p.soil, :dfluxBCdY)
+            dfluxBCdY = p.soil.dfluxBCdY
+            topBC_op = Operators.SetBoundaryOperator(
+                top = Operators.SetValue(dfluxBCdY),
+                bottom = Operators.SetValue(
+                    Geometry.Covariant3Vector(zero(FT)),
+                ),
+            )
+            # Add term from top boundary condition before applying divergence
+            # Note: need to pass 3D field on faces to `topBC_op`. Interpolating `K` to faces
+            #  for this is inefficient - we should find a better solution.
+            @. ∂ϑres∂ϑ =
+                -dtγ * (
+                    divf2c_matrix() ⋅ (
+                        MatrixFields.DiagonalMatrixRow(
+                            interpc2f_op(-p.soil.K),
+                        ) ⋅ gradc2f_matrix() ⋅ MatrixFields.DiagonalMatrixRow(
+                            ClimaLand.Soil.dψdϑ(
+                                hydrology_cm,
+                                Y.soil.ϑ_l,
+                                ν,
+                                θ_r,
+                                S_s,
+                            ),
+                        ) + MatrixFields.LowerDiagonalMatrixRow(
+                            topBC_op(
+                                Geometry.Covariant3Vector(
+                                    zero(interpc2f_op(p.soil.K)),
+                                ),
+                            ),
+                        )
+                    )
+                ) - (I,)
+        else
+            @. ∂ϑres∂ϑ =
+                -dtγ * (
+                    divf2c_matrix() ⋅
+                    MatrixFields.DiagonalMatrixRow(interpc2f_op(-p.soil.K)) ⋅
+                    gradc2f_matrix() ⋅ MatrixFields.DiagonalMatrixRow(
                         ClimaLand.Soil.dψdϑ(
                             hydrology_cm,
                             Y.soil.ϑ_l,
@@ -454,29 +495,9 @@ function ClimaLand.make_update_jacobian(model::RichardsModel{FT}) where {FT}
                             θ_r,
                             S_s,
                         ),
-                    ),
-                )
-            ),
-        )
-        # Hardcoded for single column: FIX!
-        # The derivative of the tendency may eventually live in boundary vars and be updated there. but TBD
-        z = Fields.coordinate_field(axes(Y.soil.ϑ_l)).z
-        Δz_top, Δz_bottom = ClimaLand.Domains.get_Δz(z)
-        ∂T_bc∂YN = ClimaLand.∂tendencyBC∂Y(
-            model,
-            model.boundary_conditions.top,
-            ClimaLand.TopBoundary(),
-            Δz_top,
-            Y,
-            p,
-            t,
-        )
-        #TODO: allocate space in W? See how final implementation of stencils with boundaries works out
-        N = ClimaCore.Spaces.nlevels(axes(Y.soil.ϑ_l))
-        parent(ClimaCore.Fields.level(∂ϑₜ∂ϑ.coefs.:2, N)) .=
-            parent(ClimaCore.Fields.level(∂ϑₜ∂ϑ.coefs.:2, N)) .+
-            parent(∂T_bc∂YN.soil.ϑ_l)
-
+                    )
+                ) - (I,)
+        end
     end
     return update_jacobian!
 end
@@ -504,7 +525,7 @@ end
 """
     is_saturated(Y, model::RichardsModel)
 
-A helper function which can be used to indicate whether a layer of soil is 
+A helper function which can be used to indicate whether a layer of soil is
 saturated.
 """
 function is_saturated(Y, model::RichardsModel)
diff --git a/test/shared_utilities/implicit_timestepping/richards_model.jl b/test/shared_utilities/implicit_timestepping/richards_model.jl
index c7b469053a..61fa79f832 100644
--- a/test/shared_utilities/implicit_timestepping/richards_model.jl
+++ b/test/shared_utilities/implicit_timestepping/richards_model.jl
@@ -2,7 +2,7 @@ using Test
 import ClimaComms
 @static pkgversion(ClimaComms) >= v"0.6" && ClimaComms.@import_required_backends
 using LinearAlgebra
-using ClimaCore
+import ClimaCore: MatrixFields
 import ClimaParams as CP
 using ClimaLand
 using ClimaLand.Domains: Column, HybridBox
@@ -53,10 +53,10 @@ for FT in (Float32, Float64)
             # We do not set the initial aux state here because
             # we want to test that it is updated correctly in
             # the jacobian correctly.
-            W = RichardsTridiagonalW(Y)
-            Wfact! = make_tendency_jacobian(soil)
+            jacobian = ImplicitEquationJacobian(Y)
+            jac_tendency! = make_tendency_jacobian(soil)
             dtγ = FT(1.0)
-            Wfact!(W, Y, p, dtγ, FT(0.0))
+            jac_tendency!(jacobian, Y, p, dtγ, FT(0.0))
 
             K_ic = hydraulic_conductivity(
                 hcm,
@@ -66,57 +66,74 @@ for FT in (Float32, Float64)
             dz = FT(0.01)
             dψdϑ_ic = dψdϑ(hcm, FT(0.24), ν, θ_r, S_s)
 
-            @test all(parent(W.temp2) .== FT(0.0))
-            @test all(parent(W.temp2) .== FT(0.0))
-            @test W.transform == false
-            @test typeof(W.W_column_arrays) <:
-                  Vector{LinearAlgebra.Tridiagonal{FT, Vector{FT}}}
-            @test length(W.W_column_arrays) == 1
+            @test jacobian.solver isa MatrixFields.FieldMatrixSolver
+            @test jacobian.solver.alg isa MatrixFields.BlockDiagonalSolve
+            @test jacobian.matrix.keys.values ==
+                  ((MatrixFields.@name(soil.ϑ_l), MatrixFields.@name(soil.ϑ_l)),)
+
+            jac_ϑ_l = jacobian.matrix[
+                MatrixFields.@name(soil.ϑ_l),
+                MatrixFields.@name(soil.ϑ_l)
+            ]
             if typeof(domain) <: Column
+                # Check diagonals on either side of the main diagonal
                 @test all(
-                    parent(W.∂ϑₜ∂ϑ.coefs.:1)[2:end] .≈
-                    parent(W.∂ϑₜ∂ϑ.coefs.:3)[1:(end - 1)],
+                    parent(jac_ϑ_l.entries.:1)[2:end] .≈
+                    parent(jac_ϑ_l.entries.:3)[1:(end - 1)],
                 )
-                @test Array(parent(W.∂ϑₜ∂ϑ.coefs.:1))[1] == FT(0)
-                @test Array(parent(W.∂ϑₜ∂ϑ.coefs.:3))[end] == FT(0)
+                @test Array(parent(jac_ϑ_l.entries.:1))[1] == FT(0)
+                @test Array(parent(jac_ϑ_l.entries.:3))[end] == FT(0)
                 @test all(
-                    Array(parent(W.∂ϑₜ∂ϑ.coefs.:1))[2:end] .≈
-                    K_ic / dz^2 * dψdϑ_ic,
+                    Array(parent(jac_ϑ_l.entries.:1))[2:end] .≈
+                    dtγ * (K_ic / dz^2 * dψdϑ_ic),
                 )
-                @test Array(parent(W.∂ϑₜ∂ϑ.coefs.:2))[1] .≈
-                      -K_ic / dz^2 * dψdϑ_ic
+                # Check values on main diagonal: note jac_ϑ_l is I-γJ
+                @test Array(parent(jac_ϑ_l.entries.:2))[1] .≈
+                      dtγ * (-K_ic / dz^2 * dψdϑ_ic) - I
                 @test all(
-                    Array(parent(W.∂ϑₜ∂ϑ.coefs.:2))[2:(end - 1)] .≈
-                    -2 * K_ic / dz^2 * dψdϑ_ic,
+                    Array(parent(jac_ϑ_l.entries.:2))[2:(end - 1)] .≈
+                    dtγ * (-2 * K_ic / dz^2 * dψdϑ_ic) - I,
                 )
-                @test Array(parent(W.∂ϑₜ∂ϑ.coefs.:2))[end] .≈
-                      -K_ic / dz^2 * dψdϑ_ic - K_ic / (dz * dz / 2) * dψdϑ_ic
+                @test Array(parent(jac_ϑ_l.entries.:2))[end] .≈
+                      dtγ * (
+                    -K_ic / dz^2 * dψdϑ_ic - K_ic / (dz * dz / 2) * dψdϑ_ic
+                ) - I
             elseif typeof(domain) <: HybridBox
+                # Check diagonals on either side of the main diagonal
                 @test all(
-                    Array(parent(W.∂ϑₜ∂ϑ.coefs.:1))[2:end, :, 1, 1, 1] .≈
-                    Array(parent(W.∂ϑₜ∂ϑ.coefs.:3))[1:(end - 1), :, 1, 1, 1],
+                    Array(parent(jac_ϑ_l.entries.:1))[2:end, :, 1, 1, 1] .≈
+                    Array(parent(jac_ϑ_l.entries.:3))[1:(end - 1), :, 1, 1, 1],
                 )
                 @test all(
-                    Array(parent(W.∂ϑₜ∂ϑ.coefs.:1))[1, :, 1, 1, 1] .== FT(0),
+                    Array(parent(jac_ϑ_l.entries.:1))[1, :, 1, 1, 1] .== FT(0),
                 )
                 @test all(
-                    Array(parent(W.∂ϑₜ∂ϑ.coefs.:3))[end, :, 1, 1, 1] .== FT(0),
+                    Array(parent(jac_ϑ_l.entries.:3))[end, :, 1, 1, 1] .==
+                    FT(0),
                 )
                 @test all(
-                    Array(parent(W.∂ϑₜ∂ϑ.coefs.:1))[2:end, :, 1, 1, 1] .≈
-                    K_ic / dz^2 * dψdϑ_ic,
+                    Array(parent(jac_ϑ_l.entries.:1))[2:end, :, 1, 1, 1] .≈
+                    dtγ * (K_ic / dz^2 * dψdϑ_ic),
                 )
+                # Check values on main diagonal: note jac_ϑ_l is I-γJ
                 @test all(
-                    Array(parent(W.∂ϑₜ∂ϑ.coefs.:2))[1, :, 1, 1, 1] .≈
-                    -K_ic / dz^2 * dψdϑ_ic,
+                    Array(parent(jac_ϑ_l.entries.:2))[1, :, 1, 1, 1] .≈
+                    dtγ * (-K_ic / dz^2 * dψdϑ_ic) - I,
                 )
                 @test all(
-                    Array(parent(W.∂ϑₜ∂ϑ.coefs.:2))[2:(end - 1), :, 1, 1, 1] .≈
-                    -2 * K_ic / dz^2 * dψdϑ_ic,
+                    Array(parent(jac_ϑ_l.entries.:2))[
+                        2:(end - 1),
+                        :,
+                        1,
+                        1,
+                        1,
+                    ] .≈ dtγ * (-2 * K_ic / dz^2 * dψdϑ_ic) - I,
                 )
                 @test all(
-                    Array(parent(W.∂ϑₜ∂ϑ.coefs.:2))[end, :, 1, 1, 1] .≈
-                    -K_ic / dz^2 * dψdϑ_ic - K_ic / (dz * dz / 2) * dψdϑ_ic,
+                    Array(parent(jac_ϑ_l.entries.:2))[end, :, 1, 1, 1] .≈
+                    dtγ *
+                    (-K_ic / dz^2 * dψdϑ_ic - K_ic / (dz * dz / 2) * dψdϑ_ic) -
+                    I,
                 )
             end
 
@@ -162,10 +179,10 @@ for FT in (Float32, Float64)
             # We do not set the initial aux state here because
             # we want to test that it is updated correctly in
             # the jacobian correctly.
-            W = RichardsTridiagonalW(Y)
-            Wfact! = make_tendency_jacobian(soil)
+            jacobian = ImplicitEquationJacobian(Y)
+            jacobian_tendency! = make_tendency_jacobian(soil)
             dtγ = FT(1.0)
-            Wfact!(W, Y, p, dtγ, FT(0.0))
+            jacobian_tendency!(jacobian, Y, p, dtγ, FT(0.0))
 
             K_ic = hydraulic_conductivity(
                 hcm,
@@ -174,27 +191,36 @@ for FT in (Float32, Float64)
             )
             dz = FT(0.01)
             dψdϑ_ic = dψdϑ(hcm, FT(0.24), ν, θ_r, S_s)
+            jac_ϑ_l = jacobian.matrix[
+                MatrixFields.@name(soil.ϑ_l),
+                MatrixFields.@name(soil.ϑ_l)
+            ]
             if typeof(domain) <: Column
-                @test Array(parent(W.∂ϑₜ∂ϑ.coefs.:2))[1] .≈
-                      -K_ic / dz^2 * dψdϑ_ic
+                @test Array(parent(jac_ϑ_l.entries.:2))[1] .≈
+                      dtγ * (-K_ic / dz^2 * dψdϑ_ic) - I
                 @test all(
-                    Array(parent(W.∂ϑₜ∂ϑ.coefs.:2))[2:(end - 1)] .≈
-                    -2 * K_ic / dz^2 * dψdϑ_ic,
+                    Array(parent(jac_ϑ_l.entries.:2))[2:(end - 1)] .≈
+                    dtγ * (-2 * K_ic / dz^2 * dψdϑ_ic) - I,
                 )
-                @test Array(parent(W.∂ϑₜ∂ϑ.coefs.:2))[end] .≈
-                      -K_ic / dz^2 * dψdϑ_ic
+                @test Array(parent(jac_ϑ_l.entries.:2))[end] .≈
+                      dtγ * (-K_ic / dz^2 * dψdϑ_ic) - I
             elseif typeof(domain) <: HybridBox
                 @test all(
-                    Array(parent(W.∂ϑₜ∂ϑ.coefs.:2))[1, :, 1, 1, 1] .≈
-                    -K_ic / dz^2 * dψdϑ_ic,
+                    Array(parent(jac_ϑ_l.entries.:2))[1, :, 1, 1, 1] .≈
+                    dtγ * (-K_ic / dz^2 * dψdϑ_ic) - I,
                 )
                 @test all(
-                    Array(parent(W.∂ϑₜ∂ϑ.coefs.:2))[2:(end - 1), :, 1, 1, 1] .≈
-                    -2 * K_ic / dz^2 * dψdϑ_ic,
+                    Array(parent(jac_ϑ_l.entries.:2))[
+                        2:(end - 1),
+                        :,
+                        1,
+                        1,
+                        1,
+                    ] .≈ dtγ * (-2 * K_ic / dz^2 * dψdϑ_ic) - I,
                 )
                 @test all(
-                    Array(parent(W.∂ϑₜ∂ϑ.coefs.:2))[end, :, 1, 1, 1] .≈
-                    -K_ic / dz^2 * dψdϑ_ic,
+                    Array(parent(jac_ϑ_l.entries.:2))[end, :, 1, 1, 1] .≈
+                    dtγ * (-K_ic / dz^2 * dψdϑ_ic) - I,
                 )
             end
         end