diff --git a/experiments/integrated/fluxnet/US-MOz/US-MOz_parameters.jl b/experiments/integrated/fluxnet/US-MOz/US-MOz_parameters.jl
index ae922ed5a9..3d696dfaf5 100644
--- a/experiments/integrated/fluxnet/US-MOz/US-MOz_parameters.jl
+++ b/experiments/integrated/fluxnet/US-MOz/US-MOz_parameters.jl
@@ -15,12 +15,12 @@ lat = FT(38.7441) # degree
 long = FT(-92.2000) # degree
 
 # Soil parameters
-soil_ν = FT(0.5) # m3/m3
-soil_K_sat = FT(4e-7) # m/s, matches Natan
-soil_S_s = FT(1e-3) # 1/m, guess
-soil_vg_n = FT(2.05) # unitless
-soil_vg_α = FT(0.04) # inverse meters
-θ_r = FT(0.067) # m3/m3, from Wang et al. 2021 https://doi.org/10.5194/gmd-14-6741-2021
+soil_ν = FT(0.55) # m3/m3
+soil_K_sat = FT(4e-7) # m/s
+soil_S_s = FT(1e-2) # 1/m, guess
+soil_vg_n = FT(2.0) # unitless
+soil_vg_α = FT(0.05) # inverse meters
+θ_r = FT(0.04) # m3/m3, 
 
 # Soil makeup
 ν_ss_quartz = FT(0.1)
@@ -44,7 +44,7 @@ G_Function = CLMGFunction(χl)
 ϵ_canopy = FT(0.97)
 
 # Energy Balance model
-ac_canopy = FT(7.5e3)
+ac_canopy = FT(5e2)
 
 # Conductance Model
 g1 = FT(141) # Wang et al: 141 sqrt(Pa) for Medlyn model; Natan used 300.
@@ -52,15 +52,17 @@ Drel = FT(1.6)
 g0 = FT(1e-4)
 
 #Photosynthesis model
-Vcmax25 = FT(9e-5) # from Yujie's paper 4.5e-5 , Natan used 9e-5
+Vcmax25 = FT(4.5e-5) # from Yujie's paper 4.5e-5
 
 # Plant Hydraulics and general plant parameters
+pc = FT(-2.5e6)
+sc = FT(5e-6)
 SAI = FT(1.0) # m2/m2 or: estimated from Wang et al, FT(0.00242) ?
 f_root_to_shoot = FT(3.5)
-K_sat_plant = 5e-9 # m/s # seems much too small?
+K_sat_plant = 3e-9 # m/s 
 ψ63 = FT(-4 / 0.0098) # / MPa to m, Holtzman's original parameter value is -4 MPa
 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
+a = FT(0.1 * 0.0098) # Holtzman's original parameter for the bulk modulus of elasticity
 conductivity_model =
     PlantHydraulics.Weibull{FT}(K_sat_plant, ψ63, Weibull_param)
 retention_model = PlantHydraulics.LinearRetentionCurve{FT}(a)
diff --git a/experiments/integrated/fluxnet/US-MOz/US-MOz_simulation.jl b/experiments/integrated/fluxnet/US-MOz/US-MOz_simulation.jl
index 727730ba85..e9c319a93d 100644
--- a/experiments/integrated/fluxnet/US-MOz/US-MOz_simulation.jl
+++ b/experiments/integrated/fluxnet/US-MOz/US-MOz_simulation.jl
@@ -6,7 +6,7 @@ variables for running the simulation."""
 
 # Column dimensions - separation of layers at the top and bottom of the column:
 dz_bottom = FT(1.5)
-dz_top = FT(0.025)
+dz_top = FT(0.1)
 
 # Stem and leaf compartments and their heights:
 n_stem = Int64(1)
@@ -20,4 +20,4 @@ h_leaf = FT(9.5) # m
 t0 = Float64(0)
 
 # Time step size:
-dt = Float64(450)
+dt = Float64(1800)
diff --git a/experiments/integrated/fluxnet/data_tools.jl b/experiments/integrated/fluxnet/data_tools.jl
index 5243362fbc..3fbc410fb3 100644
--- a/experiments/integrated/fluxnet/data_tools.jl
+++ b/experiments/integrated/fluxnet/data_tools.jl
@@ -109,9 +109,9 @@ function filter_column(driver_data::Matrix, column_name::String, units::String)
     # Check that the data exists and read it in if so
     col_dat, status = column_status(driver_data, column_name)
     # Set missing data threshold above which column is treated as absent
-    missing_threshold = 0.5
+    missing_threshold = 10.0
     # Set poor quality threshold above which the column undergoes no replacement using the quality flag
-    quality_threshold = 0.5
+    quality_threshold = 10.0
     # if it does not exist, exit
     if status == absent
         @info "Warning: Data for $column_name is absent"
diff --git a/experiments/integrated/fluxnet/plot_utils.jl b/experiments/integrated/fluxnet/plot_utils.jl
index 3bfe1cc05e..046c6e1692 100644
--- a/experiments/integrated/fluxnet/plot_utils.jl
+++ b/experiments/integrated/fluxnet/plot_utils.jl
@@ -1,7 +1,7 @@
 """Plotting utilities for the integrated fluxnet site experiments"""
 
 using Interpolations
-using Plots
+using CairoMakie
 using StatsBase
 
 S_PER_DAY = 86400 # Number of seconds in a day
@@ -61,17 +61,16 @@ function plot_daily_avg(
         compute_diurnal_avg(data, [0:data_dt:(num_days * S_PER_DAY);], num_days)
 
     # Plot the data diurnal cycle
-    plt = Plots.plot(size = (800, 400))
-    Plots.plot!(
-        plt,
-        0.5:0.5:24,
-        data_hh_avg,
-        label = label,
+    fig = CairoMakie.Figure(size = (800, 400))
+    ax = CairoMakie.Axis(
+        fig[1, 1],
         xlabel = "Hour of day",
         ylabel = "$var_name $(unit)",
         title = "$var_name",
     )
-    Plots.savefig(joinpath(savedir, "$(var_name)_avg.png"))
+    CairoMakie.lines!(ax, Array(0.5:0.5:24), data_hh_avg[:], label = label)
+    axislegend(ax, position = :lt)
+    CairoMakie.save(joinpath(savedir, "$(var_name)_avg.png"), fig)
 end
 
 """This function will be used to plot the comparison of the diurnal average of a 
@@ -101,17 +100,17 @@ function plot_avg_comp(
     R² = StatsBase.cor(model_hh_avg, data_hh_avg)^2
 
     # Plot the model and data diurnal cycles
-    plt = Plots.plot(size = (800, 400))
-    Plots.plot!(
-        plt,
-        0.5:0.5:24,
-        model_hh_avg,
-        label = "Model",
+    fig = CairoMakie.Figure(size = (800, 400))
+    ax = CairoMakie.Axis(
+        fig[1, 1],
         xlabel = "Hour of day",
         ylabel = "$var_name $(units)",
         title = "$var_name: RMSD = $(round(RMSD, digits = 2)), R² = $(round(R²[1][1], digits = 2))",
-        legend = :topleft,
     )
-    Plots.plot!(plt, 0.5:0.5:24, data_hh_avg, label = "Data")
-    Plots.savefig(joinpath(savedir, "$(var_name)_avg.png"))
+    CairoMakie.lines!(ax, Array(0.5:0.5:24), model_hh_avg[:], label = "Model")
+
+    CairoMakie.lines!(ax, Array(0.5:0.5:24), data_hh_avg[:], label = "Data")
+    axislegend(ax, position = :lt)
+
+    CairoMakie.save(joinpath(savedir, "$(var_name)_avg.png"), fig)
 end
diff --git a/experiments/integrated/fluxnet/run_fluxnet.jl b/experiments/integrated/fluxnet/run_fluxnet.jl
index 8fa0b9b554..d85f3a5700 100644
--- a/experiments/integrated/fluxnet/run_fluxnet.jl
+++ b/experiments/integrated/fluxnet/run_fluxnet.jl
@@ -4,7 +4,6 @@ import ClimaComms
 ClimaComms.@import_required_backends
 using ClimaCore
 import ClimaParams as CP
-using Plots
 using Statistics
 using Dates
 using Insolation
@@ -30,11 +29,11 @@ include(joinpath(climaland_dir, "experiments/integrated/fluxnet/data_tools.jl"))
 include(joinpath(climaland_dir, "experiments/integrated/fluxnet/plot_utils.jl"))
 
 # Read in the site to be run from the command line
-#if length(ARGS) < 1
-#    error("Must provide site ID on command line")
-#end
+if length(ARGS) < 1
+    error("Must provide site ID on command line")
+end
 
-site_ID = "US-MOz"#ARGS[1]
+site_ID = "US-MOz"
 
 # Read all site-specific domain parameters from the simulation file for the site
 include(
@@ -526,25 +525,18 @@ end
 # Water content in soil and snow
 # Soil water content
 # Current resolution has the first layer at 0.1 cm, the second at 5cm.
-plt1 = Plots.plot(size = (1500, 800))
-Plots.plot!(
-    plt1,
-    model_times ./ 3600 ./ 24,
-    swc,
-    label = "1.25cm",
-    xlim = [
-        minimum(model_times ./ 3600 ./ 24),
-        maximum(model_times ./ 3600 ./ 24),
-    ],
-    ylim = [0.05, 0.55],
-    xlabel = "Days",
-    ylabel = "SWC [m/m]",
-    color = "blue",
-    margin = 10Plots.mm,
+fig = Figure(size = (1500, 800), fontsize = 20)
+ax1 = Axis(fig[3, 1], xlabel = "Days", ylabel = "SWC [m/m]")
+limits!(
+    ax1,
+    minimum(model_times ./ 3600 ./ 24),
+    maximum(model_times ./ 3600 ./ 24),
+    0.05,
+    0.6,
 )
-
-plot!(
-    plt1,
+lines!(ax1, model_times ./ 3600 ./ 24, swc, label = "1.25cm", color = "blue")
+lines!(
+    ax1,
     model_times ./ 3600 ./ 24,
     si,
     color = "cyan",
@@ -552,74 +544,79 @@ plot!(
 )
 
 if drivers.SWC.status != absent
-    Plots.plot!(
-        plt1,
+    lines!(
+        ax1,
         data_times ./ 3600 ./ 24,
         drivers.SWC.values[data_id_post_spinup],
         label = "Data",
     )
 end
-plt2 = Plots.plot(size = (1500, 800))
-Plots.plot!(
-    plt2,
-    model_times ./ 3600 ./ 24,
-    SWE,
-    label = "Model",
-    xlim = [
-        minimum(model_times ./ 3600 ./ 24),
-        maximum(model_times ./ 3600 ./ 24),
-    ],
-    ylim = [0.0, 0.15],
-    xlabel = "Days",
-    ylabel = "SWE [m]",
-    color = "blue",
-    margin = 10Plots.mm,
+ax2 =
+    Axis(fig[2, 1], xlabel = "", ylabel = "SWE [m]", xticklabelsvisible = false)
+limits!(
+    ax2,
+    minimum(model_times ./ 3600 ./ 24),
+    maximum(model_times ./ 3600 ./ 24),
+    0.0,
+    maximum(SWE),
+)
+
+lines!(ax2, model_times ./ 3600 ./ 24, SWE, label = "Model", color = "blue")
+ax3 = Axis(
+    fig[1, 1],
+    xlabel = "",
+    ylabel = "Precipitation [mm/day]",
+    xticklabelsvisible = false,
 )
-plt3 = Plots.plot(
+limits!(
+    ax3,
+    minimum(model_times ./ 3600 ./ 24),
+    maximum(model_times ./ 3600 ./ 24),
+    -500,
+    0.0,
+)
+
+lines!(
+    ax3,
     data_times ./ 3600 ./ 24,
     (drivers.P.values .* (-1e3 * 24 * 3600) .* (1 .- snow_frac))[data_id_post_spinup],
     label = "Rain (data)",
-    ylabel = "Precipitation [mm/day]",
-    xlim = [
-        minimum(model_times ./ 3600 ./ 24),
-        maximum(model_times ./ 3600 ./ 24),
-    ],
-    margin = 10Plots.mm,
-    ylim = [-200, 0],
-    size = (1500, 400),
 )
-Plots.plot!(
-    plt3,
+lines!(
+    ax3,
     data_times ./ 3600 ./ 24,
     (drivers.P.values .* (-1e3 * 24 * 3600) .* snow_frac)[data_id_post_spinup],
     label = "Snow (data)",
-    ylabel = "Precipitation [mm/day]",
 )
-Plots.plot(plt3, plt2, plt1, layout = grid(3, 1, heights = [0.2, 0.4, 0.4]))
-Plots.savefig(joinpath(savedir, "ground_water_content.png"))
+CairoMakie.save(joinpath(savedir, "ground_water_content.png"), fig)
+
 
-plt1 = Plots.plot(size = (1500, 800))
-Plots.plot!(
-    plt1,
+
+fig2 = Figure(size = (1500, 800))
+ax12 = Axis(fig2[1, 1], xlabel = "Days", ylabel = "Temperature (K)")
+limits!(
+    ax12,
+    minimum(model_times ./ 3600 ./ 24),
+    maximum(model_times ./ 3600 ./ 24),
+    265,
+    315,
+)
+
+lines!(
+    ax12,
     model_times ./ 3600 ./ 24,
     soil_T,
     label = "1.25cm",
-    xlim = [
-        minimum(model_times ./ 3600 ./ 24),
-        maximum(model_times ./ 3600 ./ 24),
-    ],
-    xlabel = "Days",
-    ylabel = "Temperature (K)",
     color = "blue",
-    margin = 10Plots.mm,
 )
+
 if drivers.TS.status != absent
-    Plots.plot!(
-        plt1,
+    lines!(
+        ax12,
         data_times ./ 3600 ./ 24,
         drivers.TS.values[data_id_post_spinup],
         label = "Data",
     )
 end
 
-Plots.savefig(joinpath(savedir, "soil_temperature.png"))
+CairoMakie.save(joinpath(savedir, "soil_temperature.png"), fig2)
diff --git a/src/integrated/soil_canopy_root_interactions.jl b/src/integrated/soil_canopy_root_interactions.jl
index 7ff53313ec..a974455d3a 100644
--- a/src/integrated/soil_canopy_root_interactions.jl
+++ b/src/integrated/soil_canopy_root_interactions.jl
@@ -9,8 +9,8 @@ function update_root_extraction!(p, Y, t, land)
     z = land.soil.domain.fields.z
     (; conductivity_model) = land.canopy.hydraulics.parameters
     area_index = p.canopy.hydraulics.area_index
-    # allocates
-    above_ground_area_index =
+    above_ground_area_index = p.scratch1
+    above_ground_area_index .=
         PlantHydraulics.harmonic_mean.(
             getproperty(
                 area_index,