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/src/integrated/soil_canopy_model.jl b/src/integrated/soil_canopy_model.jl
index b74ff64274..347f269183 100644
--- a/src/integrated/soil_canopy_model.jl
+++ b/src/integrated/soil_canopy_model.jl
@@ -273,7 +273,7 @@ function make_update_boundary_fluxes(
 
         @. p.root_extraction =
             above_ground_area_index *
-            PlantHydraulics.flux(
+            PlantHydraulics.water_flux(
                 z,
                 land.canopy.hydraulics.compartment_midpoints[1],
                 p.soil.ψ,
diff --git a/src/standalone/Vegetation/Canopy.jl b/src/standalone/Vegetation/Canopy.jl
index 41dd15ea71..904d88cdb2 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 4502af2386..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,22 +337,23 @@ 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 a harmonic
-mean for mean K_sat between the two points (see CLM Technical Documentation). Previously,
-we used the arithmetic mean (Bonan, 2019; Zhu, 2008),
-but then when the soil K was very low, root extraction would
-continue. 
+mean for effective conducticity between the two layers
+(see CLM Technical Documentation).
 
-Following CLM, this should be modified for compartments of 
-differing sizes because the water will travel different path lengths 
-in each compartment. Hence we should weight each K as K -> K/path length.
+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::FT, z2::FT, ψ1::FT, ψ2::FT, K1::FT, K2::FT) where {FT}
+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)
@@ -624,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,
@@ -637,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 e65243bfeb..803ca74b48 100644
--- a/test/standalone/Vegetation/plant_hydraulics_test.jl
+++ b/test/standalone/Vegetation/plant_hydraulics_test.jl
@@ -215,7 +215,7 @@ for FT in (Float32, Float64)
             SAI,
             RAI,
         )
-        K_sat_plant = 1.0 # m/s.
+        K_sat_plant = 1.8e-8 # m/s.
         ψ63 = FT(-4 / 0.0098) # / MPa to m
         Weibull_param = FT(4) # unitless
         a = FT(0.05 * 0.0098) # 1/m
@@ -250,7 +250,7 @@ for FT in (Float32, Float64)
         )
 
         function leaf_transpiration(t)
-            T = FT(1) # m/s
+            T = FT(1e-8) # m/s
         end
 
         ψ_soil0 = FT(0.0)
@@ -285,12 +285,12 @@ 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(1.0), root = RAI, stem = SAI)
-                T0A = AI[:leaf]
+                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,
@@ -310,7 +310,7 @@ for FT in (Float32, Float64)
                             ) * AI[:stem]
                     else
                         fa =
-                            flux(
+                            water_flux(
                                 plant_hydraulics.compartment_midpoints[i - 1],
                                 plant_hydraulics.compartment_midpoints[i],
                                 Y[i - 1],