diff --git a/README.md b/README.md
index 59ad374..108bb4d 100644
--- a/README.md
+++ b/README.md
@@ -60,6 +60,7 @@ S. Gupta et al 2022, 2024](https://github.com/CliMA/ClimaArtifacts/tree/main/soi
 - [ERA5 Monthly Averages on Single Levels 1979-2024](https:////github.com/CliMA/ClimaArtifacts/tree/main/era5_monthly_averages_single_level_1979_2024)
 - [Lehmann et al. 2008 Fig. 8 evaporation data](https:////github.com/CliMA/ClimaArtifacts/tree/main/lehmann2008_evaporation)
 - [Forty years of ERA5 forcing data for ClimaLand](https:////github.com/CliMA/ClimaArtifacts/tree/main/forty_yrs_era5_land_forcing_data)
+- [Water conservation test data for RichardsModel](https:////github.com/CliMA/ClimaArtifacts/tree/main/water_conservation_test)
 
 ### Coupler/shared
 
diff --git a/water_conservation_test/Manifest.toml b/water_conservation_test/Manifest.toml
new file mode 100644
index 0000000..5ac0b9c
--- /dev/null
+++ b/water_conservation_test/Manifest.toml
@@ -0,0 +1,305 @@
+# This file is machine-generated - editing it directly is not advised
+
+julia_version = "1.11.2"
+manifest_format = "2.0"
+project_hash = "74753091693fecdd4a33be59e19da0ff479810d2"
+
+[[deps.ArgTools]]
+uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"
+version = "1.1.2"
+
+[[deps.ArtifactUtils]]
+deps = ["Downloads", "Git", "HTTP", "Pkg", "ProgressLogging", "SHA", "TOML", "gh_cli_jll"]
+git-tree-sha1 = "f7c3ca6c70ea6b035effe5428eb05231ba890c2d"
+uuid = "8b73e784-e7d8-4ea5-973d-377fed4e3bce"
+version = "0.2.4"
+
+[[deps.Artifacts]]
+uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
+version = "1.11.0"
+
+[[deps.Base64]]
+uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
+version = "1.11.0"
+
+[[deps.BitFlags]]
+git-tree-sha1 = "0691e34b3bb8be9307330f88d1a3c3f25466c24d"
+uuid = "d1d4a3ce-64b1-5f1a-9ba4-7e7e69966f35"
+version = "0.1.9"
+
+[[deps.ClimaArtifactsHelper]]
+deps = ["ArtifactUtils", "Downloads", "Pkg", "REPL", "SHA"]
+path = "../ClimaArtifactsHelper.jl"
+uuid = "6ffa2572-8378-4377-82eb-ea11db28b255"
+version = "0.0.1"
+
+[[deps.CodecZlib]]
+deps = ["TranscodingStreams", "Zlib_jll"]
+git-tree-sha1 = "bce6804e5e6044c6daab27bb533d1295e4a2e759"
+uuid = "944b1d66-785c-5afd-91f1-9de20f533193"
+version = "0.7.6"
+
+[[deps.ConcurrentUtilities]]
+deps = ["Serialization", "Sockets"]
+git-tree-sha1 = "ea32b83ca4fefa1768dc84e504cc0a94fb1ab8d1"
+uuid = "f0e56b4a-5159-44fe-b623-3e5288b988bb"
+version = "2.4.2"
+
+[[deps.Dates]]
+deps = ["Printf"]
+uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
+version = "1.11.0"
+
+[[deps.Downloads]]
+deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"]
+uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
+version = "1.6.0"
+
+[[deps.ExceptionUnwrapping]]
+deps = ["Test"]
+git-tree-sha1 = "d36f682e590a83d63d1c7dbd287573764682d12a"
+uuid = "460bff9d-24e4-43bc-9d9f-a8973cb893f4"
+version = "0.1.11"
+
+[[deps.Expat_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "e51db81749b0777b2147fbe7b783ee79045b8e99"
+uuid = "2e619515-83b5-522b-bb60-26c02a35a201"
+version = "2.6.4+1"
+
+[[deps.FileWatching]]
+uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"
+version = "1.11.0"
+
+[[deps.Git]]
+deps = ["Git_jll"]
+git-tree-sha1 = "04eff47b1354d702c3a85e8ab23d539bb7d5957e"
+uuid = "d7ba0133-e1db-5d97-8f8c-041e4b3a1eb2"
+version = "1.3.1"
+
+[[deps.Git_jll]]
+deps = ["Artifacts", "Expat_jll", "JLLWrappers", "LibCURL_jll", "Libdl", "Libiconv_jll", "OpenSSL_jll", "PCRE2_jll", "Zlib_jll"]
+git-tree-sha1 = "399f4a308c804b446ae4c91eeafadb2fe2c54ff9"
+uuid = "f8c6e375-362e-5223-8a59-34ff63f689eb"
+version = "2.47.1+0"
+
+[[deps.HTTP]]
+deps = ["Base64", "CodecZlib", "ConcurrentUtilities", "Dates", "ExceptionUnwrapping", "Logging", "LoggingExtras", "MbedTLS", "NetworkOptions", "OpenSSL", "PrecompileTools", "Random", "SimpleBufferStream", "Sockets", "URIs", "UUIDs"]
+git-tree-sha1 = "6c22309e9a356ac1ebc5c8a217045f9bae6f8d9a"
+uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3"
+version = "1.10.13"
+
+[[deps.InteractiveUtils]]
+deps = ["Markdown"]
+uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
+version = "1.11.0"
+
+[[deps.JLLWrappers]]
+deps = ["Artifacts", "Preferences"]
+git-tree-sha1 = "be3dc50a92e5a386872a493a10050136d4703f9b"
+uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
+version = "1.6.1"
+
+[[deps.LibCURL]]
+deps = ["LibCURL_jll", "MozillaCACerts_jll"]
+uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"
+version = "0.6.4"
+
+[[deps.LibCURL_jll]]
+deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
+uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"
+version = "8.6.0+0"
+
+[[deps.LibGit2]]
+deps = ["Base64", "LibGit2_jll", "NetworkOptions", "Printf", "SHA"]
+uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"
+version = "1.11.0"
+
+[[deps.LibGit2_jll]]
+deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll"]
+uuid = "e37daf67-58a4-590a-8e99-b0245dd2ffc5"
+version = "1.7.2+0"
+
+[[deps.LibSSH2_jll]]
+deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
+uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"
+version = "1.11.0+1"
+
+[[deps.Libdl]]
+uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
+version = "1.11.0"
+
+[[deps.Libiconv_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "61dfdba58e585066d8bce214c5a51eaa0539f269"
+uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531"
+version = "1.17.0+1"
+
+[[deps.Logging]]
+uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
+version = "1.11.0"
+
+[[deps.LoggingExtras]]
+deps = ["Dates", "Logging"]
+git-tree-sha1 = "f02b56007b064fbfddb4c9cd60161b6dd0f40df3"
+uuid = "e6f89c97-d47a-5376-807f-9c37f3926c36"
+version = "1.1.0"
+
+[[deps.Markdown]]
+deps = ["Base64"]
+uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
+version = "1.11.0"
+
+[[deps.MbedTLS]]
+deps = ["Dates", "MbedTLS_jll", "MozillaCACerts_jll", "NetworkOptions", "Random", "Sockets"]
+git-tree-sha1 = "c067a280ddc25f196b5e7df3877c6b226d390aaf"
+uuid = "739be429-bea8-5141-9913-cc70e7f3736d"
+version = "1.1.9"
+
+[[deps.MbedTLS_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"
+version = "2.28.6+0"
+
+[[deps.MozillaCACerts_jll]]
+uuid = "14a3606d-f60d-562e-9121-12d972cd8159"
+version = "2023.12.12"
+
+[[deps.NetworkOptions]]
+uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"
+version = "1.2.0"
+
+[[deps.OpenSSL]]
+deps = ["BitFlags", "Dates", "MozillaCACerts_jll", "OpenSSL_jll", "Sockets"]
+git-tree-sha1 = "38cb508d080d21dc1128f7fb04f20387ed4c0af4"
+uuid = "4d8831e6-92b7-49fb-bdf8-b643e874388c"
+version = "1.4.3"
+
+[[deps.OpenSSL_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "7493f61f55a6cce7325f197443aa80d32554ba10"
+uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95"
+version = "3.0.15+1"
+
+[[deps.PCRE2_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "efcefdf7-47ab-520b-bdef-62a2eaa19f15"
+version = "10.42.0+1"
+
+[[deps.Pkg]]
+deps = ["Artifacts", "Dates", "Downloads", "FileWatching", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "Random", "SHA", "TOML", "Tar", "UUIDs", "p7zip_jll"]
+uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
+version = "1.11.0"
+weakdeps = ["REPL"]
+
+    [deps.Pkg.extensions]
+    REPLExt = "REPL"
+
+[[deps.PrecompileTools]]
+deps = ["Preferences"]
+git-tree-sha1 = "5aa36f7049a63a1528fe8f7c3f2113413ffd4e1f"
+uuid = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
+version = "1.2.1"
+
+[[deps.Preferences]]
+deps = ["TOML"]
+git-tree-sha1 = "9306f6085165d270f7e3db02af26a400d580f5c6"
+uuid = "21216c6a-2e73-6563-6e65-726566657250"
+version = "1.4.3"
+
+[[deps.Printf]]
+deps = ["Unicode"]
+uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+version = "1.11.0"
+
+[[deps.ProgressLogging]]
+deps = ["Logging", "SHA", "UUIDs"]
+git-tree-sha1 = "80d919dee55b9c50e8d9e2da5eeafff3fe58b539"
+uuid = "33c8b6b6-d38a-422a-b730-caa89a2f386c"
+version = "0.1.4"
+
+[[deps.REPL]]
+deps = ["InteractiveUtils", "Markdown", "Sockets", "StyledStrings", "Unicode"]
+uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
+version = "1.11.0"
+
+[[deps.Random]]
+deps = ["SHA"]
+uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+version = "1.11.0"
+
+[[deps.SHA]]
+uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
+version = "0.7.0"
+
+[[deps.Serialization]]
+uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
+version = "1.11.0"
+
+[[deps.SimpleBufferStream]]
+git-tree-sha1 = "f305871d2f381d21527c770d4788c06c097c9bc1"
+uuid = "777ac1f9-54b0-4bf8-805c-2214025038e7"
+version = "1.2.0"
+
+[[deps.Sockets]]
+uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
+version = "1.11.0"
+
+[[deps.StyledStrings]]
+uuid = "f489334b-da3d-4c2e-b8f0-e476e12c162b"
+version = "1.11.0"
+
+[[deps.TOML]]
+deps = ["Dates"]
+uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
+version = "1.0.3"
+
+[[deps.Tar]]
+deps = ["ArgTools", "SHA"]
+uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"
+version = "1.10.0"
+
+[[deps.Test]]
+deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
+uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+version = "1.11.0"
+
+[[deps.TranscodingStreams]]
+git-tree-sha1 = "0c45878dcfdcfa8480052b6ab162cdd138781742"
+uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
+version = "0.11.3"
+
+[[deps.URIs]]
+git-tree-sha1 = "67db6cc7b3821e19ebe75791a9dd19c9b1188f2b"
+uuid = "5c2747f8-b7ea-4ff2-ba2e-563bfd36b1d4"
+version = "1.5.1"
+
+[[deps.UUIDs]]
+deps = ["Random", "SHA"]
+uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
+version = "1.11.0"
+
+[[deps.Unicode]]
+uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
+version = "1.11.0"
+
+[[deps.Zlib_jll]]
+deps = ["Libdl"]
+uuid = "83775a58-1f1d-513f-b197-d71354ab007a"
+version = "1.2.13+1"
+
+[[deps.gh_cli_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "b54ff233a6a183ad366d29fa0e3ff9bb921692db"
+uuid = "5d31d589-30fb-542f-b82d-10325e863e38"
+version = "2.35.0+0"
+
+[[deps.nghttp2_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"
+version = "1.59.0+0"
+
+[[deps.p7zip_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
+version = "17.4.0+2"
diff --git a/water_conservation_test/Project.toml b/water_conservation_test/Project.toml
new file mode 100644
index 0000000..8e4179f
--- /dev/null
+++ b/water_conservation_test/Project.toml
@@ -0,0 +1,2 @@
+[deps]
+ClimaArtifactsHelper = "6ffa2572-8378-4377-82eb-ea11db28b255"
diff --git a/water_conservation_test/README.md b/water_conservation_test/README.md
new file mode 100644
index 0000000..91f7c83
--- /dev/null
+++ b/water_conservation_test/README.md
@@ -0,0 +1,12 @@
+# Water conservation test dataset
+
+This dataset is used as a reference solution for a particular setup of the
+ClimaLand.jl soil water model (`RichardsModel`). It was created by running
+the `water_conservation.jl` experiment with an explicit solver and a very
+small timestep. To recreate it, we need to check out a commit of ClimaLand.jl
+where this type of model was stepped fully explicitly; it has since been changed
+to be stepped with a mixed implicit/explicit algorithm.
+
+Note that if the physics implemented in `RichardsModel` changes, this dataset
+will need to be regenerated so it can still be used as a comparison for the
+current code.
diff --git a/water_conservation_test/create_artifact.jl b/water_conservation_test/create_artifact.jl
new file mode 100644
index 0000000..f30981b
--- /dev/null
+++ b/water_conservation_test/create_artifact.jl
@@ -0,0 +1,20 @@
+using ClimaArtifactsHelper
+
+# Set the output directory
+output_dir = basename(@__DIR__) * "_artifact"
+if isdir(output_dir)
+    @warn "$output_dir already exists. Content will end up in the artifact and may be overwritten."
+    @warn "Abort this calculation, unless you know what you are doing."
+else
+    mkdir(output_dir)
+end
+
+# Run the bash script to create the artifact
+run(`bash run_water_conservation_exp.sh $output_dir`)
+
+@info "Data file generated!"
+
+# create_artifact_guided(
+#     output_dir;
+#     artifact_name = basename(@__DIR__),
+# )
diff --git a/water_conservation_test/run_water_conservation_exp.sh b/water_conservation_test/run_water_conservation_exp.sh
new file mode 100644
index 0000000..1cb3ddd
--- /dev/null
+++ b/water_conservation_test/run_water_conservation_exp.sh
@@ -0,0 +1,26 @@
+#!/bin/bash
+
+# This script checks out a specific commit of the ClimaLand.jl GitHub
+# repository, and uses it to run the provided `water_conservation_exp.jl`
+# experiment. That file will run the experiments with both flux and Dirichlet
+# boundary conditions, and save the solution results.
+# The output data files are saved in the paths water_conserv_exp/ref_soln_flux.csv
+# and water_conserv_exp/ref_soln_dirichlet.csv, relative to the directory
+# where this script is run.
+
+# Inform the user that the script requires one input argument, if it is not provided
+[[ $# -lt 1 ]] && echo "Error! Usage: $0 <output_dir>" && exit 0
+
+# Parse the input argument for the output directory
+output_dir=$1
+mkdir -p $output_dir
+
+# Clone the ClimaLand.jl repository
+git clone https://github.com/CliMA/ClimaLand.jl.git
+cd ClimaLand.jl
+
+# Checkout the commit that we know has what we need
+git -c advice.detachedHEAD=false checkout d7cc27a9b33227cfc1aa3745a94f3ffcc595fa5c
+
+julia --project=experiments -e 'using Pkg; Pkg.instantiate()'
+julia --project=experiments ../water_conservation_exp.jl $output_dir
diff --git a/water_conservation_test/water_conservation_exp.jl b/water_conservation_test/water_conservation_exp.jl
new file mode 100644
index 0000000..fe3b618
--- /dev/null
+++ b/water_conservation_test/water_conservation_exp.jl
@@ -0,0 +1,178 @@
+using ClimaCore
+import ClimaComms
+ClimaComms.@import_required_backends
+import SciMLBase
+import ClimaTimeSteppers as CTS
+using Plots
+using Statistics
+using DelimitedFiles
+
+using ClimaLand
+using ClimaLand.Soil
+using ClimaLand.Domains: Column
+import ClimaUtilities.OutputPathGenerator: generate_output_path
+
+rmse(v1, v2) = sqrt(mean((v1 .- v2) .^ 2))
+
+# Define simulation times
+t_start = Float64(0)
+dt = Float64(0.001)
+t_end = Float64(1e6)
+
+FT = Float64
+
+stepper = ClimaLSM.RK4()
+ode_algo = CTS.ExplicitAlgorithm(stepper)
+
+# van Genuchten parameters for clay (from Bonan 2019 supplemental program 8.2)
+ν = FT(0.495)
+K_sat = FT(0.0443 / 3600 / 100) # m/s
+vg_n = FT(1.43)
+vg_α = FT(0.026 * 100) # inverse meters
+θ_r = FT(0.124)
+S_s = FT(1e-3) #inverse meters
+hcm = vanGenuchten{FT}(; α = vg_α, n = vg_n)
+
+zmax = FT(0)
+zmin = FT(-0.5)
+nelems = 50
+
+params = Soil.RichardsParameters(;
+    ν = ν,
+    hydrology_cm = hcm,
+    K_sat = K_sat,
+    S_s = S_s,
+    θ_r = θ_r,
+)
+soil_domain = Column(; zlim = (zmin, zmax), nelements = nelems)
+sources = ()
+
+# Set flux boundary conditions (used for calculating mass balance)
+flux_in = FT(-1e-7)
+top_bc = Soil.WaterFluxBC((p, t) -> flux_in)
+flux_out = FT(0)
+bot_bc = Soil.WaterFluxBC((p, t) -> flux_out)
+
+boundary_fluxes = (; top = top_bc, bottom = bot_bc)
+
+soil = Soil.RichardsModel{FT}(;
+    parameters = params,
+    domain = soil_domain,
+    boundary_conditions = boundary_fluxes,
+    sources = sources,
+)
+
+exp_tendency! = make_exp_tendency(soil)
+set_initial_cache! = make_set_initial_cache(soil)
+imp_tendency! = make_imp_tendency(soil)
+jacobian! = make_jacobian(soil)
+
+Y, p, coords = initialize(soil)
+@. Y.soil.ϑ_l = FT(0.24)
+set_initial_cache!(p, Y, FT(0.0))
+
+jac_kwargs =
+    (; jac_prototype = ImplicitEquationJacobian(Y), Wfact = jacobian!)
+
+prob = SciMLBase.ODEProblem(
+    CTS.ClimaODEFunction(T_exp! = exp_tendency!),
+    Y,
+    (t_start, t_end),
+    p,
+)
+sol = SciMLBase.solve(prob, ode_algo; dt = dt, saveat = dt)
+
+# Save flux BC reference solution artifact
+savedir = ARGS[1]
+soln_file = joinpath(savedir, "ref_soln_flux.csv")
+open((soln_file), "w") do io
+    writedlm(io, parent(sol.u[end].soil.ϑ_l), ',')
+end
+
+# TODO remove this after running once
+    # Calculate water mass balance over entire simulation
+    mass_end = sum(sol.u[end].soil.ϑ_l)
+    mass_start = sum(sol.u[1].soil.ϑ_l)
+    t_sim = sol.t[end] - sol.t[1]
+    # Flux changes water content every timestep (assumes constant flux_in, flux_out)
+    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
+
+
+# Perform simulation with Dirichlet boundary conditions
+top_state_bc = Soil.MoistureStateBC((p, t) -> ν - 1e-3)
+flux_out = FT(0)
+bot_flux_bc = Soil.WaterFluxBC((p, t) -> flux_out)
+boundary_conds = (; top = top_state_bc, bottom = bot_flux_bc)
+
+soil_dirichlet = Soil.RichardsModel{FT}(;
+    parameters = params,
+    domain = soil_domain,
+    boundary_conditions = boundary_conds,
+    sources = sources,
+)
+
+exp_tendency! = make_exp_tendency(soil_dirichlet)
+set_initial_cache! = make_set_initial_cache(soil_dirichlet)
+imp_tendency! = make_imp_tendency(soil_dirichlet)
+jacobian! = make_jacobian(soil_dirichlet)
+update_aux! = make_update_aux(soil_dirichlet)
+
+rmses_dirichlet = Array{FT}(undef, length(dts))
+mass_errors_dirichlet = Array{FT}(undef, length(dts))
+
+Y, p, coords = initialize(soil_dirichlet)
+@. Y.soil.ϑ_l = FT(0.24)
+set_initial_cache!(p, Y, FT(0.0))
+
+jac_kwargs =
+    (; jac_prototype = ImplicitEquationJacobian(Y), Wfact = jacobian!)
+
+prob = SciMLBase.ODEProblem(
+    CTS.ClimaODEFunction(T_exp! = exp_tendency!),
+    Y,
+    (t_start, t_end),
+    p,
+)
+# TODO remove saveat, sv after running once
+saveat = Array(t_start:dt:t_end)
+sv = (;
+    t = Array{Int64}(undef, length(saveat)),
+    saveval = Array{NamedTuple}(undef, length(saveat)),
+)
+cb = ClimaLand.NonInterpSavingCallback(sv, saveat)
+sol = SciMLBase.solve(
+    prob,
+    ode_algo;
+    dt = dt,
+    callback = cb,
+    adaptive = false,
+    saveat = saveat,
+)
+
+# TODO remove this after running once
+    # Calculate water mass balance over entire simulation
+    # Because we use Backward Euler, compute fluxes at times[2:end]
+    flux_in_sim =
+        [parent(sv.saveval[k].soil.top_bc)[1] for k in 2:length(sv.saveval)]
+
+
+    mass_end = sum(sol.u[end].soil.ϑ_l)
+    mass_start = sum(sol.u[1].soil.ϑ_l)
+    t_sim = sol.t[end] - sol.t[1]
+    mass_change_exp = -(sum(flux_in_sim) * dt - flux_out * t_sim)
+    mass_change_actual = mass_end - mass_start
+    relerr = abs(mass_change_actual - mass_change_exp) / mass_change_exp
+    @assert relerr < 1e11 * eps(FT)
+    mass_errors_dirichlet[i] = relerr
+
+# Save Dirichlet BC reference solution artifact
+soln_file_dirichlet = joinpath(savedir, "ref_soln_dirichlet.csv")
+open((soln_file_dirichlet), "w") do io
+    writedlm(io, parent(sol.u[end].soil.ϑ_l), ',')
+    writedlm(io, parent(sol.u[end].soil.ϑ_l), ',')
+end