diff --git a/.buildkite/longruns_gpu/pipeline.yml b/.buildkite/longruns_gpu/pipeline.yml
index 568e2130c5..815901f0ae 100644
--- a/.buildkite/longruns_gpu/pipeline.yml
+++ b/.buildkite/longruns_gpu/pipeline.yml
@@ -3,8 +3,6 @@ agents:
   slurm_mem: 8G
   modules: climacommon/2024_10_09
 
-timeout_in_minutes: 1440
-
 steps:
   - label: "init :GPU:"
     key: "init_gpu_env"
@@ -34,10 +32,12 @@ steps:
       - label: ":snow_capped_mountain: Snowy Land"
         command:
           - julia --color=yes --project=.buildkite experiments/long_runs/snowy_land.jl
-        artifact_paths: "snowy_land_longrun_gpu/*png"
+        artifact_paths:
+          - "snowy_land_longrun_gpu/*png"
+          - "snowy_land_longrun_gpu/*pdf"
         agents:
           slurm_gpus: 1
-          slurm_time: 03:30:00
+          slurm_time: 15:00:00
         env:
           CLIMACOMMS_DEVICE: "CUDA"
 
@@ -47,7 +47,7 @@ steps:
         artifact_paths: "land_longrun_gpu/*pdf"
         agents:
           slurm_gpus: 1
-          slurm_time: 03:00:00
+          slurm_time: 15:00:00
         env:
           CLIMACOMMS_DEVICE: "CUDA"
 
@@ -67,7 +67,7 @@ steps:
         artifact_paths: "soil_longrun_gpu/*pdf"
         agents:
           slurm_gpus: 1
-          slurm_time: 01:30:00
+          slurm_time: 15:00:00
         env:
           CLIMACOMMS_DEVICE: "CUDA"
 
diff --git a/docs/Manifest-v1.11.toml b/docs/Manifest-v1.11.toml
index b0540ba25b..6c38dbb708 100644
--- a/docs/Manifest-v1.11.toml
+++ b/docs/Manifest-v1.11.toml
@@ -2,7 +2,7 @@
 
 julia_version = "1.11.1"
 manifest_format = "2.0"
-project_hash = "a5dc443666e7b6f001f45671e08cb0a420c01a66"
+project_hash = "15fe7e962a3bea15852e19f92138be4db325cca0"
 
 [[deps.ADTypes]]
 git-tree-sha1 = "30bb95a372787af850addf28ac937f1be7b79173"
@@ -26,16 +26,6 @@ git-tree-sha1 = "574baf8110975760d391c710b6341da1afa48d8c"
 uuid = "a4c015fc-c6ff-483c-b24f-f7ea428134e9"
 version = "0.0.1"
 
-[[deps.About]]
-deps = ["InteractiveUtils", "JuliaSyntaxHighlighting", "PrecompileTools", "StyledStrings"]
-git-tree-sha1 = "efbf5b623b7ee2a41ce5aed6299aa62ab7f2d5b9"
-uuid = "69d22d85-9f48-4c46-bbbe-7ad8341ff72a"
-version = "1.0.1"
-weakdeps = ["Pkg"]
-
-    [deps.About.extensions]
-    PkgExt = "Pkg"
-
 [[deps.AbstractFFTs]]
 deps = ["LinearAlgebra"]
 git-tree-sha1 = "d92ad398961a3ed262d8bf04a1a2b8340f915fef"
@@ -1469,17 +1459,6 @@ git-tree-sha1 = "af433a10f3942e882d3c671aacb203e006a5808f"
 uuid = "9c1d0b0a-7046-5b2e-a33f-ea22f176ac7e"
 version = "0.2.1+0"
 
-[[deps.JuliaSyntax]]
-git-tree-sha1 = "937da4713526b96ac9a178e2035019d3b78ead4a"
-uuid = "70703baa-626e-46a2-a12c-08ffd08c73b4"
-version = "0.4.10"
-
-[[deps.JuliaSyntaxHighlighting]]
-deps = ["JuliaSyntax", "StyledStrings"]
-git-tree-sha1 = "19ecee1ea81c60156486a92b062e443b6bba60b7"
-uuid = "ac6e5ff7-fb65-4e79-a425-ec3bc9c03011"
-version = "0.1.0"
-
 [[deps.JuliaVariables]]
 deps = ["MLStyle", "NameResolution"]
 git-tree-sha1 = "49fb3cb53362ddadb4415e9b73926d6b40709e70"
diff --git a/docs/Manifest.toml b/docs/Manifest.toml
index 837977b6e3..41bfaf5000 100644
--- a/docs/Manifest.toml
+++ b/docs/Manifest.toml
@@ -2,7 +2,7 @@
 
 julia_version = "1.10.5"
 manifest_format = "2.0"
-project_hash = "17901eb92a5662b053574516009767eae929e671"
+project_hash = "0271f9782f48e923fcefcee2c42babb40262588c"
 
 [[deps.ADTypes]]
 git-tree-sha1 = "30bb95a372787af850addf28ac937f1be7b79173"
@@ -26,16 +26,6 @@ git-tree-sha1 = "574baf8110975760d391c710b6341da1afa48d8c"
 uuid = "a4c015fc-c6ff-483c-b24f-f7ea428134e9"
 version = "0.0.1"
 
-[[deps.About]]
-deps = ["InteractiveUtils", "JuliaSyntaxHighlighting", "PrecompileTools", "StyledStrings"]
-git-tree-sha1 = "efbf5b623b7ee2a41ce5aed6299aa62ab7f2d5b9"
-uuid = "69d22d85-9f48-4c46-bbbe-7ad8341ff72a"
-version = "1.0.1"
-weakdeps = ["Pkg"]
-
-    [deps.About.extensions]
-    PkgExt = "Pkg"
-
 [[deps.AbstractFFTs]]
 deps = ["LinearAlgebra"]
 git-tree-sha1 = "d92ad398961a3ed262d8bf04a1a2b8340f915fef"
@@ -1461,17 +1451,6 @@ git-tree-sha1 = "af433a10f3942e882d3c671aacb203e006a5808f"
 uuid = "9c1d0b0a-7046-5b2e-a33f-ea22f176ac7e"
 version = "0.2.1+0"
 
-[[deps.JuliaSyntax]]
-git-tree-sha1 = "937da4713526b96ac9a178e2035019d3b78ead4a"
-uuid = "70703baa-626e-46a2-a12c-08ffd08c73b4"
-version = "0.4.10"
-
-[[deps.JuliaSyntaxHighlighting]]
-deps = ["JuliaSyntax", "StyledStrings"]
-git-tree-sha1 = "19ecee1ea81c60156486a92b062e443b6bba60b7"
-uuid = "ac6e5ff7-fb65-4e79-a425-ec3bc9c03011"
-version = "0.1.0"
-
 [[deps.JuliaVariables]]
 deps = ["MLStyle", "NameResolution"]
 git-tree-sha1 = "49fb3cb53362ddadb4415e9b73926d6b40709e70"
@@ -2742,12 +2721,6 @@ weakdeps = ["Adapt", "GPUArraysCore", "SparseArrays", "StaticArrays"]
     StructArraysSparseArraysExt = "SparseArrays"
     StructArraysStaticArraysExt = "StaticArrays"
 
-[[deps.StyledStrings]]
-deps = ["PrecompileTools", "TOML"]
-git-tree-sha1 = "711c9650010c95814911c2005ea04e70e70e65ed"
-uuid = "f489334b-da3d-4c2e-b8f0-e476e12c162b"
-version = "1.0.3"
-
 [[deps.SuiteSparse]]
 deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
 uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"
diff --git a/docs/Project.toml b/docs/Project.toml
index 15d5463fb5..4e397fa817 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,5 +1,4 @@
 [deps]
-About = "69d22d85-9f48-4c46-bbbe-7ad8341ff72a"
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
diff --git a/experiments/long_runs/land_region.jl b/experiments/long_runs/land_region.jl
index 9656b58264..e38b5e9f05 100644
--- a/experiments/long_runs/land_region.jl
+++ b/experiments/long_runs/land_region.jl
@@ -1,6 +1,7 @@
-# # Global run of land model
+# # Regional run of full land model
 
-# The code sets up and runs the soil/canopy model for 6 hours on a small region of the globe in Southern California,
+# The code sets up and runs the soil/canopy/snow model for 6 hours on a small
+# region of the globe in Southern California,
 # using ERA5 data. In this simulation, we have
 # turned lateral flow off because horizontal boundary conditions and the
 # land/sea mask are not yet supported by ClimaCore.
@@ -8,10 +9,10 @@
 # Simulation Setup
 # Number of spatial elements: 10x10 in horizontal, 15 in vertical
 # Soil depth: 50 m
-# Simulation duration: 365 d
-# Timestep: 900 s
+# Simulation duration: 4 years
+# Timestep: 450 s
 # Timestepper: ARS111
-# Fixed number of iterations: 1
+# Fixed number of iterations: 3
 # Jacobian update: every new timestep
 # Atmos forcing update: every 3 hours
 import SciMLBase
@@ -32,6 +33,7 @@ import ClimaUtilities.ClimaArtifacts: @clima_artifact
 import ClimaParams as CP
 
 using ClimaLand
+using ClimaLand.Snow
 using ClimaLand.Soil
 using ClimaLand.Canopy
 import ClimaLand
@@ -284,15 +286,21 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (10, 10, 15))
         parameters = shared_params,
         domain = ClimaLand.obtain_surface_domain(domain),
     )
+    # Snow model
+    snow_parameters = SnowParameters{FT}(Δt; earth_param_set = earth_param_set)
+    snow_args = (;
+        parameters = snow_parameters,
+        domain = ClimaLand.obtain_surface_domain(domain),
+    )
+    snow_model_type = Snow.SnowModel
 
-    # Integrated plant hydraulics and soil model
     land_input = (
         atmos = atmos,
         radiation = radiation,
         runoff = runoff_model,
         soil_organic_carbon = Csom,
     )
-    land = SoilCanopyModel{FT}(;
+    land = LandModel{FT}(;
         soilco2_type = soilco2_type,
         soilco2_args = soilco2_args,
         land_args = land_input,
@@ -301,11 +309,17 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (10, 10, 15))
         canopy_component_types = canopy_component_types,
         canopy_component_args = canopy_component_args,
         canopy_model_args = canopy_model_args,
+        snow_args = snow_args,
+        snow_model_type = snow_model_type,
     )
 
     Y, p, cds = initialize(land)
 
     init_soil(ν, θ_r) = θ_r + (ν - θ_r) / 2
+
+    Y.snow.S .= 0.0
+    Y.snow.U .= 0.0
+
     Y.soil.ϑ_l .= init_soil.(ν, θ_r)
     Y.soil.θ_i .= FT(0.0)
     T = FT(276.85)
@@ -381,10 +395,10 @@ function setup_and_solve_problem(; greet = false)
 
     t0 = 0.0
     tf = 60 * 60.0 * 24 * 365
-    Δt = 900.0
+    Δt = 450.0
     nelements = (10, 10, 15)
     if greet
-        @info "Run: Regional Soil-Canopy Model"
+        @info "Run: Regional Soil-Canopy-Snow Model"
         @info "Resolution: $nelements"
         @info "Timestep: $Δt s"
         @info "Duration: $(tf - t0) s"
diff --git a/experiments/long_runs/snowy_land.jl b/experiments/long_runs/snowy_land.jl
index 61bbd67e42..21ebcd8790 100644
--- a/experiments/long_runs/snowy_land.jl
+++ b/experiments/long_runs/snowy_land.jl
@@ -13,7 +13,7 @@
 # Timestep: 450 s
 # Timestepper: ARS111
 # Fixed number of iterations: 3
-# Jacobian update: every Newton iteration
+# Jacobian update: every new Newton iteration
 # Atmos forcing update: every 3 hours
 import SciMLBase
 import ClimaComms
@@ -29,6 +29,7 @@ import ClimaUtilities
 
 import ClimaUtilities.TimeVaryingInputs:
     TimeVaryingInput, LinearInterpolation, PeriodicCalendar
+import ClimaUtilities.SpaceVaryingInputs: SpaceVaryingInput
 import ClimaUtilities.ClimaArtifacts: @clima_artifact
 import ClimaParams as CP
 
@@ -45,6 +46,8 @@ import GeoMakie
 using Dates
 import NCDatasets
 
+using Poppler_jll: pdfunite
+
 const FT = Float64;
 time_interpolation_method = LinearInterpolation(PeriodicCalendar())
 context = ClimaComms.context()
@@ -386,7 +389,7 @@ end
 function setup_and_solve_problem(; greet = false)
 
     t0 = 0.0
-    tf = 60 * 60.0 * 24 * 365
+    tf = 60 * 60.0 * 24 * 365 * 1
     Δt = 450.0
     nelements = (101, 15)
     if greet
@@ -460,3 +463,36 @@ for (group_id, group) in
 
     CairoMakie.save(joinpath(root_path, "$(group_name).png"), fig)
 end
+
+short_names = ["gpp", "swc", "et", "ct", "swe", "si"]
+
+mktempdir(root_path) do tmpdir
+    for short_name in short_names
+        var = get(simdir; short_name)
+        times = [ClimaAnalysis.times(var)[end]]
+        for t in times
+            fig = CairoMakie.Figure(size = (600, 400))
+            kwargs = ClimaAnalysis.has_altitude(var) ? Dict(:z => 1) : Dict()
+            viz.heatmap2D_on_globe!(
+                fig,
+                ClimaAnalysis.slice(var, time = t; kwargs...),
+                mask = viz.oceanmask(),
+                more_kwargs = Dict(
+                    :mask => ClimaAnalysis.Utils.kwargs(color = :white),
+                    :plot => ClimaAnalysis.Utils.kwargs(rasterize = true),
+                ),
+            )
+            CairoMakie.save(joinpath(tmpdir, "$(short_name)_$t.pdf"), fig)
+        end
+    end
+    figures = readdir(tmpdir, join = true)
+    pdfunite() do unite
+        run(
+            Cmd([
+                unite,
+                figures...,
+                joinpath(root_path, "last_year_figures.pdf"),
+            ]),
+        )
+    end
+end
diff --git a/experiments/long_runs/soil.jl b/experiments/long_runs/soil.jl
index 532e7ddf8b..e778799508 100644
--- a/experiments/long_runs/soil.jl
+++ b/experiments/long_runs/soil.jl
@@ -9,7 +9,7 @@
 # Number of spatial elements: 101 1in horizontal, 15 in vertical
 # Soil depth: 50 m
 # Simulation duration: 365 d
-# Timestep: 900 s
+# Timestep: 450 s
 # Timestepper: ARS111
 # Fixed number of iterations: 1
 # Jacobian update: every new timestep
@@ -148,44 +148,9 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         soil_args...,
     )
 
-
     Y, p, cds = initialize(soil)
-    z = ClimaCore.Fields.coordinate_field(cds.subsurface).z
-    lat = ClimaCore.Fields.coordinate_field(cds.subsurface).lat
-    # This function approximates the hydrostatic equilibrium solution in
-    # the vadose and unsaturated regimes by solving for ∂(ψ+z)/∂z = 0,
-    # assuming the transition between the two is at a coordinate of z_∇.
-
-    # The approximation arises because the porosity, residual water fraction,
-    # and van Genuchtem parameters are spatially varying but treated constant
-    # in solving for equilibrium. Finally, we make a plausible but total guess
-    # for the water table depth using the TOPMODEL parameters.
-    function hydrostatic_profile(
-        lat::FT,
-        z::FT,
-        ν::FT,
-        θ_r::FT,
-        α::FT,
-        n::FT,
-        S_s::FT,
-        fmax,
-    ) where {FT}
-        m = 1 - 1 / n
-        zmin = FT(-50.0)
-        zmax = FT(0.0)
-
-        z_∇ = FT(zmin / 5.0 + (zmax - zmin) / 2.5 * (fmax - 0.35) / 0.7)
-        if z > z_∇
-            S = FT((FT(1) + (α * (z - z_∇))^n)^(-m))
-            ϑ_l = S * (ν - θ_r) + θ_r
-        else
-            ϑ_l = -S_s * (z - z_∇) + ν
-        end
-        return FT(ϑ_l)
-    end
-    vg_α = hydrology_cm.α
-    vg_n = hydrology_cm.n
-    Y.soil.ϑ_l .= hydrostatic_profile.(lat, z, ν, θ_r, vg_α, vg_n, S_s, f_max)
+    init_soil(ν, θ_r) = θ_r + (ν - θ_r) / 2
+    Y.soil.ϑ_l .= init_soil.(ν, θ_r)
     Y.soil.θ_i .= FT(0.0)
     T = FT(276.85)
     ρc_s =
@@ -256,7 +221,7 @@ function setup_and_solve_problem(; greet = false)
 
     t0 = 0.0
     tf = 60 * 60.0 * 24 * 365
-    Δt = 900.0
+    Δt = 450.0
     nelements = (101, 15)
     if greet
         @info "Run: Global Soil Model"
diff --git a/src/standalone/Snow/snow_parameterizations.jl b/src/standalone/Snow/snow_parameterizations.jl
index 0217b53d5e..7bc9d56199 100644
--- a/src/standalone/Snow/snow_parameterizations.jl
+++ b/src/standalone/Snow/snow_parameterizations.jl
@@ -30,19 +30,20 @@ end
         p,
     ) where {FT}
 
-Returns the surface height of the `Snow` model; this returns the depth
-of the snow and hence implicitly treats the surface (ground) elevation as
-at zero.
+Returns the surface height of the `Snow` model; surface (ground) elevation
+and ignores snow depth (CLM).
 
 Once topography or land ice is incorporated, this will need to change to
-z_sfc + land_ice_depth + snow_depth. Note that land ice can 
+z_sfc + land_ice_depth. Note that land ice can
 be ~1-3 km thick on Greenland/
 
-In order to compute surface fluxes, this cannot be large than the 
-height of the atmosphere measurement location (z_atmos > z_land always). 
+In order to compute surface fluxes, this cannot be large than the
+height of the atmosphere measurement location (z_atmos > z_land always).
+
+This assumes that the surface elevation is zero.
 """
 function ClimaLand.surface_height(model::SnowModel{FT}, Y, p) where {FT}
-    return p.snow.z
+    return FT(0)
 end
 
 
@@ -67,7 +68,7 @@ end
 """
     ClimaLand.surface_temperature(model::SnowModel, Y, p)
 
-a helper function which returns the surface temperature for the snow 
+a helper function which returns the surface temperature for the snow
 model, which is stored in the aux state.
 """
 function ClimaLand.surface_temperature(model::SnowModel, Y, p, t)
@@ -318,7 +319,7 @@ end
 """
     energy_from_q_l_and_swe(S::FT, q_l::FT, parameters) where {FT}
 
-A helper function for compute the snow energy per unit area, given snow 
+A helper function for compute the snow energy per unit area, given snow
 water equivalent S, liquid fraction q_l,  and snow model parameters.
 
 Note that liquid water can only exist at the freezing point in this model,
@@ -340,7 +341,7 @@ end
 """
     energy_from_T_and_swe(S::FT, T::FT, parameters) where {FT}
 
-A helper function for compute the snow energy per unit area, given snow 
+A helper function for compute the snow energy per unit area, given snow
 water equivalent S, bulk temperature T, and snow model parameters.
 
 If T = T_freeze, we return the energy as if q_l = 0.
diff --git a/src/standalone/Vegetation/canopy_boundary_fluxes.jl b/src/standalone/Vegetation/canopy_boundary_fluxes.jl
index 781bf3f1cc..309dfbdb0d 100644
--- a/src/standalone/Vegetation/canopy_boundary_fluxes.jl
+++ b/src/standalone/Vegetation/canopy_boundary_fluxes.jl
@@ -31,7 +31,7 @@ to as ``boundary conditions", at the surface and
 bottom of the canopy, for water and energy.
 
 These fluxes include turbulent surface fluxes
-computed with Monin-Obukhov theory, radiative fluxes, 
+computed with Monin-Obukhov theory, radiative fluxes,
 and root extraction.
 $(DocStringExtensions.FIELDS)
 """
@@ -61,7 +61,7 @@ end
 
 An outer constructor for `AtmosDrivenCanopyBC` which is
 intended for use as a default when running standlone canopy
-models. 
+models.
 
 This is also checks the logic that:
 - If the `ground` field is Prescribed, :soil should not be a prognostic_land_component
@@ -129,9 +129,11 @@ end
 
 A helper function which returns the surface height for the canopy
 model, which is stored in the parameter struct.
+
+This assumes that the surface elevation is zero.
 """
-function ClimaLand.surface_height(model::CanopyModel, _...)
-    return model.hydraulics.compartment_surfaces[1]
+function ClimaLand.surface_height(model::CanopyModel{FT}, _...) where {FT}
+    return FT(0)
 end
 
 function make_update_boundary_fluxes(canopy::CanopyModel)
@@ -382,7 +384,7 @@ function canopy_turbulent_fluxes_at_a_point(
             states,
             scheme,
         ) * r_ae / (canopy_r_b + r_ae)
-    # The above follows from CLM 5, tech note Equation 5.88, setting H_v = SH and solving to remove T_s, ignoring difference between cp in atmos and above canopy. 
+    # The above follows from CLM 5, tech note Equation 5.88, setting H_v = SH and solving to remove T_s, ignoring difference between cp in atmos and above canopy.
     LH = _LH_v0 * E
 
     # Derivatives