Chnages to parameters and runoff, scf parameterization

experiments/integrated/fluxnet/US-Ha1/US-Ha1_simulation.jl

@@ -7,7 +7,10 @@ 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_tuple = (dz_bottom, dz_top)
+nelements = 20
+zmin = FT(-10)
+zmax = FT(0)
 # Stem and leaf compartments and their heights:
 n_stem = Int64(1)
 n_leaf = Int64(1)

experiments/integrated/fluxnet/US-MOz/US-MOz_parameters.jl

@@ -27,7 +27,7 @@ soil_vg_α = FT(0.05) # inverse meters
 ν_ss_om = FT(0.1)
 ν_ss_gravel = FT(0.0);
 z_0m_soil = FT(0.01)
-z_0b_soil = FT(0.001)
+z_0b_soil = FT(0.01)
 soil_ϵ = FT(0.98)
 soil_α_PAR = FT(0.2)
 soil_α_NIR = FT(0.2)
@@ -52,14 +52,14 @@ Drel = FT(1.6)
 g0 = FT(1e-4)
 
 #Photosynthesis model
-Vcmax25 = FT(4.5e-5) # from Yujie's paper 4.5e-5
+Vcmax25 = FT(6e-5) # from Yujie's paper 4.5e-5
 
 # Plant Hydraulics and general plant parameters
-pc = FT(-2.5e6)
+pc = FT(-2.0e6)
 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 = 3e-9 # m/s 
+K_sat_plant = 7e-8 # 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.1 * 0.0098) # Holtzman's original parameter for the bulk modulus of elasticity

experiments/integrated/fluxnet/US-MOz/US-MOz_simulation.jl

@@ -7,7 +7,10 @@ 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.1)
-
+dz_tuple = (dz_bottom, dz_top)
+nelements = 20
+zmin = FT(-10)
+zmax = FT(0)
 # Stem and leaf compartments and their heights:
 n_stem = Int64(1)
 n_leaf = Int64(1)
@@ -20,4 +23,4 @@ h_leaf = FT(9.5) # m
 t0 = Float64(0)
 
 # Time step size:
-dt = Float64(1800)
+dt = Float64(900)

experiments/integrated/fluxnet/US-NR1/US-NR1_simulation.jl

@@ -7,7 +7,10 @@ variables for running the simulation."""
 # Column dimensions - separation of layers at the top and bottom of the column:
 dz_bottom = FT(1.25)
 dz_top = FT(0.05)
-
+dz_tuple = (dz_bottom, dz_top)
+nelements = 20
+zmin = FT(-10)
+zmax = FT(0)
 # Stem and leaf compartments and their heights:
 n_stem = Int64(1)
 n_leaf = Int64(1)

experiments/integrated/fluxnet/US-Var/US-Var_parameters.jl

@@ -13,7 +13,7 @@ fluxtower site."""
 data_link = "https://caltech.box.com/shared/static/54huhy74kxnn23i6w1s54vddo2j4hl71.csv"
 
 # Height of sensor
-atmos_h = FT(2)
+atmos_h = FT(2) # from BADM
 
 # Timezone (offset from UTC in hrs)
 time_offset = 8
@@ -23,22 +23,22 @@ lat = FT(38.4133) # degree
 long = FT(-120.9508) # degree
 
 # Soil parameters
-soil_ν = FT(0.45) # m3/m3
-soil_K_sat = FT(0.45 / 3600 / 100) # m/s,
+soil_ν = FT(0.45) # m3/m3, Bonan Table 8.3
+soil_K_sat = FT(0.45 / 3600 / 100) # m/s, Bonan Table 8.3
 soil_S_s = FT(1e-3) # 1/m, guess
-soil_vg_n = FT(2.0) # unitless
-soil_vg_α = FT(2.0) # inverse meters
-θ_r = FT(0.067) # m3/m3, 
+soil_vg_n = FT(2.0) # unitless, near Bonan Table 8.3
+soil_vg_α = FT(2.0) # inverse meters, Bonan Table 8.3
+θ_r = FT(0.067) # m3/m3, near Bonan Table 8.3
 
 # Soil heat transfer parameters; not needed for hydrology only test
-ν_ss_quartz = FT(0.38)
-ν_ss_om = FT(0.0)
-ν_ss_gravel = FT(0.1);
+ν_ss_quartz = FT(0.3) # Xu and Baldocchi (2003)
+ν_ss_om = FT(0.02) # Xu and Baldocchi (2003)
+ν_ss_gravel = FT(0.0);
 z_0m_soil = FT(0.01)
-z_0b_soil = FT(0.001)
+z_0b_soil = FT(0.01)
 soil_ϵ = FT(0.98)
-soil_α_PAR = FT(0.3)
-soil_α_NIR = FT(0.4)
+soil_α_PAR = FT(0.2)
+soil_α_NIR = FT(0.2)
 
 # TwoStreamModel parameters
 Ω = FT(1.0)
@@ -57,23 +57,25 @@ Drel = FT(1.6)
 g0 = FT(1e-4)
 
 #Photosynthesis model
-Vcmax25 = FT(4.225e-5) # CLM C3 grass, Slevin et al. 2015
+Vcmax25 = FT(2 * 4.225e-5) # 2x CLM C3 grass, Slevin et al. 2015
 
 # Energy Balance model
 ac_canopy = FT(745)
 
 # Plant Hydraulics and general plant parameters
+pc = FT(-3e5)
+sc = FT(1e-3)
 SAI = FT(0)
 f_root_to_shoot = FT(1.0)
-K_sat_plant = 5e-9 # m/s # seems much too small?
+K_sat_plant = 2e-8 # m/s
 ψ63 = FT(-2.7 / 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
 conductivity_model =
     PlantHydraulics.Weibull{FT}(K_sat_plant, ψ63, Weibull_param)
 retention_model = PlantHydraulics.LinearRetentionCurve{FT}(a)
-plant_ν = FT(8.93e-4)
+plant_ν = FT(8.93e-3)
 plant_S_s = FT(1e-2 * 0.0098) # m3/m3/MPa to m3/m3/m
-rooting_depth = FT(2.6) # from Bonan Table 2.3
+rooting_depth = FT(0.3) # based off of soil depth, Xu and Baldocchi
 z0_m = FT(0.13) * h_canopy
 z0_b = FT(0.1) * z0_m

experiments/integrated/fluxnet/US-Var/US-Var_simulation.jl

@@ -4,14 +4,11 @@ variables for running the simulation."""
 
 # DOMAIN SETUP:
 
-# Column dimensions - separation of layers at the top and bottom of the column:
-dz_bottom = FT(1.0)
-dz_top = FT(0.05)
-
+# Column dimensions
 # Stem and leaf compartments and their heights:
 n_stem = Int64(0)
 n_leaf = Int64(1)
-h_leaf = FT(0.7) # m
+h_leaf = FT(0.5) # m, Xu and Baldocchi, 2003
 h_stem = FT(0) # m
 
 # TIME STEPPING:
@@ -20,4 +17,9 @@ h_stem = FT(0) # m
 t0 = Float64(0)
 
 # Time step size:
-dt = Float64(1800)
+dt = Float64(900)
+# For soil column
+nelements = 14
+zmin = FT(-0.5) #m, Xu and Baldocchi, 2003
+zmax = FT(0)
+dz_tuple = nothing

experiments/integrated/fluxnet/fluxnet_domain.jl

@@ -4,17 +4,11 @@ given in the {site-ID}_simulation.jl file in each site directory."""
 
 # Domain setup
 # For soil column
-nelements = 20
-zmin = FT(-10)
-zmax = FT(0)
 h_canopy = h_stem + h_leaf
 compartment_midpoints =
     n_stem > 0 ? [h_stem / 2, h_stem + h_leaf / 2] : [h_leaf / 2]
 compartment_surfaces = n_stem > 0 ? [zmax, h_stem, h_canopy] : [zmax, h_leaf]
 
-land_domain = Column(;
-    zlim = (zmin, zmax),
-    nelements = nelements,
-    dz_tuple = (dz_bottom, dz_top),
-)
+land_domain =
+    Column(; zlim = (zmin, zmax), nelements = nelements, dz_tuple = dz_tuple)
 canopy_domain = ClimaLand.Domains.obtain_surface_domain(land_domain)

experiments/integrated/fluxnet/plot_utils.jl

@@ -68,7 +68,13 @@ function plot_daily_avg(
         ylabel = "$var_name $(unit)",
         title = "$var_name",
     )
-    CairoMakie.lines!(ax, Array(0.5:0.5:24), data_hh_avg[:], label = label)
+    CairoMakie.lines!(
+        ax,
+        Array(0.5:0.5:24),
+        data_hh_avg[:],
+        label = label,
+        color = "blue",
+    )
     axislegend(ax, position = :lt)
     CairoMakie.save(joinpath(savedir, "$(var_name)_avg.png"), fig)
 end
@@ -107,10 +113,88 @@ function plot_avg_comp(
         ylabel = "$var_name $(units)",
         title = "$var_name: RMSD = $(round(RMSD, digits = 2)), R² = $(round(R²[1][1], digits = 2))",
     )
-    CairoMakie.lines!(ax, Array(0.5:0.5:24), model_hh_avg[:], label = "Model")
+    CairoMakie.lines!(
+        ax,
+        Array(0.5:0.5:24),
+        model_hh_avg[:],
+        label = "Model",
+        color = "blue",
+    )
 
-    CairoMakie.lines!(ax, Array(0.5:0.5:24), data_hh_avg[:], label = "Data")
+    CairoMakie.lines!(
+        ax,
+        Array(0.5:0.5:24),
+        data_hh_avg[:],
+        label = "Data",
+        color = "yellow",
+    )
     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 monthly average of a 
+variable between the model and the data. Saves the plot to the directory
+specified by savedir.
+"""
+function plot_monthly_avg_comp(
+    var_name::String,
+    ref_time,
+    model::Vector,
+    model_times::Vector,
+    data::Vector,
+    data_times::Vector,
+    units::String,
+    savedir::String,
+)
+    model_avg = compute_monthly_avg(model, model_times)
+    data_avg = compute_monthly_avg(data, data_times)
+
+    # Plot the model and data monthly cycles
+    fig = CairoMakie.Figure(size = (800, 400))
+    ax = CairoMakie.Axis(
+        fig[1, 1],
+        xlabel = "Month of Year",
+        ylabel = "$var_name $(units)",
+    )
+    CairoMakie.lines!(
+        ax,
+        Array(1:1:12),
+        model_avg[:],
+        label = "Model",
+        color = "blue",
+    )
+
+    CairoMakie.lines!(
+        ax,
+        Array(1:1:12),
+        data_avg[:],
+        label = "Data",
+        color = "yellow",
+    )
+    axislegend(ax, position = :lt)
+
+    CairoMakie.save(joinpath(savedir, "$(var_name)_monthly_avg.png"), fig)
+end
+
+"""
+    compute_monthly_avg(data::Vector, times::Vector, ref_date::Dates.DateTime)
+
+Computes the average per month of the `data` measured at `times`, where `times` is in units of seconds past the reference date `ref_date`.
+"""
+function compute_monthly_avg(
+    data::Vector,
+    times::Vector,
+    ref_date::Dates.DateTime,
+)
+    months = Dates.month.(ref_date .+ Dates.Second.(times))
+    # group by month
+    unique_months = unique(months)
+    output = zeros(length(unique_months))
+    for i in unique_months
+        mask = months .== i
+        output[i] = mean(data[mask])
+    end
+    return output, months
+end

experiments/integrated/fluxnet/run_fluxnet.jl

@@ -346,7 +346,7 @@ swc, soil_T, si = [
     ClimaLand.Diagnostics.diagnostic_as_vectors(
         d_writer,
         diag_name;
-        layer = 20, #surface layer
+        layer = nelements, #surface layer
     )[2][(N_spinup_days * 24):end] for diag_name in hourly_diag_name_2D
 ]
 dt_save = 3600.0 # hourly diagnostics

experiments/integrated/performance/conservation/ozark_conservation_setup.jl

@@ -14,7 +14,10 @@ h_leaf = FT(9.5) # m
 # TIME STEPPING:
 t0 = Float64(120 * 3600 * 24)# start mid year
 dt = Float64(150)
-
+dz_tuple = (dz_bottom, dz_top)
+nelements = 20
+zmin = FT(-10)
+zmax = FT(0)
 timestepper = CTS.ARS111()
 # Select conv. condition based on float type due to different precision
 err = (FT == Float64) ? 1e-8 : 1e-4

src/diagnostics/default_diagnostics.jl

@@ -318,6 +318,7 @@ function default_diagnostics(
             "et",
             "sr",
             "swe",
+            "snd",
             "rn",
             "lhf",
             "shf",
@@ -338,6 +339,7 @@ function default_diagnostics(
             "er",
             "sr",
             "swe",
+            "snd",
         ]
     end
 

src/diagnostics/define_diagnostics.jl

@@ -866,4 +866,15 @@ function define_diagnostics!(land_model)
             compute_snow_water_equivalent!(out, Y, p, t, land_model),
     )
 
+    # Snow depth
+    add_diagnostic_variable!(
+        short_name = "snd",
+        long_name = "Snow depth",
+        standard_name = "snow_depth",
+        units = "m",
+        comments = "The snow depth",
+        compute! = (out, Y, p, t) ->
+            compute_snow_depth!(out, Y, p, t, land_model),
+    )
+
 end

src/diagnostics/land_compute_methods.jl

@@ -432,6 +432,7 @@ end
 @diagnostic_compute "soil_internal_energy" Union{SoilCanopyModel, LandModel} Y.soil.ρe_int
 
 @diagnostic_compute "snow_water_equivalent" LandModel Y.snow.S
+@diagnostic_compute "snow_depth" LandModel p.snow.z_snow
 
 ### EnergyHydrology ###
 

src/standalone/Snow/Snow.jl

@@ -372,7 +372,7 @@ function ClimaLand.make_update_aux(model::SnowModel{FT}) where {FT}
         @. p.snow.energy_runoff =
             p.snow.water_runoff *
             volumetric_internal_energy_liq(p.snow.T, parameters)
-        @. p.snow.snow_cover_fraction = snow_cover_fraction(Y.snow.S)
+        @. p.snow.snow_cover_fraction = snow_cover_fraction(p.snow.z_snow)
     end
 end
 

src/standalone/Snow/snow_parameterizations.jl

@@ -13,15 +13,16 @@ export snow_surface_temperature,
     phase_change_flux
 
 """
-    snow_cover_fraction(x::FT; α = FT(1e-3))::FT where {FT}
+    snow_cover_fraction(x::FT; z0 = FT(1e-1), α = FT(2))::FT where {FT}
 
-Returns the snow cover fraction, assuming it is a heaviside
-function at 1e-3 meters.
-
-In the future we can play around with other forms.
+Returns the snow cover fraction from snow depth `z`, from 
+Wu, Tongwen, and Guoxiong Wu. "An empirical formula to compute 
+snow cover fraction in GCMs." Advances in Atmospheric Sciences 
+21 (2004): 529-535.
 """
-function snow_cover_fraction(x::FT; α = FT(1e-3))::FT where {FT}
-    return heaviside(x - α)
+function snow_cover_fraction(z::FT; z0 = FT(1e-1), α = FT(2))::FT where {FT}
+    z̃ = z / z0
+    return min(α * z̃ / (z̃ + 1), 1)
 end
 
 """