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/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..78392d98a0 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
@@ -150,8 +150,11 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
 
 
     Y, p, cds = initialize(soil)
-    z = ClimaCore.Fields.coordinate_field(cds.subsurface).z
-    lat = ClimaCore.Fields.coordinate_field(cds.subsurface).lat
+    init_soil(ν, θ_r) = θ_r + (ν - θ_r) / 2
+    Y.soil.ϑ_l .= init_soil.(ν, θ_r)
+    Y.soil.θ_i .= FT(0.0)
+    #    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_∇.
@@ -160,33 +163,35 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
     # 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)
+    #=
+        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_∇) + ν
+            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
-        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)
-    Y.soil.θ_i .= FT(0.0)
+        vg_α = hydrology_cm.α
+        vg_n = hydrology_cm.n
+        Y.soil.ϑ_l .= hydrostatic_profile.(lat, z, ν, θ_r, vg_α, vg_n, S_s, f_max)
+        Y.soil.θ_i .= FT(0.0)
+    =#
     T = FT(276.85)
     ρc_s =
         Soil.volumetric_heat_capacity.(
@@ -256,7 +261,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"