diff --git a/docs/src/APIs/canopy/PlantHydraulics.md b/docs/src/APIs/canopy/PlantHydraulics.md
index 4e7ab468d3..cabe6ce005 100644
--- a/docs/src/APIs/canopy/PlantHydraulics.md
+++ b/docs/src/APIs/canopy/PlantHydraulics.md
@@ -17,7 +17,7 @@ ClimaLand.PlantHydraulics.augmented_liquid_fraction
 ClimaLand.PlantHydraulics.water_retention_curve
 ClimaLand.PlantHydraulics.inverse_water_retention_curve
 ClimaLand.PlantHydraulics.root_water_flux_per_ground_area!
-ClimaLand.PlantHydraulics.flux
+ClimaLand.PlantHydraulics.water_flux
 ClimaLand.PlantHydraulics.hydraulic_conductivity
 ```
 
diff --git a/experiments/long_runs/land.jl b/experiments/long_runs/land.jl
index 466d75200d..aa6e2024e6 100644
--- a/experiments/long_runs/land.jl
+++ b/experiments/long_runs/land.jl
@@ -28,7 +28,8 @@ import ClimaAnalysis
 import ClimaAnalysis.Visualize as viz
 import ClimaUtilities
 
-import ClimaUtilities.TimeVaryingInputs: TimeVaryingInput
+import ClimaUtilities.TimeVaryingInputs:
+    TimeVaryingInput, LinearInterpolation, PeriodicCalendar
 import ClimaUtilities.SpaceVaryingInputs: SpaceVaryingInput
 import ClimaUtilities.Regridders: InterpolationsRegridder
 import ClimaUtilities.ClimaArtifacts: @clima_artifact
@@ -47,7 +48,7 @@ using Dates
 import NCDatasets
 
 const FT = Float64;
-
+time_interpolation_method = LinearInterpolation(PeriodicCalendar())
 regridder_type = :InterpolationsRegridder
 context = ClimaComms.context()
 device = ClimaComms.device()
@@ -83,6 +84,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         reference_date = start_date,
         regridder_type,
         file_reader_kwargs = (; preprocess_func = (data) -> -data / 3600,),
+        method = time_interpolation_method,
     )
 
     snow_precip = TimeVaryingInput(
@@ -92,6 +94,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         reference_date = start_date,
         regridder_type,
         file_reader_kwargs = (; preprocess_func = (data) -> -data / 3600,),
+        method = time_interpolation_method,
     )
 
     u_atmos = TimeVaryingInput(
@@ -100,6 +103,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         surface_space;
         reference_date = start_date,
         regridder_type,
+        method = time_interpolation_method,
     )
     q_atmos = TimeVaryingInput(
         joinpath(era5_artifact_path, "era5_2021_0.9x1.25_clima.nc"),
@@ -107,6 +111,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         surface_space;
         reference_date = start_date,
         regridder_type,
+        method = time_interpolation_method,
     )
     P_atmos = TimeVaryingInput(
         joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"),
@@ -114,6 +119,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         surface_space;
         reference_date = start_date,
         regridder_type,
+        method = time_interpolation_method,
     )
 
     T_atmos = TimeVaryingInput(
@@ -122,6 +128,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         surface_space;
         reference_date = start_date,
         regridder_type,
+        method = time_interpolation_method,
     )
     h_atmos = FT(10)
 
@@ -145,6 +152,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         reference_date = start_date,
         regridder_type,
         file_reader_kwargs = (; preprocess_func = (data) -> data / 3600,),
+        method = time_interpolation_method,
     )
     LW_d = TimeVaryingInput(
         joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"),
@@ -153,6 +161,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         reference_date = start_date,
         regridder_type,
         file_reader_kwargs = (; preprocess_func = (data) -> data / 3600,),
+        method = time_interpolation_method,
     )
 
     function zenith_angle(
@@ -489,6 +498,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         file_reader_kwargs = (;
             preprocess_func = (data) -> data > 0.05 ? data : 0.0,
         ),
+        method = time_interpolation_method,
     )
     ai_parameterization =
         Canopy.PrescribedSiteAreaIndex{FT}(LAIfunction, SAI, RAI)
@@ -613,6 +623,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         start_date;
         output_writer = nc_writer,
         output_vars = :short,
+        average_period = :monthly,
     )
 
     diagnostic_handler =
@@ -627,7 +638,7 @@ end
 function setup_and_solve_problem(; greet = false)
 
     t0 = 0.0
-    tf = 60 * 60.0 * 24 * 14
+    tf = 60 * 60.0 * 24 * 365
     Δt = 900.0
     nelements = (101, 15)
     if greet
diff --git a/experiments/long_runs/soil.jl b/experiments/long_runs/soil.jl
index c603ebfd30..e21413ec6b 100644
--- a/experiments/long_runs/soil.jl
+++ b/experiments/long_runs/soil.jl
@@ -28,7 +28,8 @@ using ClimaAnalysis
 import ClimaAnalysis.Visualize as viz
 using ClimaUtilities
 
-import ClimaUtilities.TimeVaryingInputs: TimeVaryingInput
+import ClimaUtilities.TimeVaryingInputs:
+    TimeVaryingInput, LinearInterpolation, PeriodicCalendar
 import ClimaUtilities.SpaceVaryingInputs: SpaceVaryingInput
 import ClimaUtilities.Regridders: InterpolationsRegridder
 import ClimaUtilities.ClimaArtifacts: @clima_artifact
@@ -45,7 +46,7 @@ using Dates
 import NCDatasets
 
 const FT = Float64;
-
+time_interpolation_method = LinearInterpolation(PeriodicCalendar())
 regridder_type = :InterpolationsRegridder
 context = ClimaComms.context()
 device = ClimaComms.device()
@@ -81,6 +82,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         reference_date = start_date,
         regridder_type,
         file_reader_kwargs = (; preprocess_func = (data) -> -data / 3600,),
+        method = time_interpolation_method,
     )
 
     snow_precip = TimeVaryingInput(
@@ -90,6 +92,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         reference_date = start_date,
         regridder_type,
         file_reader_kwargs = (; preprocess_func = (data) -> -data / 3600,),
+        method = time_interpolation_method,
     )
 
     u_atmos = TimeVaryingInput(
@@ -98,6 +101,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         surface_space;
         reference_date = start_date,
         regridder_type,
+        method = time_interpolation_method,
     )
     q_atmos = TimeVaryingInput(
         joinpath(era5_artifact_path, "era5_2021_0.9x1.25_clima.nc"),
@@ -105,6 +109,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         surface_space;
         reference_date = start_date,
         regridder_type,
+        method = time_interpolation_method,
     )
     P_atmos = TimeVaryingInput(
         joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"),
@@ -112,6 +117,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         surface_space;
         reference_date = start_date,
         regridder_type,
+        method = time_interpolation_method,
     )
 
     T_atmos = TimeVaryingInput(
@@ -120,6 +126,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         surface_space;
         reference_date = start_date,
         regridder_type,
+        method = time_interpolation_method,
     )
     h_atmos = FT(10)
 
@@ -143,6 +150,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         reference_date = start_date,
         regridder_type,
         file_reader_kwargs = (; preprocess_func = (data) -> data / 3600,),
+        method = time_interpolation_method,
     )
     LW_d = TimeVaryingInput(
         joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"),
@@ -151,6 +159,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         reference_date = start_date,
         regridder_type,
         file_reader_kwargs = (; preprocess_func = (data) -> data / 3600,),
+        method = time_interpolation_method,
     )
 
     function zenith_angle(
@@ -430,7 +439,7 @@ end
 function setup_and_solve_problem(; greet = false)
 
     t0 = 0.0
-    tf = 60 * 60.0 * 24 * 340
+    tf = 60 * 60.0 * 24 * 365
     Δt = 900.0
     nelements = (101, 15)
     if greet
diff --git a/src/integrated/soil_canopy_model.jl b/src/integrated/soil_canopy_model.jl
index b7c40688bc..347f269183 100644
--- a/src/integrated/soil_canopy_model.jl
+++ b/src/integrated/soil_canopy_model.jl
@@ -263,18 +263,22 @@ function make_update_boundary_fluxes(
             area_index,
             land.canopy.hydraulics.compartment_labels[1],
         )
+        # Note that we model the flux between each soil layer and the canopy as:
+        # Flux = -K_eff x [(ψ_canopy - ψ_soil)/(z_canopy - z_soil) + 1], where
+        # K_eff = K_soil K_canopy /(K_canopy + K_soil)
+
+        # Note that in `PrescribedSoil` mode, we compute the flux using K_soil = K_plant(ψ_soil)
+        # and K_canopy = K_plant(ψ_canopy). In `PrognosticSoil` mode here, we compute the flux using
+        # K_soil = K_soil(ψ_soil) and K_canopy = K_plant(ψ_canopy).
 
         @. p.root_extraction =
             above_ground_area_index *
-            PlantHydraulics.flux(
+            PlantHydraulics.water_flux(
                 z,
                 land.canopy.hydraulics.compartment_midpoints[1],
                 p.soil.ψ,
                 p.canopy.hydraulics.ψ.:1,
-                PlantHydraulics.hydraulic_conductivity(
-                    conductivity_model,
-                    p.soil.ψ,
-                ),
+                p.soil.K,
                 PlantHydraulics.hydraulic_conductivity(
                     conductivity_model,
                     p.canopy.hydraulics.ψ.:1,
diff --git a/src/standalone/Vegetation/Canopy.jl b/src/standalone/Vegetation/Canopy.jl
index 1b0db859f5..d81b495c25 100644
--- a/src/standalone/Vegetation/Canopy.jl
+++ b/src/standalone/Vegetation/Canopy.jl
@@ -524,7 +524,7 @@ function ClimaLand.make_update_aux(
             # Compute the flux*area between the current compartment `i`
             # and the compartment above.
             @. fa.:($$i) =
-                PlantHydraulics.flux(
+                PlantHydraulics.water_flux(
                     hydraulics.compartment_midpoints[i],
                     hydraulics.compartment_midpoints[ip1],
                     ψ.:($$i),
diff --git a/src/standalone/Vegetation/PlantHydraulics.jl b/src/standalone/Vegetation/PlantHydraulics.jl
index 26727dc5d8..47ae13a88a 100644
--- a/src/standalone/Vegetation/PlantHydraulics.jl
+++ b/src/standalone/Vegetation/PlantHydraulics.jl
@@ -23,7 +23,7 @@ import ClimaLand:
     name
 export PlantHydraulicsModel,
     AbstractPlantHydraulicsModel,
-    flux,
+    water_flux,
     effective_saturation,
     augmented_liquid_fraction,
     water_retention_curve,
@@ -328,7 +328,7 @@ end
 
 
 """
-    flux(
+    water_flux(
         z1,
         z2,
         ψ1,
@@ -337,18 +337,25 @@ end
         K2,
     ) where {FT}
 
-Computes the water flux given the absolute potential (pressure/(ρg))
- at the center of the two compartments z1 and z2,
-and the conductivity along
-the flow path between these two points.
+Computes the water flux given the absolute potential ψ (pressure/(ρg))
+ and the conductivity K (m/s) at the center of the two layers
+with midpoints z1 and z2.
 
-We currently assuming an arithmetic
-mean for mean K_sat between the two points (Bonan, 2019; Zhu, 2008)
-to take into account the change in K_sat halfway
-between z1 and z2; this is incorrect for compartments of differing sizes.
+We currently assuming a harmonic
+mean for effective conducticity between the two layers
+(see CLM Technical Documentation).
+
+To account for different path lengths in the two compartments Δz1 and
+Δz2, we would require the following conductance k (1/s)
+k_eff = K1/Δz1*K2/Δz2/(K1/Δz1+K2/Δz2) 
+and a water flux of
+F = -k_eff * (ψ1 +z1 - ψ2 - z2) (m/s).
+
+This currently assumes the path lengths are equal.
 """
-function flux(z1, z2, ψ1, ψ2, K1, K2)
-    flux = -(K1 + K2) / 2 * ((ψ2 - ψ1) / (z2 - z1) + 1)
+function water_flux(z1::FT, z2::FT, ψ1::FT, ψ2::FT, K1::FT, K2::FT) where {FT}
+    K_eff = K1 * K2 / max(K1 + K2, eps(FT))
+    flux = -K_eff * ((ψ2 - ψ1) / (z2 - z1) + 1)
     return flux # (m/s)
 end
 
@@ -580,8 +587,16 @@ end
     ) where {FT}
 
 A method which computes the water flux between the soil and the stem, via the roots,
-and multiplied by the RAI, in the case of a model running without an integrated
-soil model.
+and multiplied by the RAI, in the case of a model running without a prognostic
+soil model:
+
+Flux = -K_eff x [(ψ_stem - ψ_soil)/(z_stem - z_soil) + 1], where
+K_eff = K_soil K_stem /(K_stem + K_soil)
+
+Note that in `PrescribedSoil` mode, we compute the flux using K_soil = K_plant(ψ_soil)
+and K_stem = K_plant(ψ_stem). In `PrognosticSoil` mode, we compute the flux using
+K_soil = K_soil(ψ_soil) and K_stem = K_plant(ψ_stem). The latter is a better model, but
+our `PrescribedSoil` struct does not store K_soil, only ψ_soil.
 
 The returned flux is per unit ground area. This assumes that the stem compartment
 is the first element of `Y.canopy.hydraulics.ϑ_l`.
@@ -610,7 +625,7 @@ function root_water_flux_per_ground_area!(
 
         if i != n_root_layers
             @. fa +=
-                flux(
+                water_flux(
                     root_depths[i],
                     model.compartment_midpoints[1],
                     ψ_soil,
@@ -623,7 +638,7 @@ function root_water_flux_per_ground_area!(
                 above_ground_area_index
         else
             @. fa +=
-                flux(
+                water_flux(
                     root_depths[i],
                     model.compartment_midpoints[1],
                     ψ_soil,
diff --git a/test/standalone/Vegetation/plant_hydraulics_test.jl b/test/standalone/Vegetation/plant_hydraulics_test.jl
index 6c3d380af5..afbc7457c1 100644
--- a/test/standalone/Vegetation/plant_hydraulics_test.jl
+++ b/test/standalone/Vegetation/plant_hydraulics_test.jl
@@ -284,13 +284,13 @@ for FT in (Float32, Float64)
             )
             # Set system to hydrostatic equilibrium state by setting fluxes to zero, and setting LHS of both ODEs to 0
             function initial_compute_exp_tendency!(F, Y)
-                AI = (; leaf = LAI(0.0), root = RAI, stem = SAI)
+                AI = (; leaf = LAI(1.0), root = RAI, stem = SAI)
                 T0A = FT(1e-8) * AI[:leaf]
                 for i in 1:(n_leaf + n_stem)
                     if i == 1
                         fa =
                             sum(
-                                flux.(
+                                water_flux.(
                                     root_depths,
                                     plant_hydraulics.compartment_midpoints[i],
                                     ψ_soil0,
@@ -307,26 +307,23 @@ for FT in (Float32, Float64)
                                     vcat(root_depths, [0.0])[2:end] -
                                     vcat(root_depths, [0.0])[1:(end - 1)]
                                 ),
-                            ) * AI[:root]
+                            ) * AI[:stem]
                     else
                         fa =
-                            flux(
+                            water_flux(
                                 plant_hydraulics.compartment_midpoints[i - 1],
                                 plant_hydraulics.compartment_midpoints[i],
                                 Y[i - 1],
                                 Y[i],
                                 PlantHydraulics.hydraulic_conductivity(
                                     conductivity_model,
-                                    ψ_soil0,
+                                    Y[i - 1],
                                 ),
                                 PlantHydraulics.hydraulic_conductivity(
                                     conductivity_model,
                                     Y[i],
                                 ),
-                            ) * (
-                                AI[plant_hydraulics.compartment_labels[1]] +
-                                AI[plant_hydraulics.compartment_labels[2]]
-                            ) / 2
+                            ) * AI[plant_hydraulics.compartment_labels[i]]
                     end
                     F[i] = fa - T0A
                 end
@@ -335,7 +332,9 @@ for FT in (Float32, Float64)
             soln = nlsolve(
                 initial_compute_exp_tendency!,
                 Vector{FT}(-0.03:0.01:0.07);
-                ftol = 1e-11,
+                ftol = eps(FT),
+                method = :newton,
+                iterations = 20,
             )
 
             S_l =
@@ -372,12 +371,8 @@ for FT in (Float32, Float64)
             canopy_exp_tendency! = make_exp_tendency(model)
             canopy_exp_tendency!(dY, Y, p, 0.0)
 
-            m = similar(dY.canopy.hydraulics.ϑ_l.:1)
-            m .= FT(0)
-            for i in 1:(n_stem + n_leaf)
-                @. m += sqrt(dY.canopy.hydraulics.ϑ_l.:($$i)^2.0)
-            end
-            @test maximum(parent(m)) < 200 * eps(FT) # starts in equilibrium
+            tolerance = sqrt(eps(FT))
+            @test all(parent(dY.canopy.hydraulics.ϑ_l) .< tolerance) # starts in equilibrium
 
 
             # repeat using the plant hydraulics model directly
@@ -394,13 +389,12 @@ for FT in (Float32, Float64)
             standalone_exp_tendency! =
                 make_compute_exp_tendency(model.hydraulics, model)
             standalone_exp_tendency!(standalone_dY, Y, p, 0.0)
-
-            m = similar(dY.canopy.hydraulics.ϑ_l.:1)
-            m .= FT(0)
-            for i in 1:(n_stem + n_leaf)
-                @. m += sqrt(dY.canopy.hydraulics.ϑ_l.:($$i)^2.0)
-            end
-            @test maximum(parent(m)) < 10^2.1 * eps(FT) # starts in equilibrium
+            @test all(
+                parent(
+                    standalone_dY.canopy.hydraulics.ϑ_l .-
+                    dY.canopy.hydraulics.ϑ_l,
+                ) .≈ FT(0),
+            )
         end
     end