diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml
index b1ae24d70e..8185ffecca 100644
--- a/.buildkite/pipeline.yml
+++ b/.buildkite/pipeline.yml
@@ -127,6 +127,14 @@ steps:
         agents:
           slurm_mem: 16G
 
+      - label: "Canopy Implicit Stepping CPU"
+        command: "julia --color=yes --project=.buildkite experiments/standalone/Vegetation/timestep_test.jl"
+        artifact_paths: "experiments/standalone/Vegetation/timestep_test/*"
+
+      - label: "SoilCanopy Implicit Stepping CPU"
+        command: "julia --color=yes --project=.buildkite experiments/integrated/performance/integrated_timestep_test.jl"
+        artifact_paths: "experiments/integrated/performance/integrated_timestep_test/*"
+
   - group: "Experiments on GPU"
     steps:
       - label: "Richards Runoff GPU"
diff --git a/docs/src/APIs/canopy/Photosynthesis.md b/docs/src/APIs/canopy/Photosynthesis.md
index 6361344835..5220e3f07f 100644
--- a/docs/src/APIs/canopy/Photosynthesis.md
+++ b/docs/src/APIs/canopy/Photosynthesis.md
@@ -20,8 +20,6 @@ ClimaLand.Canopy.intercellular_co2
 ClimaLand.Canopy.co2_compensation
 ClimaLand.Canopy.rubisco_assimilation
 ClimaLand.Canopy.light_assimilation
-ClimaLand.Canopy.C3
-ClimaLand.Canopy.C4
 ClimaLand.Canopy.max_electron_transport
 ClimaLand.Canopy.electron_transport
 ClimaLand.Canopy.net_photosynthesis
diff --git a/experiments/benchmarks/land.jl b/experiments/benchmarks/land.jl
index 07682b3f23..a615fa49a3 100644
--- a/experiments/benchmarks/land.jl
+++ b/experiments/benchmarks/land.jl
@@ -606,7 +606,7 @@ function setup_simulation(; greet = false)
     Δt = 900.0
     nelements = (101, 15)
     if greet
-        @info "Run: Global RichardsModel"
+        @info "Run: Global Land Model"
         @info "Resolution: $nelements"
         @info "Timestep: $Δt s"
         @info "Duration: $(tf - t0) s"
@@ -615,12 +615,12 @@ function setup_simulation(; greet = false)
     prob, cb = setup_prob(t0, tf, Δt; nelements)
 
     # Define timestepper and ODE algorithm
-    stepper = CTS.ARS343()
+    stepper = CTS.ARS111()
     ode_algo = CTS.IMEXAlgorithm(
         stepper,
         CTS.NewtonsMethod(
-            max_iters = 1,
-            update_j = CTS.UpdateEvery(CTS.NewTimeStep),
+            max_iters = 6,
+            update_j = CTS.UpdateEvery(CTS.NewNewtonIteration),
         ),
     )
     return prob, ode_algo, Δt, cb
@@ -632,7 +632,7 @@ SciMLBase.solve(prob, ode_algo; dt = Δt, callback = cb)
 
 @info "Starting profiling"
 # Stop when we profile for MAX_PROFILING_TIME_SECONDS or MAX_PROFILING_SAMPLES
-MAX_PROFILING_TIME_SECONDS = 1000
+MAX_PROFILING_TIME_SECONDS = 500
 MAX_PROFILING_SAMPLES = 100
 time_now = time()
 timings_s = Float64[]
diff --git a/experiments/integrated/fluxnet/US-Ha1/US-Ha1_simulation.jl b/experiments/integrated/fluxnet/US-Ha1/US-Ha1_simulation.jl
index 6afcd4670d..0241f7f925 100644
--- a/experiments/integrated/fluxnet/US-Ha1/US-Ha1_simulation.jl
+++ b/experiments/integrated/fluxnet/US-Ha1/US-Ha1_simulation.jl
@@ -20,7 +20,4 @@ h_stem = FT(14) # m
 t0 = Float64(120 * 3600 * 24)# start mid year to avoid snow
 
 # Time step size:
-dt = Float64(40)
-
-# Number of timesteps between saving output:
-n = 45
+dt = Float64(450)
diff --git a/experiments/integrated/fluxnet/US-MOz/US-MOz_simulation.jl b/experiments/integrated/fluxnet/US-MOz/US-MOz_simulation.jl
index b229ab5024..6299ab1d63 100644
--- a/experiments/integrated/fluxnet/US-MOz/US-MOz_simulation.jl
+++ b/experiments/integrated/fluxnet/US-MOz/US-MOz_simulation.jl
@@ -20,7 +20,4 @@ h_leaf = FT(9.5) # m
 t0 = Float64(120 * 3600 * 24)# start mid year to avoid snow
 
 # Time step size:
-dt = Float64(450)
-
-# Number of timesteps between saving output:
-n = 4
+dt = Float64(900)
diff --git a/experiments/integrated/fluxnet/US-NR1/US-NR1_simulation.jl b/experiments/integrated/fluxnet/US-NR1/US-NR1_simulation.jl
index e27c0b4827..644b836bfa 100644
--- a/experiments/integrated/fluxnet/US-NR1/US-NR1_simulation.jl
+++ b/experiments/integrated/fluxnet/US-NR1/US-NR1_simulation.jl
@@ -20,7 +20,4 @@ h_stem = FT(7.5) # m
 t0 = Float64(120 * 3600 * 24)# start mid year to avoid snow
 
 # Time step size:
-dt = Float64(40)
-
-# Number of timesteps between saving output:
-n = 45
+dt = Float64(450)
diff --git a/experiments/integrated/fluxnet/US-Var/US-Var_simulation.jl b/experiments/integrated/fluxnet/US-Var/US-Var_simulation.jl
index 1ba85752a1..563cf22bc6 100644
--- a/experiments/integrated/fluxnet/US-Var/US-Var_simulation.jl
+++ b/experiments/integrated/fluxnet/US-Var/US-Var_simulation.jl
@@ -20,7 +20,4 @@ h_stem = FT(0) # m
 t0 = Float64(21 * 3600 * 24)# start day 21 of the year
 
 # Time step size:
-dt = Float64(20)
-
-# Number of timesteps between saving output:
-n = 45
+dt = Float64(900)
diff --git a/experiments/integrated/fluxnet/fluxnet_simulation.jl b/experiments/integrated/fluxnet/fluxnet_simulation.jl
index 818e9af174..cdf16b8c90 100644
--- a/experiments/integrated/fluxnet/fluxnet_simulation.jl
+++ b/experiments/integrated/fluxnet/fluxnet_simulation.jl
@@ -6,14 +6,15 @@ N_spinup_days = 30
 N_days = N_spinup_days + 30
 tf = Float64(t0 + 3600 * 24 * N_days)
 t_spinup = Float64(t0 + N_spinup_days * 3600 * 24)
-saveat = Array(t_spinup:(n * dt):tf)
+dt_save = 3600.0
+saveat = Array(t_spinup:dt_save:tf)
 
 # Set up timestepper
-timestepper = CTS.ARS343();
+timestepper = CTS.ARS111();
 ode_algo = CTS.IMEXAlgorithm(
     timestepper,
     CTS.NewtonsMethod(
-        max_iters = 1,
+        max_iters = 6,
         update_j = CTS.UpdateEvery(CTS.NewNewtonIteration),
     ),
 );
diff --git a/experiments/integrated/fluxnet/ozark_pft.jl b/experiments/integrated/fluxnet/ozark_pft.jl
index 4d7be37565..a9618244e1 100644
--- a/experiments/integrated/fluxnet/ozark_pft.jl
+++ b/experiments/integrated/fluxnet/ozark_pft.jl
@@ -348,7 +348,7 @@ num_days = N_days - N_spinup_days
 
 # Time series of model and data outputs
 data_times = [0:DATA_DT:(num_days * S_PER_DAY);]
-model_times = [0:(n * dt):(num_days * S_PER_DAY);]
+model_times = [0:dt_save:(num_days * S_PER_DAY);]
 
 # Plot model diurnal cycles without data comparisons
 # Autotrophic Respiration
@@ -357,7 +357,15 @@ AR =
         parent(sv.saveval[k].canopy.autotrophic_respiration.Ra)[1] for
         k in 1:length(sv.saveval)
     ] .* 1e6
-plot_daily_avg("AutoResp", AR, dt * n, num_days, "μmol/m^2/s", savedir, "Model")
+plot_daily_avg(
+    "AutoResp",
+    AR,
+    dt_save,
+    num_days,
+    "μmol/m^2/s",
+    savedir,
+    "Model",
+)
 
 # Plot all comparisons of model diurnal cycles to data diurnal cycles
 # GPP
@@ -371,7 +379,7 @@ if drivers.GPP.status == absent
     plot_daily_avg(
         "GPP",
         model_GPP,
-        dt * n,
+        dt_save,
         num_days,
         "μmol/m^2/s",
         savedir,
@@ -383,7 +391,7 @@ else
     plot_avg_comp(
         "GPP",
         model_GPP,
-        dt * n,
+        dt_save,
         GPP_data,
         FT(DATA_DT),
         num_days,
@@ -398,7 +406,7 @@ if drivers.SW_OUT.status == absent
     plot_daily_avg(
         "SW OUT",
         SW_out_model,
-        dt * n,
+        dt_save,
         num_days,
         "w/m^2",
         savedir,
@@ -411,7 +419,7 @@ else
     plot_avg_comp(
         "SW OUT",
         SW_out_model,
-        dt * n,
+        dt_save,
         SW_out_data,
         FT(DATA_DT),
         num_days,
@@ -426,7 +434,7 @@ if drivers.LW_OUT.status == absent
     plot_daily_avg(
         "LW OUT",
         LW_out_model,
-        dt * n,
+        dt_save,
         num_days,
         "w/m^2",
         savedir,
@@ -439,7 +447,7 @@ else
     plot_avg_comp(
         "LW OUT",
         LW_out_model,
-        dt * n,
+        dt_save,
         LW_out_data,
         FT(DATA_DT),
         num_days,
@@ -463,7 +471,15 @@ E =
     ] .* (1e3 * 24 * 3600)
 ET_model = T .+ E
 if drivers.LE.status == absent
-    plot_daily_avg("ET", ET_model, dt * n, num_days, "mm/day", savedir, "Model")
+    plot_daily_avg(
+        "ET",
+        ET_model,
+        dt_save,
+        num_days,
+        "mm/day",
+        savedir,
+        "Model",
+    )
 else
     measured_T =
         drivers.LE.values ./ (LP.LH_v0(earth_param_set) * 1000) .*
@@ -472,7 +488,7 @@ else
     plot_avg_comp(
         "ET",
         ET_model,
-        dt * n,
+        dt_save,
         ET_data,
         FT(DATA_DT),
         num_days,
@@ -492,7 +508,7 @@ if drivers.H.status == absent
     plot_daily_avg(
         "SHF",
         SHF_model,
-        dt * n,
+        dt_save,
         num_days,
         "w/m^2",
         savedir,
@@ -503,7 +519,7 @@ else
     plot_avg_comp(
         "SHF",
         SHF_model,
-        dt * n,
+        dt_save,
         SHF_data,
         FT(DATA_DT),
         N_days - N_spinup_days,
@@ -520,13 +536,13 @@ LHF_canopy =
     [parent(sv.saveval[k].canopy.energy.lhf)[1] for k in 1:length(sol.t)]
 LHF_model = LHF_soil + LHF_canopy
 if drivers.LE.status == absent
-    plot_daily_avg("LHF", LHF_model, dt * n, num_days, "w/m^2", savedir)
+    plot_daily_avg("LHF", LHF_model, dt_save, num_days, "w/m^2", savedir)
 else
     LHF_data = drivers.LE.values[Int64(t_spinup ÷ DATA_DT):Int64(tf ÷ DATA_DT)]
     plot_avg_comp(
         "LHF",
         LHF_model,
-        dt * n,
+        dt_save,
         LHF_data,
         FT(DATA_DT),
         N_days - N_spinup_days,
diff --git a/experiments/integrated/fluxnet/run_fluxnet.jl b/experiments/integrated/fluxnet/run_fluxnet.jl
index e41095c316..31899854b4 100644
--- a/experiments/integrated/fluxnet/run_fluxnet.jl
+++ b/experiments/integrated/fluxnet/run_fluxnet.jl
@@ -292,7 +292,7 @@ num_days = N_days - N_spinup_days
 
 # Time series of model and data outputs
 data_times = [0:DATA_DT:(num_days * S_PER_DAY);]
-model_times = [0:(n * dt):(num_days * S_PER_DAY);]
+model_times = [0:dt_save:(num_days * S_PER_DAY);]
 
 # Plot model diurnal cycles without data comparisons
 # Autotrophic Respiration
@@ -301,7 +301,15 @@ AR =
         parent(sv.saveval[k].canopy.autotrophic_respiration.Ra)[1] for
         k in 1:length(sv.saveval)
     ] .* 1e6
-plot_daily_avg("AutoResp", AR, dt * n, num_days, "μmol/m^2/s", savedir, "Model")
+plot_daily_avg(
+    "AutoResp",
+    AR,
+    dt_save,
+    num_days,
+    "μmol/m^2/s",
+    savedir,
+    "Model",
+)
 
 # Plot all comparisons of model diurnal cycles to data diurnal cycles
 # GPP
@@ -315,7 +323,7 @@ if drivers.GPP.status == absent
     plot_daily_avg(
         "GPP",
         model_GPP,
-        dt * n,
+        dt_save,
         num_days,
         "μmol/m^2/s",
         savedir,
@@ -327,7 +335,7 @@ else
     plot_avg_comp(
         "GPP",
         model_GPP,
-        dt * n,
+        dt_save,
         GPP_data,
         FT(DATA_DT),
         num_days,
@@ -342,7 +350,7 @@ if drivers.SW_OUT.status == absent
     plot_daily_avg(
         "SW OUT",
         SW_out_model,
-        dt * n,
+        dt_save,
         num_days,
         "w/m^2",
         savedir,
@@ -355,7 +363,7 @@ else
     plot_avg_comp(
         "SW OUT",
         SW_out_model,
-        dt * n,
+        dt_save,
         SW_out_data,
         FT(DATA_DT),
         num_days,
@@ -370,7 +378,7 @@ if drivers.LW_OUT.status == absent
     plot_daily_avg(
         "LW OUT",
         LW_out_model,
-        dt * n,
+        dt_save,
         num_days,
         "w/m^2",
         savedir,
@@ -383,7 +391,7 @@ else
     plot_avg_comp(
         "LW OUT",
         LW_out_model,
-        dt * n,
+        dt_save,
         LW_out_data,
         FT(DATA_DT),
         num_days,
@@ -407,7 +415,15 @@ E =
     ] .* (1e3 * 24 * 3600)
 ET_model = T .+ E
 if drivers.LE.status == absent
-    plot_daily_avg("ET", ET_model, dt * n, num_days, "mm/day", savedir, "Model")
+    plot_daily_avg(
+        "ET",
+        ET_model,
+        dt_save,
+        num_days,
+        "mm/day",
+        savedir,
+        "Model",
+    )
 else
     measured_T =
         drivers.LE.values ./ (LP.LH_v0(earth_param_set) * 1000) .*
@@ -416,7 +432,7 @@ else
     plot_avg_comp(
         "ET",
         ET_model,
-        dt * n,
+        dt_save,
         ET_data,
         FT(DATA_DT),
         num_days,
@@ -436,7 +452,7 @@ if drivers.H.status == absent
     plot_daily_avg(
         "SHF",
         SHF_model,
-        dt * n,
+        dt_save,
         num_days,
         "w/m^2",
         savedir,
@@ -447,7 +463,7 @@ else
     plot_avg_comp(
         "SHF",
         SHF_model,
-        dt * n,
+        dt_save,
         SHF_data,
         FT(DATA_DT),
         N_days - N_spinup_days,
@@ -464,13 +480,13 @@ LHF_canopy =
     [parent(sv.saveval[k].canopy.energy.lhf)[1] for k in 1:length(sol.t)]
 LHF_model = LHF_soil + LHF_canopy
 if drivers.LE.status == absent
-    plot_daily_avg("LHF", LHF_model, dt * n, num_days, "w/m^2", savedir)
+    plot_daily_avg("LHF", LHF_model, dt_save, num_days, "w/m^2", savedir)
 else
     LHF_data = drivers.LE.values[Int64(t_spinup ÷ DATA_DT):Int64(tf ÷ DATA_DT)]
     plot_avg_comp(
         "LHF",
         LHF_model,
-        dt * n,
+        dt_save,
         LHF_data,
         FT(DATA_DT),
         N_days - N_spinup_days,
diff --git a/experiments/integrated/performance/integrated_timestep_test.jl b/experiments/integrated/performance/integrated_timestep_test.jl
new file mode 100644
index 0000000000..f0022ae62a
--- /dev/null
+++ b/experiments/integrated/performance/integrated_timestep_test.jl
@@ -0,0 +1,506 @@
+# # Timestep test for integrated SoilCanopyModel
+
+# This experiment runs the standalone canopy model for 12.5 hours on a spherical domain
+# and outputs useful plots related to the timestepping and stability of
+# the simulation using different timesteps. This information can be used to
+# determine the optimal timestep for a this simulation setup.
+
+# The simulation is run multiple times: first, with a small reference timestep,
+# then with increasing timestep sizes. The goal here is to see which timesteps
+# are small enough to produce a stable simulation, and which timesteps are too
+# large to produce accurate results.
+
+# The runs with larger timesteps are compared with the reference timestep run
+# using multiple metrics, including:
+# - mean, 99%th, and 95th% error over the course of the simulation
+# - T at the end of the simulation (used to test convergence with each dt)
+# - T throughout the entire simulation
+
+# Note that the simulation is run for 12.5 hours to allow us to test the
+# convergence behavior of the solver.
+# Convergence behavior must be assessed when T is actively changing, rather
+# than when it is in equilibrium. T changes in the time period after a driver
+# update, which occurs every 3 hours (so exactly at 12 hours). Running for
+# just longer than 12 hours allows us to observe convergence behavior at a time
+# when T is actively changing. See the discussion in
+# https://github.com/CliMA/ClimaLand.jl/pull/675 for more information.
+
+# Simulation Setup
+# Number of spatial elements: 10 in horizontal, 5 in vertical
+# Soil depth: 5 m
+# Simulation duration: 12.5 hours
+# Timestep: variable, ranging from 0.6s to 3600s
+# Timestepper: ARS111
+# Fixed number of iterations: 6
+# Jacobian update: Every Newton iteration
+# Atmospheric and radiation driver updates: Every 3 hours
+
+import SciMLBase
+import ClimaTimeSteppers as CTS
+using ClimaCore
+using Dates
+using Insolation
+using Statistics
+import StatsBase: percentile
+using CairoMakie
+import ClimaComms
+@static pkgversion(ClimaComms) >= v"0.6" && ClimaComms.@import_required_backends
+import ClimaUtilities.TimeVaryingInputs: TimeVaryingInput
+
+using ClimaLand
+using ClimaLand.Domains
+using ClimaLand.Soil
+using ClimaLand.Soil.Biogeochemistry
+using ClimaLand.Canopy
+using ClimaLand.Canopy.PlantHydraulics
+import ClimaLand
+import ClimaLand.Parameters as LP
+import ClimaParams
+
+
+PROFILING = false
+try
+    import Profile, ProfileCanvas
+    global PROFILING = true
+    @info "ProfileCanvas found, running with profiler"
+catch
+end
+
+function set_initial_conditions(land, t0)
+    Y, p, cds = initialize(land)
+    FT = eltype(Y.soil.ϑ_l)
+    set_initial_cache! = make_set_initial_cache(land)
+
+    Y.soil.ϑ_l = FT(0.3)
+    Y.soil.θ_i = FT(0.0)
+    T_0 = FT(297.5)
+    ρc_s =
+        volumetric_heat_capacity.(
+            Y.soil.ϑ_l,
+            Y.soil.θ_i,
+            land.soil.parameters.ρc_ds,
+            land.soil.parameters.earth_param_set,
+        )
+    Y.soil.ρe_int =
+        volumetric_internal_energy.(
+            Y.soil.θ_i,
+            ρc_s,
+            T_0,
+            land.soil.parameters.earth_param_set,
+        )
+
+    Y.soilco2.C = FT(0.000412) # set to atmospheric co2, mol co2 per mol air
+
+    ψ_stem_0 = FT(-1e5 / 9800)
+    ψ_leaf_0 = FT(-2e5 / 9800)
+    canopy_params = land.canopy.hydraulics.parameters
+    S_l_ini =
+        inverse_water_retention_curve.(
+            canopy_params.retention_model,
+            [ψ_stem_0, ψ_leaf_0],
+            canopy_params.ν,
+            canopy_params.S_s,
+        )
+
+    for i in 1:2
+        Y.canopy.hydraulics.ϑ_l.:($i) .=
+            augmented_liquid_fraction.(canopy_params.ν, S_l_ini[i])
+    end
+
+    Y.canopy.energy.T = FT(297.5)
+    set_initial_cache!(p, Y, t0)
+    return Y, p
+end
+
+context = ClimaComms.context()
+device_suffix =
+    typeof(ClimaComms.context().device) <: ClimaComms.CPUSingleThreaded ?
+    "cpu" : "gpu"
+const FT = Float64
+earth_param_set = LP.LandParameters(FT)
+
+# Set the model domain
+dz_bottom = FT(2.0)
+dz_top = FT(0.2)
+soil_depth = FT(5)
+z_sfc = FT(0)
+n_stem = Int64(1)
+n_leaf = Int64(1)
+h_stem = FT(9) # m
+h_leaf = FT(9.5) # m
+
+h_canopy = h_stem + h_leaf
+compartment_midpoints = [h_stem / 2, h_stem + h_leaf / 2]
+compartment_surfaces = [z_sfc, h_stem, h_canopy]
+
+land_domain = ClimaLand.Domains.SphericalShell(;
+    radius = FT(6.3781e6),
+    depth = soil_depth,
+    nelements = (10, 5),
+    npolynomial = 1,
+    dz_tuple = FT.((dz_bottom, dz_top)),
+);
+canopy_domain = ClimaLand.Domains.obtain_surface_domain(land_domain)
+sfc_cds = ClimaCore.Fields.coordinate_field(land_domain.space.surface)
+
+# First pick the start date and time of the simulation, since time varying input depends on that
+t0 = Float64(0) # start at start date
+start_date = DateTime("202001010000", "yyyymmddHHMM")
+# Time varying input
+LAIfunction = TimeVaryingInput((t) -> 2.0)
+# Atmospheric and radiative forcing
+precip = TimeVaryingInput((t) -> -1.0e-7)
+atmos_q = TimeVaryingInput((t) -> 0.002)
+atmos_T = TimeVaryingInput((t) -> 298.0)
+atmos_p = TimeVaryingInput((t) -> 101320)
+atmos_u = TimeVaryingInput((t) -> 3.0)
+LW_IN = TimeVaryingInput((t) -> 5.67e-8 * 298.0^4)
+SW_IN = TimeVaryingInput((t) -> 500.0)
+snow_precip = TimeVaryingInput((t) -> 0.0)
+atmos_h = FT(32)
+# Construct the drivers. For the start date we will use the UTC time at the
+# start of the simulation
+atmos = ClimaLand.PrescribedAtmosphere(
+    precip,
+    snow_precip,
+    atmos_T,
+    atmos_u,
+    atmos_q,
+    atmos_p,
+    start_date,
+    atmos_h,
+    earth_param_set;
+)
+function zenith_angle(
+    t,
+    start_date;
+    cd_field = sfc_cds,
+    insol_params::Insolation.Parameters.InsolationParameters{FT} = earth_param_set.insol_params,
+) where {FT}
+    # This should be time in UTC
+    current_datetime = start_date + Dates.Second(round(t))
+
+    # Orbital Data uses Float64, so we need to convert to our sim FT
+    d, δ, η_UTC =
+        FT.(
+            Insolation.helper_instantaneous_zenith_angle(
+                current_datetime,
+                start_date,
+                insol_params,
+            )
+        )
+
+    Insolation.instantaneous_zenith_angle.(
+        d,
+        δ,
+        η_UTC,
+        sfc_cds.long,
+        sfc_cds.lat,
+    ).:1
+end
+
+radiation = ClimaLand.PrescribedRadiativeFluxes(
+    FT,
+    SW_IN,
+    LW_IN,
+    start_date,
+    θs = zenith_angle,
+)
+
+# Model parameters
+# Soil parameters
+soil_ν = FT(0.5) # m3/m3
+soil_K_sat = FT(4e-7) # m/
+soil_S_s = FT(1e-3) # 1/m
+soil_vg_n = FT(2.05) # unitless
+soil_vg_α = FT(0.04) # inverse meters
+θ_r = FT(0.067)
+ν_ss_quartz = FT(0.1)
+ν_ss_om = FT(0.1)
+ν_ss_gravel = FT(0.0);
+# Energy Balance model
+ac_canopy = FT(2.5e3)
+# Conductance Model
+g1 = FT(141)
+#Photosynthesis model
+Vcmax25 = FT(9e-5)
+# Plant Hydraulics and general plant parameters
+K_sat_plant = 5e-9
+ψ63 = FT(-4 / 0.0098)
+Weibull_param = FT(4)
+a = FT(0.05 * 0.0098)
+conductivity_model =
+    PlantHydraulics.Weibull{FT}(K_sat_plant, ψ63, Weibull_param)
+retention_model = PlantHydraulics.LinearRetentionCurve{FT}(a)
+plant_ν = FT(1.44e-4)
+SAI = FT(1.0)
+maxLAI = FT(2.0)
+f_root_to_shoot = FT(3.5)
+RAI = maxLAI * f_root_to_shoot
+capacity = plant_ν * maxLAI * h_leaf * FT(1000)
+plant_S_s = FT(1e-2 * 0.0098) # m3/m3/MPa to m3/m3/m
+z0_m = FT(0.13) * h_canopy
+z0_b = FT(0.1) * z0_m
+
+# Set up model
+# Now we set up the model. For the soil model, we pick
+# a model type and model args:
+soil_domain = land_domain
+soil_ps = Soil.EnergyHydrologyParameters(
+    FT;
+    ν = soil_ν,
+    ν_ss_om,
+    ν_ss_quartz,
+    ν_ss_gravel,
+    hydrology_cm = vanGenuchten{FT}(; α = soil_vg_α, n = soil_vg_n),
+    K_sat = soil_K_sat,
+    S_s = soil_S_s,
+    θ_r,
+)
+
+soil_args = (domain = soil_domain, parameters = soil_ps)
+soil_model_type = Soil.EnergyHydrology{FT}
+
+# Soil microbes model
+soilco2_type = Soil.Biogeochemistry.SoilCO2Model{FT}
+
+soilco2_ps = SoilCO2ModelParameters(FT)
+
+# soil microbes args
+Csom = ClimaLand.PrescribedSoilOrganicCarbon{FT}(TimeVaryingInput((t) -> 5))
+
+# Set the soil CO2 BC to being atmospheric CO2
+soilco2_top_bc = Soil.Biogeochemistry.AtmosCO2StateBC()
+soilco2_bot_bc = Soil.Biogeochemistry.SoilCO2StateBC((p, t) -> 0.0)
+soilco2_sources = (MicrobeProduction{FT}(),)
+
+soilco2_boundary_conditions = (; top = soilco2_top_bc, bottom = soilco2_bot_bc)
+
+soilco2_args = (;
+    boundary_conditions = soilco2_boundary_conditions,
+    sources = soilco2_sources,
+    domain = soil_domain,
+    parameters = soilco2_ps,
+)
+
+# Now we set up the canopy model, which we set up by component:
+# Component Types
+canopy_component_types = (;
+    autotrophic_respiration = Canopy.AutotrophicRespirationModel{FT},
+    radiative_transfer = Canopy.TwoStreamModel{FT},
+    photosynthesis = Canopy.FarquharModel{FT},
+    conductance = Canopy.MedlynConductanceModel{FT},
+    hydraulics = Canopy.PlantHydraulicsModel{FT},
+    energy = Canopy.BigLeafEnergyModel{FT},
+)
+# Individual Component arguments
+energy_args = (parameters = Canopy.BigLeafEnergyParameters{FT}(ac_canopy),)
+
+# Set up autotrophic respiration
+autotrophic_respiration_args =
+    (; parameters = AutotrophicRespirationParameters(FT))
+# Set up radiative transfer
+radiative_transfer_args = (; parameters = TwoStreamParameters(FT))
+# Set up conductance
+conductance_args = (; parameters = MedlynConductanceParameters(FT; g1))
+# Set up photosynthesis
+is_c3 = FT(1) # set the photosynthesis mechanism to C3
+photosynthesis_args =
+    (; parameters = FarquharParameters(FT, is_c3; Vcmax25 = Vcmax25))
+# Set up plant hydraulics
+ai_parameterization = PrescribedSiteAreaIndex{FT}(LAIfunction, SAI, RAI)
+function root_distribution(z::T) where {T}
+    return T(1.0 / 0.5 * exp(z / 0.5)) # 1/m
+end
+plant_hydraulics_ps = PlantHydraulics.PlantHydraulicsParameters(;
+    ai_parameterization = ai_parameterization,
+    ν = plant_ν,
+    S_s = plant_S_s,
+    root_distribution = root_distribution,
+    conductivity_model = conductivity_model,
+    retention_model = retention_model,
+)
+plant_hydraulics_args = (
+    parameters = plant_hydraulics_ps,
+    n_stem = n_stem,
+    n_leaf = n_leaf,
+    compartment_midpoints = compartment_midpoints,
+    compartment_surfaces = compartment_surfaces,
+)
+
+# Canopy component args
+canopy_component_args = (;
+    autotrophic_respiration = autotrophic_respiration_args,
+    radiative_transfer = radiative_transfer_args,
+    photosynthesis = photosynthesis_args,
+    conductance = conductance_args,
+    hydraulics = plant_hydraulics_args,
+    energy = energy_args,
+)
+# Other info needed
+shared_params = SharedCanopyParameters{FT, typeof(earth_param_set)}(
+    z0_m,
+    z0_b,
+    earth_param_set,
+)
+
+canopy_model_args = (; parameters = shared_params, domain = canopy_domain)
+
+# Integrated plant hydraulics and soil model
+land_input = (
+    atmos = atmos,
+    radiation = radiation,
+    soil_organic_carbon = Csom,
+    runoff = ClimaLand.Soil.Runoff.SurfaceRunoff(),
+)
+land = SoilCanopyModel{FT}(;
+    soilco2_type = soilco2_type,
+    soilco2_args = soilco2_args,
+    land_args = land_input,
+    soil_model_type = soil_model_type,
+    soil_args = soil_args,
+    canopy_component_types = canopy_component_types,
+    canopy_component_args = canopy_component_args,
+    canopy_model_args = canopy_model_args,
+)
+
+# Define explicit and implicit tendencies, and the jacobian
+exp_tendency! = make_exp_tendency(land)
+imp_tendency! = make_imp_tendency(land);
+jacobian! = make_jacobian(land);
+
+# Timestepping information
+N_hours = 12.5
+tf = Float64(t0 + N_hours * 3600.0)
+sim_time = round((tf - t0) / 3600, digits = 2) # simulation length in hours
+
+timestepper = CTS.ARS111()
+ode_algo = CTS.IMEXAlgorithm(
+    timestepper,
+    CTS.NewtonsMethod(
+        max_iters = 6,
+        update_j = CTS.UpdateEvery(CTS.NewNewtonIteration),
+    ),
+);
+
+# Set up simulation callbacks
+drivers = ClimaLand.get_drivers(land)
+updatefunc = ClimaLand.make_update_drivers(drivers)
+
+# Choose timesteps and set up arrays to store information
+ref_dt = 50.0
+dts = [ref_dt, 100.0, 200.0, 450.0, 900.0, 1800.0, 3600.0]
+
+ref_T = []
+mean_err = []
+p95_err = []
+p99_err = []
+sol_err = []
+T_states = []
+times = []
+ΔT = FT(0)
+for dt in dts
+    @info dt
+    # Initialize model and set initial conditions before each simulation
+    Y, p = set_initial_conditions(land, t0)
+    jac_kwargs = (;
+        jac_prototype = ClimaLand.ImplicitEquationJacobian(Y),
+        Wfact = jacobian!,
+    )
+
+    # Create callback for driver updates
+    saveat = vcat(Array(t0:(3 * 3600):tf), tf)
+    updateat = vcat(Array(t0:(3 * 3600):tf), tf)
+    driver_cb = ClimaLand.DriverUpdateCallback(updateat, updatefunc)
+
+    # Solve simulation
+    prob = SciMLBase.ODEProblem(
+        CTS.ClimaODEFunction(
+            T_exp! = exp_tendency!,
+            T_imp! = SciMLBase.ODEFunction(imp_tendency!; jac_kwargs...),
+            dss! = ClimaLand.dss!,
+        ),
+        Y,
+        (t0, tf),
+        p,
+    )
+    @time sol = SciMLBase.solve(
+        prob,
+        ode_algo;
+        dt = dt,
+        callback = driver_cb,
+        saveat = saveat,
+    )
+
+    # Save results for comparison
+    if dt == ref_dt
+        global ref_T =
+            [parent(sol.u[k].canopy.energy.T)[1] for k in 1:length(sol.t)]
+    else
+        T = [parent(sol.u[k].canopy.energy.T)[1] for k in 1:length(sol.t)]
+
+        ΔT = abs.(T .- ref_T)
+        push!(mean_err, FT(mean(ΔT)))
+        push!(p95_err, FT(percentile(ΔT, 95)))
+        push!(p99_err, FT(percentile(ΔT, 99)))
+        push!(sol_err, ΔT[end])
+        push!(T_states, T)
+        push!(times, sol.t)
+    end
+end
+
+savedir = joinpath(
+    pkgdir(ClimaLand),
+    "experiments",
+    "integrated",
+    "performance",
+    "integrated_timestep_test",
+);
+!ispath(savedir) && mkpath(savedir)
+
+# Create plot with statistics
+# Compare T state computed with small vs large dt
+fig = Figure()
+ax = Axis(
+    fig[1, 1],
+    xlabel = "Timestep (sec)",
+    ylabel = "Temperature (K)",
+    xscale = log10,
+    yscale = log10,
+    title = "Error in T over $(sim_time) hour sim, dts in [$(dts[1]), $(dts[end])]",
+)
+dts = dts[2:end]
+lines!(ax, dts, FT.(mean_err), label = "Mean Error")
+lines!(ax, dts, FT.(p95_err), label = "95th% Error")
+lines!(ax, dts, FT.(p99_err), label = "99th% Error")
+axislegend(ax)
+save(joinpath(savedir, "errors.png"), fig)
+
+# Create convergence plot
+fig2 = Figure()
+ax2 = Axis(
+    fig2[1, 1],
+    xlabel = "log(dt)",
+    ylabel = "log(|T[end] - T_ref[end]|)",
+    xscale = log10,
+    yscale = log10,
+    title = "Convergence of T; sim time = $(sim_time) hours, dts in [$(dts[1]), $(dts[end])]",
+)
+scatter!(ax2, dts, FT.(sol_err))
+lines!(ax2, dts, dts / 1e3)
+save(joinpath(savedir, "convergence.png"), fig2)
+
+# Plot T throughout full simulation for each run
+fig3 = Figure()
+ax3 = Axis(
+    fig3[1, 1],
+    xlabel = "time (hr)",
+    ylabel = "T (K)",
+    title = "T throughout simulation; length = $(sim_time) hours, dts in [$(dts[1]), $(dts[end])]",
+)
+times = times ./ 3600.0 # hours
+for i in 1:length(times)
+    lines!(ax3, times[i], T_states[i], label = "dt $(dts[i]) s")
+end
+axislegend(ax3, position = :rt)
+save(joinpath(savedir, "states.png"), fig3)
diff --git a/experiments/integrated/performance/profile_allocations.jl b/experiments/integrated/performance/profile_allocations.jl
index 2ec19f9192..1d0b9af6b5 100644
--- a/experiments/integrated/performance/profile_allocations.jl
+++ b/experiments/integrated/performance/profile_allocations.jl
@@ -332,15 +332,15 @@ imp_tendency! = make_imp_tendency(land);
 jacobian! = make_jacobian(land);
 
 # Set up timestepping and simulation callbacks
-dt = Float64(150)
+dt = Float64(900)
 tf = Float64(t0 + 100dt)
 saveat = Array(t0:dt:tf)
 updateat = Array(t0:dt:tf)
-timestepper = CTS.ARS343()
+timestepper = CTS.ARS111()
 ode_algo = CTS.IMEXAlgorithm(
     timestepper,
     CTS.NewtonsMethod(
-        max_iters = 1,
+        max_iters = 6,
         update_j = CTS.UpdateEvery(CTS.NewNewtonIteration),
     ),
 );
diff --git a/experiments/long_runs/land.jl b/experiments/long_runs/land.jl
index 060f0af4bf..643cb0c9ec 100644
--- a/experiments/long_runs/land.jl
+++ b/experiments/long_runs/land.jl
@@ -10,9 +10,9 @@
 # Soil depth: 50 m
 # Simulation duration: 365 d
 # Timestep: 900 s
-# Timestepper: ARS343
-# Fixed number of iterations: 1
-# Jacobian update: every new timestep
+# Timestepper: ARS111
+# Fixed number of iterations: 6
+# Jacobian update: every new Newton iteration
 # Atmos forcing update: every 3 hours
 import SciMLBase
 import ClimaComms
@@ -363,7 +363,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
     τ_NIR_leaf = FT(0.25)
 
     # Energy Balance model
-    ac_canopy = FT(2.5e4) # this will likely be 10x smaller!
+    ac_canopy = FT(2.5e3)
 
     #clm_data is used for g1 and vcmax maps
     clm_artifact_path = ClimaLand.Artifacts.clm_data_folder_path(; context)
@@ -631,7 +631,7 @@ end
 function setup_and_solve_problem(; greet = false)
 
     t0 = 0.0
-    tf = 60 * 60.0 * 24 * 365
+    tf = 60 * 60.0 * 24 * 7 # keep short until it runs! * 365
     Δt = 900.0
     nelements = (101, 15)
     if greet
@@ -648,7 +648,7 @@ function setup_and_solve_problem(; greet = false)
     ode_algo = CTS.IMEXAlgorithm(
         stepper,
         CTS.NewtonsMethod(
-            max_iters = 3,
+            max_iters = 6,
             update_j = CTS.UpdateEvery(CTS.NewNewtonIteration),
         ),
     )
diff --git a/experiments/standalone/Vegetation/timestep_test.jl b/experiments/standalone/Vegetation/timestep_test.jl
new file mode 100644
index 0000000000..072099249c
--- /dev/null
+++ b/experiments/standalone/Vegetation/timestep_test.jl
@@ -0,0 +1,336 @@
+# # Timestep test for standalone CanopyModel
+
+# This experiment runs the standalone canopy model for 20 days at a single point
+# and outputs useful plots related to the timestepping and stability of
+# the simulation using different timesteps. This was initially implemented to
+# test the behavior of stepping canopy temperature (T) implicitly.
+
+# The simulation is run multiple times: first, with a small reference timestep,
+# then with increasing timestep sizes. The goal here is to see which timesteps
+# are small enough to produce a stable simulation, and which timesteps are too
+# large to produce accurate results.
+
+# The runs with larger timesteps are compared with the reference timestep run
+# using multiple metrics, including:
+# - mean, 99%th, and 95th% error over the course of the simulation
+# - T at the end of the simulation (used to test convergence with each dt)
+# - T throughout the entire simulation
+
+# Note that the simulation is run for 20 days + 80 seconds. This somewhat
+# unusual length is necessary to test the convergence behavior of the solver.
+# Convergence behavior must be assessed when T is actively changing, rather
+# than when it is in equilibrium. T changes in the time period after a driver
+# update, which occurs every 3 hours (so exactly at 20 days). Running for
+# 20 days + 80 seconds allows us to look at results over a longer simulation,
+# as well as convergence behavior. See the discussion in
+# https://github.com/CliMA/ClimaLand.jl/pull/675 for more information.
+
+# Simulation Setup
+# Space: single point
+# Simulation duration: 20 days + 80 seconds (= 480.02 hours)
+# Timestep: variable, ranging from 0.6s to 3600s
+# Timestepper: ARS111
+# Fixed number of iterations: 6
+# Jacobian update: Every Newton iteration
+# Atmospheric and radiation driver updates: Every 3 hours
+
+import SciMLBase
+import ClimaComms
+@static pkgversion(ClimaComms) >= v"0.6" && ClimaComms.@import_required_backends
+using CairoMakie
+using Statistics
+import StatsBase: percentile
+using Dates
+using Insolation
+using StaticArrays
+
+# Load CliMA Packages and ClimaLand Modules:
+
+using ClimaCore
+import ClimaParams as CP
+import ClimaTimeSteppers as CTS
+using ClimaLand
+using ClimaLand.Domains: Point
+using ClimaLand.Canopy
+using ClimaLand.Canopy.PlantHydraulics
+import ClimaLand
+import ClimaLand.Parameters as LP
+const FT = Float64;
+earth_param_set = LP.LandParameters(FT);
+f_root_to_shoot = FT(3.5)
+plant_ν = FT(2.46e-4) # kg/m^2
+n_stem = Int64(1)
+n_leaf = Int64(1)
+h_leaf = FT(9.5)
+h_stem = FT(9)
+compartment_midpoints = [h_stem / 2, h_stem + h_leaf / 2]
+compartment_surfaces = [FT(0), h_stem, h_stem + h_leaf]
+land_domain = Point(; z_sfc = FT(0.0))
+include(
+    joinpath(pkgdir(ClimaLand), "experiments/integrated/fluxnet/data_tools.jl"),
+);
+time_offset = 7
+lat = FT(38.7441) # degree
+long = FT(-92.2000) # degree
+atmos_h = FT(32)
+site_ID = "US-MOz"
+data_link = "https://caltech.box.com/shared/static/7r0ci9pacsnwyo0o9c25mhhcjhsu6d72.csv"
+
+include(
+    joinpath(
+        pkgdir(ClimaLand),
+        "experiments/integrated/fluxnet/met_drivers_FLUXNET.jl",
+    ),
+);
+
+z0_m = FT(2)
+z0_b = FT(0.2)
+
+shared_params = SharedCanopyParameters{FT, typeof(earth_param_set)}(
+    z0_m,
+    z0_b,
+    earth_param_set,
+);
+ψ_soil0 = FT(0.0)
+
+soil_driver = PrescribedSoil(
+    FT;
+    root_depths = SVector{10, FT}(-(10:-1:1.0) ./ 10.0 * 2.0 .+ 0.2 / 2.0),
+    ψ = t -> ψ_soil0,
+    α_PAR = FT(0.2),
+    α_NIR = FT(0.4),
+    T = t -> 298.0,
+    ϵ = FT(0.99),
+);
+
+rt_params = TwoStreamParameters(
+    FT;
+    G_Function = ConstantGFunction(FT(0.5)),
+    α_PAR_leaf = FT(0.1),
+    α_NIR_leaf = FT(0.45),
+    τ_PAR_leaf = FT(0.05),
+    τ_NIR_leaf = FT(0.25),
+    Ω = FT(0.69),
+    λ_γ_PAR = FT(5e-7),
+    λ_γ_NIR = FT(1.65e-6),
+)
+rt_model = TwoStreamModel{FT}(rt_params);
+
+cond_params = MedlynConductanceParameters(FT; g1 = FT(141.0))
+stomatal_model = MedlynConductanceModel{FT}(cond_params);
+
+is_c3 = FT(1)
+photo_params = FarquharParameters(FT, is_c3; Vcmax25 = FT(5e-5))
+photosynthesis_model = FarquharModel{FT}(photo_params);
+
+AR_params = AutotrophicRespirationParameters(FT)
+AR_model = AutotrophicRespirationModel{FT}(AR_params);
+
+f_root_to_shoot = FT(3.5)
+SAI = FT(1.0)
+RAI = FT(3f_root_to_shoot)
+ai_parameterization = PrescribedSiteAreaIndex{FT}(LAIfunction, SAI, RAI)
+rooting_depth = FT(1.0)
+function root_distribution(z::T; rooting_depth = rooting_depth) where {T}
+    return T(1.0 / rooting_depth) * exp(z / T(rooting_depth)) # 1/m
+end
+
+K_sat_plant = FT(1.8e-6)
+ψ63 = FT(-4 / 0.0098)
+Weibull_param = FT(4)
+a = FT(0.05 * 0.0098)
+
+conductivity_model =
+    PlantHydraulics.Weibull{FT}(K_sat_plant, ψ63, Weibull_param)
+
+retention_model = PlantHydraulics.LinearRetentionCurve{FT}(a);
+
+ν = FT(0.7)
+S_s = FT(1e-2 * 0.0098)
+
+plant_hydraulics_ps = PlantHydraulics.PlantHydraulicsParameters(;
+    ai_parameterization = ai_parameterization,
+    ν = ν,
+    S_s = S_s,
+    root_distribution = root_distribution,
+    conductivity_model = conductivity_model,
+    retention_model = retention_model,
+);
+
+
+plant_hydraulics = PlantHydraulics.PlantHydraulicsModel{FT}(;
+    parameters = plant_hydraulics_ps,
+    n_stem = n_stem,
+    n_leaf = n_leaf,
+    compartment_surfaces = compartment_surfaces,
+    compartment_midpoints = compartment_midpoints,
+);
+ac_canopy = FT(1e3)
+energy_model = ClimaLand.Canopy.BigLeafEnergyModel{FT}(
+    BigLeafEnergyParameters{FT}(ac_canopy),
+)
+
+canopy = ClimaLand.Canopy.CanopyModel{FT}(;
+    parameters = shared_params,
+    domain = land_domain,
+    autotrophic_respiration = AR_model,
+    radiative_transfer = rt_model,
+    photosynthesis = photosynthesis_model,
+    conductance = stomatal_model,
+    energy = energy_model,
+    hydraulics = plant_hydraulics,
+    soil_driver = soil_driver,
+    atmos = atmos,
+    radiation = radiation,
+);
+
+exp_tendency! = make_exp_tendency(canopy)
+imp_tendency! = make_imp_tendency(canopy)
+jacobian! = make_jacobian(canopy)
+drivers = ClimaLand.get_drivers(canopy)
+updatefunc = ClimaLand.make_update_drivers(drivers)
+
+ψ_leaf_0 = FT(-2e5 / 9800)
+ψ_stem_0 = FT(-1e5 / 9800)
+S_l_ini =
+    inverse_water_retention_curve.(
+        retention_model,
+        [ψ_stem_0, ψ_leaf_0],
+        ν,
+        S_s,
+    )
+
+seconds_per_day = 3600 * 24.0
+t0 = 150seconds_per_day
+N_days = 20.0
+tf = t0 + N_days * seconds_per_day + 80
+sim_time = round((tf - t0) / 3600, digits = 2) # simulation length in hours
+set_initial_cache! = make_set_initial_cache(canopy)
+
+timestepper = CTS.ARS111();
+ode_algo = CTS.IMEXAlgorithm(
+    timestepper,
+    CTS.NewtonsMethod(
+        max_iters = 6,
+        update_j = CTS.UpdateEvery(CTS.NewNewtonIteration),
+    ),
+);
+
+ref_dt = 6.0
+dts = [ref_dt, 12.0, 48.0, 225.0, 450.0, 900.0, 1800.0, 3600.0]
+
+ref_T = []
+mean_err = []
+p95_err = []
+p99_err = []
+sol_err = []
+T_states = []
+times = []
+ΔT = FT(0)
+for dt in dts
+    @info dt
+
+    # Initialize model before each simulation
+    Y, p, coords = ClimaLand.initialize(canopy)
+    jac_kwargs = (;
+        jac_prototype = ClimaLand.ImplicitEquationJacobian(Y),
+        Wfact = jacobian!,
+    )
+
+    # Set initial conditions
+    Y.canopy.hydraulics.ϑ_l.:1 .= augmented_liquid_fraction.(ν, S_l_ini[1])
+    Y.canopy.hydraulics.ϑ_l.:2 .= augmented_liquid_fraction.(ν, S_l_ini[2])
+    evaluate!(Y.canopy.energy.T, atmos.T, t0)
+    set_initial_cache!(p, Y, t0)
+
+    # Create callback for driver updates
+    saveat = vcat(Array(t0:(3 * 3600):tf), tf)
+    updateat = vcat(Array(t0:(3 * 3600):tf), tf)
+    cb = ClimaLand.DriverUpdateCallback(updateat, updatefunc)
+
+    # Solve simulation
+    prob = SciMLBase.ODEProblem(
+        CTS.ClimaODEFunction(
+            T_exp! = exp_tendency!,
+            T_imp! = SciMLBase.ODEFunction(imp_tendency!; jac_kwargs...),
+            dss! = ClimaLand.dss!,
+        ),
+        Y,
+        (t0, tf),
+        p,
+    )
+    @time sol =
+        SciMLBase.solve(prob, ode_algo; dt = dt, callback = cb, saveat = saveat)
+
+    # Save results for comparison
+    if dt == ref_dt
+        global ref_T =
+            [parent(sol.u[k].canopy.energy.T)[1] for k in 1:length(sol.t)]
+    else
+        T = [parent(sol.u[k].canopy.energy.T)[1] for k in 1:length(sol.t)]
+
+        global ΔT = abs.(T .- ref_T)
+        push!(mean_err, FT(mean(ΔT)))
+        push!(p95_err, FT(percentile(ΔT, 95)))
+        push!(p99_err, FT(percentile(ΔT, 99)))
+        push!(sol_err, ΔT[end])
+        push!(T_states, T)
+        push!(times, sol.t)
+    end
+end
+
+savedir = joinpath(
+    pkgdir(ClimaLand),
+    "experiments",
+    "standalone",
+    "Vegetation",
+    "timestep_test",
+);
+!ispath(savedir) && mkpath(savedir)
+
+# Create plot with statistics
+# Compare T state computed with small vs large dt
+fig = Figure()
+ax = Axis(
+    fig[1, 1],
+    xlabel = "Timestep (minutes)",
+    ylabel = "Temperature (K)",
+    xscale = log10,
+    yscale = log10,
+    title = "Error in T over $(sim_time / 24) day sim, dts in [$(dts[1]), $(dts[end])]",
+)
+dts = dts[2:end] ./ 60
+lines!(ax, dts, FT.(mean_err), label = "Mean Error")
+lines!(ax, dts, FT.(p95_err), label = "95th% Error")
+lines!(ax, dts, FT.(p99_err), label = "99th% Error")
+axislegend(ax)
+save(joinpath(savedir, "errors.png"), fig)
+
+# Create convergence plot
+fig2 = Figure()
+ax2 = Axis(
+    fig2[1, 1],
+    xlabel = "log(dt)",
+    ylabel = "log(|T[end] - T_ref[end]|)",
+    xscale = log10,
+    yscale = log10,
+    title = "Convergence of T; sim time = $(sim_time / 24) days, dts in [$(dts[1]), $(dts[end])]",
+)
+scatter!(ax2, dts, FT.(sol_err))
+lines!(ax2, dts, dts / 1e3)
+save(joinpath(savedir, "convergence.png"), fig2)
+
+# Plot T throughout full simulation for each run
+fig3 = Figure()
+ax3 = Axis(
+    fig3[1, 1],
+    xlabel = "time (hr)",
+    ylabel = "T (K)",
+    title = "T throughout simulation; length = $(sim_time / 24) days, dts in [$(dts[1]), $(dts[end])]",
+)
+times = times ./ 3600.0 # hours
+for i in 1:length(times)
+    lines!(ax3, times[i], T_states[i], label = "dt $(dts[i]) s")
+end
+axislegend(ax3, position = :rt)
+save(joinpath(savedir, "states.png"), fig3)
diff --git a/src/shared_utilities/implicit_timestepping.jl b/src/shared_utilities/implicit_timestepping.jl
index 8ed224df1c..905ff081b4 100644
--- a/src/shared_utilities/implicit_timestepping.jl
+++ b/src/shared_utilities/implicit_timestepping.jl
@@ -1,5 +1,6 @@
 using ClimaCore.MatrixFields
 import ClimaCore.MatrixFields: @name, ⋅
+using ClimaCore: Spaces
 import LinearAlgebra
 import LinearAlgebra: I
 
@@ -107,53 +108,57 @@ to the `explicit_names` tuple.
 """
 function ImplicitEquationJacobian(Y::ClimaCore.Fields.FieldVector)
     FT = eltype(Y)
-    center_space = axes(Y.soil.ϑ_l)
-
-    # Construct 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)
-
     # Only add jacobian blocks for fields that are in Y for this model
-    is_in_Y(name) = MatrixFields.has_field(Y, name)
+    is_in_Y(var) = MatrixFields.has_field(Y, var)
 
     # Define the implicit and explicit variables of any model we use
-    implicit_names = (@name(soil.ϑ_l), @name(soil.ρe_int))
-    explicit_names = (
+    implicit_vars =
+        (@name(soil.ϑ_l), @name(soil.ρe_int), @name(canopy.energy.T))
+    explicit_vars = (
         @name(soilco2.C),
         @name(soil.θ_i),
         @name(canopy.hydraulics.ϑ_l),
-        @name(canopy.energy.T),
         @name(snow.S),
         @name(snow.U)
     )
 
     # Filter out the variables that are not in this model's state, `Y`
-    available_implicit_names =
-        MatrixFields.unrolled_filter(is_in_Y, implicit_names)
-    available_explicit_names =
-        MatrixFields.unrolled_filter(is_in_Y, explicit_names)
-
+    available_implicit_vars =
+        MatrixFields.unrolled_filter(is_in_Y, implicit_vars)
+    available_explicit_vars =
+        MatrixFields.unrolled_filter(is_in_Y, explicit_vars)
+
+    get_jac_type(
+        space::Union{
+            Spaces.FiniteDifferenceSpace,
+            Spaces.ExtrudedFiniteDifferenceSpace,
+        },
+        FT,
+    ) = MatrixFields.TridiagonalMatrixRow{FT}
+    get_jac_type(
+        space::Union{Spaces.PointSpace, Spaces.SpectralElementSpace2D},
+        FT,
+    ) = MatrixFields.DiagonalMatrixRow{FT}
+
+    get_j_field(space, FT) = zeros(get_jac_type(space, FT), space)
+
+    implicit_blocks = MatrixFields.unrolled_map(
+        var ->
+            (var, var) =>
+                get_j_field(axes(MatrixFields.get_field(Y, var)), FT),
+        available_implicit_vars,
+    )
     # For explicitly-stepped variables, use the negative identity matrix
     # Note: We have to use FT(-1) * I instead of -I because inv(-1) == -1.0,
     # which means that multiplying inv(-1) by a Float32 will yield a Float64.
-    identity_blocks = MatrixFields.unrolled_map(
-        name -> (name, name) => FT(-1) * I,
-        available_explicit_names,
+    explicit_blocks = MatrixFields.unrolled_map(
+        var -> (var, var) => FT(-1) * I,
+        available_explicit_vars,
     )
 
-    # For implicitly-stepped variables, use a tridiagonal matrix
-    tridiagonal_blocks = MatrixFields.unrolled_map(
-        name ->
-            (name, name) => Fields.Field(
-                tridiag_type,
-                axes(MatrixFields.get_field(Y, name)),
-            ),
-        available_implicit_names,
-    )
 
-    matrix = MatrixFields.FieldMatrix(identity_blocks..., tridiagonal_blocks...)
+
+    matrix = MatrixFields.FieldMatrix(implicit_blocks..., explicit_blocks...)
 
     # Set up block diagonal solver for block Jacobian
     alg = MatrixFields.BlockDiagonalSolve()
diff --git a/src/standalone/Vegetation/Canopy.jl b/src/standalone/Vegetation/Canopy.jl
index d81b495c25..3a50f72d52 100644
--- a/src/standalone/Vegetation/Canopy.jl
+++ b/src/standalone/Vegetation/Canopy.jl
@@ -3,11 +3,13 @@ using DocStringExtensions
 using Thermodynamics
 using ClimaLand
 using ClimaCore
+using ClimaCore.MatrixFields
+import ClimaCore.MatrixFields: @name, ⋅
+import LinearAlgebra: I
 using ClimaLand: AbstractRadiativeDrivers, AbstractAtmosphericDrivers
 import ..Parameters as LP
 
 import ClimaLand:
-    AbstractExpModel,
     name,
     prognostic_vars,
     prognostic_types,
@@ -20,6 +22,8 @@ import ClimaLand:
     make_update_boundary_fluxes,
     make_update_aux,
     make_compute_exp_tendency,
+    make_compute_imp_tendency,
+    make_compute_jacobian,
     get_drivers
 
 using ClimaLand.Domains: Point, Plane, SphericalSurface
@@ -54,7 +58,7 @@ struct SharedCanopyParameters{FT <: AbstractFloat, PSE}
 end
 
 """
-     CanopyModel{FT, AR, RM, PM, SM, PHM, EM, SM, A, R, S, PS, D} <: AbstractExpModel{FT}
+     CanopyModel{FT, AR, RM, PM, SM, PHM, EM, SM, A, R, S, PS, D} <: ClimaLand.AbstractImExModel{FT}
 
 The model struct for the canopy, which contains
 - the canopy model domain (a point for site-level simulations, or
@@ -98,7 +102,7 @@ treated differently.
 $(DocStringExtensions.FIELDS)
 """
 struct CanopyModel{FT, AR, RM, PM, SM, PHM, EM, SIFM, A, R, S, PS, D} <:
-       AbstractExpModel{FT}
+       ClimaLand.AbstractImExModel{FT}
     "Autotrophic respiration model, a canopy component model"
     autotrophic_respiration::AR
     "Radiative transfer model, a canopy component model"
@@ -627,13 +631,72 @@ function make_compute_exp_tendency(
     end
     return compute_exp_tendency!
 end
-function ClimaLand.get_drivers(model::CanopyModel)
-    return (model.atmos, model.radiation)
+
+"""
+    make_compute_imp_tendency(canopy::CanopyModel)
+
+Creates and returns the compute_imp_tendency! for the `CanopyModel`.
+"""
+function make_compute_imp_tendency(
+    canopy::CanopyModel{
+        FT,
+        <:AutotrophicRespirationModel,
+        <:Union{BeerLambertModel, TwoStreamModel},
+        <:Union{FarquharModel, OptimalityFarquharModel},
+        <:MedlynConductanceModel,
+        <:PlantHydraulicsModel,
+        <:Union{PrescribedCanopyTempModel, BigLeafEnergyModel},
+    },
+) where {FT}
+    components = canopy_components(canopy)
+    compute_imp_tendency_list = map(
+        x -> make_compute_imp_tendency(getproperty(canopy, x), canopy),
+        components,
+    )
+    function compute_imp_tendency!(dY, Y, p, t)
+        for f! in compute_imp_tendency_list
+            f!(dY, Y, p, t)
+        end
+
+    end
+    return compute_imp_tendency!
+end
+
+"""
+    ClimaLand.make_compute_jacobian(canopy::CanopyModel)
+
+Creates and returns the compute_jacobian! for the `CanopyModel`.
+"""
+function ClimaLand.make_compute_jacobian(
+    canopy::CanopyModel{
+        FT,
+        <:AutotrophicRespirationModel,
+        <:Union{BeerLambertModel, TwoStreamModel},
+        <:Union{FarquharModel, OptimalityFarquharModel},
+        <:MedlynConductanceModel,
+        <:PlantHydraulicsModel,
+        <:Union{PrescribedCanopyTempModel, BigLeafEnergyModel},
+    },
+) where {FT}
+    components = canopy_components(canopy)
+    update_jacobian_list = map(
+        x -> make_compute_jacobian(getproperty(canopy, x), canopy),
+        components,
+    )
+    function compute_jacobian!(W, Y, p, dtγ, t)
+        for f! in update_jacobian_list
+            f!(W, Y, p, dtγ, t)
+        end
+
+    end
+    return compute_jacobian!
 end
 
 
+function ClimaLand.get_drivers(model::CanopyModel)
+    return (model.atmos, model.radiation)
+end
 include("./canopy_boundary_fluxes.jl")
-
 #Make the canopy model broadcastable
 Base.broadcastable(C::CanopyModel) = tuple(C)
 end
diff --git a/src/standalone/Vegetation/canopy_boundary_fluxes.jl b/src/standalone/Vegetation/canopy_boundary_fluxes.jl
index 1e20b98c0a..49c0a5ca84 100644
--- a/src/standalone/Vegetation/canopy_boundary_fluxes.jl
+++ b/src/standalone/Vegetation/canopy_boundary_fluxes.jl
@@ -174,7 +174,8 @@ function canopy_boundary_fluxes!(
     shf .= canopy_tf.shf
     lhf .= canopy_tf.lhf
     r_ae .= canopy_tf.r_ae
-
+    p.canopy.energy.∂LHF∂qc .= canopy_tf.∂LHF∂qc
+    p.canopy.energy.∂SHF∂Tc .= canopy_tf.∂SHF∂Tc
     # Transpiration is per unit ground area, not leaf area (mult by LAI)
     fa.:($i_end) .= PlantHydraulics.transpiration_per_ground_area(
         canopy.hydraulics.transpiration,
@@ -335,7 +336,7 @@ function canopy_turbulent_fluxes_at_a_point(
     T_in::FT = Thermodynamics.air_temperature(thermo_params, ts_in)
     ΔT = T_in - T_sfc
     r_ae::FT = 1 / (conditions.Ch * SurfaceFluxes.windspeed(sc))
-    ρ_air::FT = Thermodynamics.air_density(thermo_params, ts_in)
+    ρ_sfc::FT = Thermodynamics.air_density(thermo_params, ts_sfc)
     ustar::FT = conditions.ustar
     r_b::FT = FT(1 / 0.01 * (ustar / 0.04)^(-1 / 2)) # CLM 5, tech note Equation 5.122
     leaf_r_b = r_b / LAI
@@ -343,7 +344,17 @@ function canopy_turbulent_fluxes_at_a_point(
     E0::FT = SurfaceFluxes.evaporation(surface_flux_params, sc, conditions.Ch)
     E = E0 * r_ae / (leaf_r_b + leaf_r_stomata + r_ae) # CLM 5, tech note Equation 5.101, and fig 5.2b, assuming all sunlit, f_wet = 0
     Ẽ = E / _ρ_liq
-    H = -ρ_air * cp_m * ΔT / (canopy_r_b + r_ae) # CLM 5, tech note Equation 5.88, setting H_v = H and solving to remove T_s
+    H = -ρ_sfc * cp_m * ΔT / (canopy_r_b + r_ae) # CLM 5, tech note Equation 5.88, setting H_v = H and solving to remove T_s
     LH = _LH_v0 * E
-    return (lhf = LH, shf = H, vapor_flux = Ẽ, r_ae = r_ae)
+    # We ignore ∂r_ae/∂T_sfc, ∂u*/∂T_sfc
+    ∂LHF∂qc = ρ_sfc * _LH_v0 / (leaf_r_b + leaf_r_stomata + r_ae) # Note that our estimate of ρ_sfc depends on T_sfc, but we ignore the derivative here.
+    ∂SHF∂Tc = ρ_sfc * cp_m / (canopy_r_b + r_ae) # Same approximation re: ρ_sfc(T_sfc)
+    return (
+        lhf = LH,
+        shf = H,
+        vapor_flux = Ẽ,
+        r_ae = r_ae,
+        ∂LHF∂qc = ∂LHF∂qc,
+        ∂SHF∂Tc = ∂SHF∂Tc,
+    )
 end
diff --git a/src/standalone/Vegetation/canopy_energy.jl b/src/standalone/Vegetation/canopy_energy.jl
index c3950225c7..96d4016941 100644
--- a/src/standalone/Vegetation/canopy_energy.jl
+++ b/src/standalone/Vegetation/canopy_energy.jl
@@ -10,12 +10,28 @@ abstract type AbstractCanopyEnergyModel{FT} <: AbstractCanopyComponent{FT} end
 ClimaLand.name(model::AbstractCanopyEnergyModel) = :energy
 
 
-ClimaLand.auxiliary_vars(model::AbstractCanopyEnergyModel) =
-    (:shf, :lhf, :fa_energy_roots, :r_ae)
+ClimaLand.auxiliary_vars(model::AbstractCanopyEnergyModel) = (
+    :shf,
+    :lhf,
+    :fa_energy_roots,
+    :r_ae,
+    :∂SHF∂Tc,
+    :∂LHF∂qc,
+    :∂LW_n∂Tc,
+    :∂qc∂Tc,
+)
 ClimaLand.auxiliary_types(model::AbstractCanopyEnergyModel{FT}) where {FT} =
-    (FT, FT, FT, FT)
-ClimaLand.auxiliary_domain_names(model::AbstractCanopyEnergyModel) =
-    (:surface, :surface, :surface, :surface)
+    (FT, FT, FT, FT, FT, FT, FT, FT)
+ClimaLand.auxiliary_domain_names(model::AbstractCanopyEnergyModel) = (
+    :surface,
+    :surface,
+    :surface,
+    :surface,
+    :surface,
+    :surface,
+    :surface,
+    :surface,
+)
 
 """
     PrescribedCanopyTempModel{FT} <: AbstractCanopyEnergyModel{FT}
@@ -87,11 +103,11 @@ where the canopy temperature is modeled prognostically.
 canopy_temperature(model::BigLeafEnergyModel, canopy, Y, p, t) =
     Y.canopy.energy.T
 
-function make_compute_exp_tendency(
+function make_compute_imp_tendency(
     model::BigLeafEnergyModel{FT},
     canopy,
 ) where {FT}
-    function compute_exp_tendency!(dY, Y, p, t)
+    function compute_imp_tendency!(dY, Y, p, t)
         area_index = p.canopy.hydraulics.area_index
         ac_canopy = model.parameters.ac_canopy
         # Energy Equation:
@@ -117,7 +133,7 @@ function make_compute_exp_tendency(
                 p.canopy.energy.lhf - p.canopy.energy.fa_energy_roots
             ) / c_per_ground_area
     end
-    return compute_exp_tendency!
+    return compute_imp_tendency!
 end
 
 """
@@ -154,3 +170,67 @@ function root_energy_flux_per_ground_area!(
 ) where {FT}
     fa_energy .= FT(0)
 end
+
+function ClimaLand.make_compute_jacobian(
+    model::BigLeafEnergyModel{FT},
+    canopy,
+) where {FT}
+    function compute_jacobian!(jacobian::ImplicitEquationJacobian, Y, p, dtγ, t)
+        (; matrix) = jacobian
+
+        # The derivative of the residual with respect to the prognostic variable
+        ∂Tres∂T = matrix[@name(canopy.energy.T), @name(canopy.energy.T)]
+        ∂LHF∂qc = p.canopy.energy.∂LHF∂qc
+        ∂SHF∂Tc = p.canopy.energy.∂SHF∂Tc
+        ∂LW_n∂Tc = p.canopy.energy.∂LW_n∂Tc
+        ∂qc∂Tc = p.canopy.energy.∂qc∂Tc
+        ϵ_c = p.canopy.radiative_transfer.ϵ
+        area_index = p.canopy.hydraulics.area_index
+        ac_canopy = model.parameters.ac_canopy
+        earth_param_set = canopy.parameters.earth_param_set
+        _T_freeze = LP.T_freeze(earth_param_set)
+        _σ = LP.Stefan(earth_param_set)
+        @. ∂LW_n∂Tc = -2 * 4 * _σ * ϵ_c * Y.canopy.energy.T^3 # ≈ ϵ_ground = 1
+        @. ∂qc∂Tc = partial_q_sat_partial_T_liq(
+            p.drivers.P,
+            Y.canopy.energy.T - _T_freeze,
+        )# use atmos air pressure as approximation for surface air pressure
+        @. ∂Tres∂T =
+            dtγ * MatrixFields.DiagonalMatrixRow(
+                (∂LW_n∂Tc - ∂SHF∂Tc - ∂LHF∂qc * ∂qc∂Tc) /
+                (ac_canopy * max(area_index.leaf + area_index.stem, eps(FT))),
+            ) - (I,)
+    end
+    return compute_jacobian!
+end
+
+"""
+    partial_q_sat_partial_T_liq(P::FT, T::FT) where {FT}
+
+Computes the quantity ∂q_sat∂T at temperature T and pressure P,
+over liquid water. The temperature must be in Celsius.
+
+Uses the polynomial approximation from Flatau et al. (1992).
+"""
+function partial_q_sat_partial_T_liq(P::FT, T::FT) where {FT}
+    esat = FT(
+        6.11213476e2 +
+        4.44007856e1 * T +
+        1.43064234 * T^2 +
+        2.64461437e-2 * T^3 +
+        3.05903558e-4 * T^4 +
+        1.96237241e-6 * T^5 +
+        8.92344772e-9 * T^6 - 3.73208410e-11 * T^7 + 2.09339997e-14 * T^8,
+    )
+    desatdT = FT(
+        4.44017302e1 +
+        2.86064092 * T +
+        7.94683137e-2 * T^2 +
+        1.21211669e-3 * T^3 +
+        1.03354611e-5 * T^4 +
+        4.04125005e-8 * T^5 - 7.88037859e-11 * T^6 - 1.14596802e-12 * T^7 +
+        3.81294516e-15 * T^8,
+    )
+
+    return FT(0.622) * P / (P - FT(0.378) * esat)^2 * desatdT
+end
diff --git a/src/standalone/Vegetation/component_models.jl b/src/standalone/Vegetation/component_models.jl
index 3697e0ce7c..ffda43199c 100644
--- a/src/standalone/Vegetation/component_models.jl
+++ b/src/standalone/Vegetation/component_models.jl
@@ -123,11 +123,50 @@ may be needed and passed in (via the `canopy` model itself).
 The right hand side for the entire canopy model can make use of
 these functions for the individual components.
 """
-function ClimaLand.make_compute_exp_tendency(::AbstractCanopyComponent, canopy)
-    function compute_exp_tendency!(dY, Y, p, t) end
+function ClimaLand.make_compute_exp_tendency(
+    component::AbstractCanopyComponent,
+    canopy,
+)
+    function compute_exp_tendency!(dY, Y, p, t)
+        vars = prognostic_vars(component)
+        dY_canopy = getproperty(dY, name(canopy))
+        if !isempty(vars)
+            getproperty(dY_canopy, name(component)) .= 0
+        end
+    end
     return compute_exp_tendency!
 end
 
+"""
+     ClimaLand.make_compute_imp_tendency(component::AbstractCanopyComponent, canopy)
+
+Creates the compute_imp_tendency!(dY,Y,p,t) function for the canopy `component`.
+
+Since component models are not standalone models, other information
+may be needed and passed in (via the `canopy` model itself).
+The right hand side for the entire canopy model can make use of
+these functions for the individual components.
+"""
+function ClimaLand.make_compute_imp_tendency(
+    component::AbstractCanopyComponent,
+    canopy,
+)
+    function compute_imp_tendency!(dY, Y, p, t)
+        vars = prognostic_vars(component)
+        dY_canopy = getproperty(dY, name(canopy))
+        if !isempty(vars)
+            getproperty(dY_canopy, name(component)) .= 0
+        end
+
+    end
+    return compute_imp_tendency!
+end
+
+
+function make_compute_jacobian(component::AbstractCanopyComponent, canopy)
+    function compute_jacobian!(W, Y, p, dtγ, t) end
+    return compute_jacobian!
+end
 
 """
      set_canopy_prescribed_field!(component::AbstractCanopyComponent,
diff --git a/test/shared_utilities/variable_types.jl b/test/shared_utilities/variable_types.jl
index 06b2b2e9ee..53fa98f1d4 100644
--- a/test/shared_utilities/variable_types.jl
+++ b/test/shared_utilities/variable_types.jl
@@ -78,8 +78,14 @@ end
 
 @testset "Default canopy component" begin
     FT = Float64
+    struct DefaultCanopy{FT} <: AbstractImExModel{FT} end
+    ClimaLand.name(::DefaultCanopy) where {FT} = :canopy
+    dc = DefaultCanopy{FT}()
+
     struct DefaultCanopyComponent{FT} <: AbstractCanopyComponent{FT} end
+    ClimaLand.name(::DefaultCanopyComponent) where {FT} = :component
     dcc = DefaultCanopyComponent{FT}()
+
     @test ClimaLand.prognostic_vars(dcc) == ()
     @test ClimaLand.prognostic_types(dcc) == ()
     @test ClimaLand.prognostic_domain_names(dcc) == ()
@@ -88,10 +94,9 @@ end
     @test ClimaLand.auxiliary_types(dcc) == ()
     @test ClimaLand.auxiliary_domain_names(dcc) == ()
 
-    x = [0, 1, 2, 3]
-    dcc_compute_exp_tendency! = make_compute_exp_tendency(dcc, nothing)
+    x = [(canopy = (component = [1],),), 1, 2, 3]
+    dcc_compute_exp_tendency! = make_compute_exp_tendency(dcc, dc)
     @test dcc_compute_exp_tendency!(x[1], x[2], x[3], x[4]) == nothing
-    @test x == [0, 1, 2, 3]
 
     @test_throws MethodError make_compute_imp_tendency(dcc)
 end
diff --git a/test/standalone/Vegetation/canopy_model.jl b/test/standalone/Vegetation/canopy_model.jl
index b00bec43e6..60eca1e96e 100644
--- a/test/standalone/Vegetation/canopy_model.jl
+++ b/test/standalone/Vegetation/canopy_model.jl
@@ -13,7 +13,7 @@ using ClimaLand.Canopy
 using ClimaLand.Canopy.PlantHydraulics
 using ClimaLand.Domains: Point
 import Insolation
-
+using ClimaCore.MatrixFields: @name
 import ClimaLand
 import ClimaLand.Parameters as LP
 import ClimaParams
@@ -42,7 +42,6 @@ import ClimaParams
             RTparams = BeerLambertParameters(FT)
             photosynthesis_params = FarquharParameters(FT, is_c3; Vcmax25)
             stomatal_g_params = MedlynConductanceParameters(FT; g1)
-
             AR_model = AutotrophicRespirationModel{FT}(AR_params)
             stomatal_model = MedlynConductanceModel{FT}(stomatal_g_params)
             photosynthesis_model = FarquharModel{FT}(photosynthesis_params)
@@ -247,6 +246,14 @@ import ClimaParams
             parent(Y.canopy.hydraulics.ϑ_l) .= plant_ν
             set_initial_cache! = make_set_initial_cache(canopy)
             exp_tendency! = make_exp_tendency(canopy)
+            imp_tendency! = ClimaLand.make_imp_tendency(canopy)
+            jacobian! = ClimaLand.make_jacobian(canopy)
+            # set up jacobian info
+            jac_kwargs = (;
+                jac_prototype = ImplicitEquationJacobian(Y),
+                Wfact = jacobian!,
+            )
+
             t0 = FT(0.0)
             dY = similar(Y)
             set_initial_cache!(p, Y, t0)
@@ -662,6 +669,7 @@ end
                 θs = zenith_angle,
             )
 
+
             # Plant Hydraulics
             RAI = FT(1)
             SAI = FT(0)
@@ -787,6 +795,14 @@ end
             t0 = FT(0.0)
             set_initial_cache!(p, Y, t0)
             exp_tendency! = make_exp_tendency(canopy)
+            imp_tendency! = ClimaLand.make_imp_tendency(canopy)
+            jacobian! = ClimaLand.make_jacobian(canopy)
+            # set up jacobian info
+            jac_kwargs = (;
+                jac_prototype = ImplicitEquationJacobian(Y),
+                Wfact = jacobian!,
+            )
+
             dY = similar(Y)
             exp_tendency!(dY, Y, p, t0)
             turb_fluxes =
@@ -819,6 +835,278 @@ end
     end
 end
 
+@testset "Jacobian for Temperature" begin
+    for FT in (Float32, Float64)
+        domain = Point(; z_sfc = FT(0.0))
+
+        g1 = FT(790)
+        Vcmax25 = FT(9e-5)
+        is_c3 = FT(1)
+        RTparams = BeerLambertParameters(FT)
+        photosynthesis_params = FarquharParameters(FT, is_c3; Vcmax25)
+        stomatal_g_params = MedlynConductanceParameters(FT; g1)
+
+        stomatal_model = MedlynConductanceModel{FT}(stomatal_g_params)
+        photosynthesis_model = FarquharModel{FT}(photosynthesis_params)
+        rt_model = BeerLambertModel{FT}(RTparams)
+        energy_model = BigLeafEnergyModel{FT}(BigLeafEnergyParameters{FT}())
+        earth_param_set = LP.LandParameters(FT)
+        LAI = FT(8.0) # m2 [leaf] m-2 [ground]
+        z_0m = FT(2.0) # m, Roughness length for momentum - value from tall forest ChatGPT
+        z_0b = FT(0.1) # m, Roughness length for scalars - value from tall forest ChatGPT
+        h_int = FT(30.0) # m, "where measurements would be taken at a typical flux tower of a 20m canopy"
+        shared_params = SharedCanopyParameters{FT, typeof(earth_param_set)}(
+            z_0m,
+            z_0b,
+            earth_param_set,
+        )
+        lat = FT(0.0) # degree
+        long = FT(-180) # degree
+
+        function zenith_angle(
+            t,
+            start_date;
+            latitude = lat,
+            longitude = long,
+            insol_params = earth_param_set.insol_params,
+        )
+            current_datetime = start_date + Dates.Second(round(t))
+            d, δ, η_UTC =
+                FT.(
+                    Insolation.helper_instantaneous_zenith_angle(
+                        current_datetime,
+                        start_date,
+                        insol_params,
+                    )
+                )
+            return Insolation.instantaneous_zenith_angle(
+                d,
+                δ,
+                η_UTC,
+                longitude,
+                latitude,
+            )[1]
+        end
+
+        function shortwave_radiation(
+            t;
+            latitude = lat,
+            longitude = long,
+            insol_params = earth_param_set.insol_params,
+        )
+            return 1000 # W/m^2
+        end
+
+        function longwave_radiation(t)
+            return 200 # W/m^2
+        end
+
+        u_atmos = t -> 10 #m.s-1
+
+        liquid_precip = (t) -> 0 # m
+        snow_precip = (t) -> 0 # m
+        T_atmos = t -> 290 # Kelvin
+        q_atmos = t -> 0.001 # kg/kg
+        P_atmos = t -> 1e5 # Pa
+        h_atmos = h_int # m
+        c_atmos = (t) -> 4.11e-4 # mol/mol
+        start_date = DateTime(2005)
+        atmos = PrescribedAtmosphere(
+            TimeVaryingInput(liquid_precip),
+            TimeVaryingInput(snow_precip),
+            TimeVaryingInput(T_atmos),
+            TimeVaryingInput(u_atmos),
+            TimeVaryingInput(q_atmos),
+            TimeVaryingInput(P_atmos),
+            start_date,
+            h_atmos,
+            earth_param_set;
+            c_co2 = TimeVaryingInput(c_atmos),
+        )
+        radiation = PrescribedRadiativeFluxes(
+            FT,
+            TimeVaryingInput(shortwave_radiation),
+            TimeVaryingInput(longwave_radiation),
+            start_date;
+            θs = zenith_angle,
+        )
+
+        # Plant Hydraulics
+        RAI = FT(1)
+        SAI = FT(0)
+        lai_fun = t -> LAI
+        ai_parameterization = PlantHydraulics.PrescribedSiteAreaIndex{FT}(
+            TimeVaryingInput(lai_fun),
+            SAI,
+            RAI,
+        )
+        K_sat_plant = FT(1.8e-8) # m/s
+        ψ63 = FT(-4 / 0.0098) # / MPa to m, Holtzman's original parameter value
+        Weibull_param = FT(4) # unitless, Holtzman's original c param value
+        a = FT(0.05 * 0.0098) # Holtzman's original parameter for the bulk modulus of elasticity
+        plant_ν = FT(0.7) # m3/m3
+        plant_S_s = FT(1e-2 * 0.0098) # m3/m3/MPa to m3/m3/m
+        conductivity_model =
+            PlantHydraulics.Weibull{FT}(K_sat_plant, ψ63, Weibull_param)
+        retention_model = PlantHydraulics.LinearRetentionCurve{FT}(a)
+        root_depths = SVector{10, FT}(-(10:-1:1.0) ./ 10.0 * 2.0 .+ 0.2 / 2.0) # 1st element is the deepest root depth
+        function root_distribution(z::T) where {T}
+            return T(1.0 / 0.5) * exp(z / T(0.5)) # (1/m)
+        end
+        param_set = PlantHydraulics.PlantHydraulicsParameters(;
+            ai_parameterization = ai_parameterization,
+            ν = plant_ν,
+            S_s = plant_S_s,
+            root_distribution = root_distribution,
+            conductivity_model = conductivity_model,
+            retention_model = retention_model,
+        )
+        Δz = FT(1.0) # height of compartments
+        n_stem = Int64(0) # number of stem elements
+        n_leaf = Int64(1) # number of leaf elements
+        compartment_centers =
+            FT.(
+                Vector(
+                    range(
+                        start = Δz / 2,
+                        step = Δz,
+                        stop = Δz * (n_stem + n_leaf) - (Δz / 2),
+                    ),
+                ),
+            )
+        compartment_faces =
+            FT.(
+                Vector(
+                    range(
+                        start = 0.0,
+                        step = Δz,
+                        stop = Δz * (n_stem + n_leaf),
+                    ),
+                )
+            )
+
+        ψ_soil0 = FT(0.0)
+        T_soil0 = FT(290)
+        soil_driver = PrescribedSoil(
+            root_depths,
+            (t) -> ψ_soil0,
+            (t) -> T_soil0,
+            FT(0.2),
+            FT(0.4),
+            FT(0.98),
+        )
+
+        plant_hydraulics = PlantHydraulics.PlantHydraulicsModel{FT}(;
+            parameters = param_set,
+            n_stem = n_stem,
+            n_leaf = n_leaf,
+            compartment_surfaces = compartment_faces,
+            compartment_midpoints = compartment_centers,
+        )
+        autotrophic_parameters = AutotrophicRespirationParameters(FT)
+        autotrophic_respiration_model =
+            AutotrophicRespirationModel{FT}(autotrophic_parameters)
+
+        canopy = ClimaLand.Canopy.CanopyModel{FT}(;
+            parameters = shared_params,
+            domain = domain,
+            radiative_transfer = rt_model,
+            photosynthesis = photosynthesis_model,
+            conductance = stomatal_model,
+            autotrophic_respiration = autotrophic_respiration_model,
+            energy = energy_model,
+            hydraulics = plant_hydraulics,
+            soil_driver = soil_driver,
+            atmos = atmos,
+            radiation = radiation,
+        )
+
+        Y, p, coords = ClimaLand.initialize(canopy)
+
+        Y.canopy.hydraulics .= plant_ν
+        Y.canopy.energy.T = FT(289)
+
+        set_initial_cache! = make_set_initial_cache(canopy)
+        t0 = FT(0.0)
+        imp_tendency! = ClimaLand.make_imp_tendency(canopy)
+        jacobian! = ClimaLand.make_jacobian(canopy)
+
+        set_initial_cache!(p, Y, t0)
+        T_sfc = ClimaLand.surface_temperature(canopy, Y, p, t0)
+        ρ_sfc =
+            ClimaLand.surface_air_density(canopy.atmos, canopy, Y, p, t0, T_sfc)
+        q_sfc = ClimaLand.surface_specific_humidity(canopy, Y, p, T_sfc, ρ_sfc)
+        dY = similar(Y)
+        imp_tendency!(dY, Y, p, t0)
+        jac = ImplicitEquationJacobian(Y)
+        jacobian!(jac, Y, p, FT(1), t0)
+        jac_value =
+            parent(jac.matrix[@name(canopy.energy.T), @name(canopy.energy.T)])
+        ΔT = FT(0.01)
+
+        Y_2 = deepcopy(Y)
+        Y_2.canopy.energy.T = FT(289 + ΔT)
+        p_2 = deepcopy(p)
+        set_initial_cache!(p_2, Y_2, t0)
+        T_sfc2 = ClimaLand.surface_temperature(canopy, Y_2, p_2, t0)
+        ρ_sfc2 = ClimaLand.surface_air_density(
+            canopy.atmos,
+            canopy,
+            Y_2,
+            p_2,
+            t0,
+            T_sfc2,
+        )
+        q_sfc2 = ClimaLand.surface_specific_humidity(
+            canopy,
+            Y_2,
+            p_2,
+            T_sfc2,
+            ρ_sfc2,
+        )
+        dY_2 = similar(Y_2)
+        imp_tendency!(dY_2, Y_2, p_2, t0)
+
+        finitediff_LW =
+            (
+                p_2.canopy.radiative_transfer.LW_n .-
+                p.canopy.radiative_transfer.LW_n
+            ) ./ ΔT
+        estimated_LW = p.canopy.energy.∂LW_n∂Tc
+        @test parent(abs.(finitediff_LW .- estimated_LW) ./ finitediff_LW)[1] <
+              0.01
+
+        finitediff_SHF = (p_2.canopy.energy.shf .- p.canopy.energy.shf) ./ ΔT
+        estimated_SHF = p.canopy.energy.∂SHF∂Tc
+        @test parent(abs.(finitediff_SHF .- estimated_SHF) ./ finitediff_SHF)[1] <
+              0.05
+
+        finitediff_LHF = (p_2.canopy.energy.lhf .- p.canopy.energy.lhf) ./ ΔT
+        estimated_LHF = p.canopy.energy.∂LHF∂qc .* p.canopy.energy.∂qc∂Tc
+        @test parent(abs.(finitediff_LHF .- estimated_LHF) ./ finitediff_LHF)[1] <
+              0.5
+
+        # Error in `q` derivative is large
+        finitediff_q = (q_sfc2 .- q_sfc) ./ ΔT
+        estimated_q = p.canopy.energy.∂qc∂Tc
+        @test parent(abs.(finitediff_q .- estimated_q) ./ finitediff_q)[1] <
+              0.25
+
+        # Im not sure why this is not smaller! There must be an error in ∂LHF∂qc also.
+        estimated_LHF_with_correct_q = p.canopy.energy.∂LHF∂qc .* finitediff_q
+        @test parent(
+            abs.(finitediff_LHF .- estimated_LHF_with_correct_q) ./
+            finitediff_LHF,
+        )[1] < 0.5
+
+        # Recall jac = ∂Ṫ∂T - 1 [dtγ = 1]
+        ∂Ṫ∂T = jac_value .+ 1
+        @test (abs.(
+            parent((dY_2.canopy.energy.T .- dY.canopy.energy.T) ./ ΔT) - ∂Ṫ∂T
+        ) / ∂Ṫ∂T)[1] < 0.25 # Error propagates here from ∂LHF∂T
+    end
+end
+
 
 @testset "Zero LAI;" begin
     for FT in (Float32, Float64)