Evaporation scheme described in the paper

NEWS.md

@@ -3,10 +3,12 @@ ClimaLand.jl Release Notes
 
 main
 -------
-- Use MODIS LAI by default in experiments and longruns
-  PR[#973](https://github.com/CliMA/ClimaLand.jl/pull/973)
+- Evaporation scheme from Lehmann et al 2008/2018
+  PR[#1023](https://github.com/CliMA/ClimaLand.jl/pull/1023)
 - Revert PR993 changes to runoff
   PR[#1021](https://github.com/CliMA/ClimaLand.jl/pull/1021)
+- Use MODIS LAI by default in experiments and longruns
+  PR[#973](https://github.com/CliMA/ClimaLand.jl/pull/973)
 
 v0.15.9
 -------

docs/tutorials/standalone/Soil/evaporation.jl

@@ -62,7 +62,7 @@ e = rh * esat
 q = FT(0.622 * e / (101325 - 0.378 * e))
 precip = (t) -> 0.0
 T_atmos = (t) -> T_air
-u_atmos = (t) -> 1.0
+u_atmos = (t) -> 0.44
 q_atmos = (t) -> q
 h_atmos = FT(0.1)
 P_atmos = (t) -> 101325
@@ -92,20 +92,20 @@ boundary_fluxes = (;
 # Define the parameters
 # n and alpha estimated by matching vG curve.
 K_sat = FT(225.1 / 3600 / 24 / 1000)
-vg_n = FT(10.0)
-vg_α = FT(6.0)
+vg_n = FT(8.91)#FT(6.34) optimized values by root finding
+vg_α = FT(3) #FT(7.339) optimized values by root finding
 hcm = vanGenuchten{FT}(; α = vg_α, n = vg_n)
 ν = FT(0.43)
-θ_r = FT(0.045)
+θ_r = FT(0.043)
 S_s = FT(1e-3)
 ν_ss_om = FT(0.0)
 ν_ss_quartz = FT(1.0)
 ν_ss_gravel = FT(0.0)
 emissivity = FT(1.0)
 PAR_albedo = FT(0.2)
 NIR_albedo = FT(0.4)
-z_0m = FT(1e-3)
-z_0b = FT(1e-4)
+z_0m = FT(1e-2)
+z_0b = FT(1e-2)
 d_ds = FT(0.01)
 params = ClimaLand.Soil.EnergyHydrologyParameters(
     FT;
@@ -126,24 +126,7 @@ params = ClimaLand.Soil.EnergyHydrologyParameters(
     d_ds,
 );
 
-# Domain - single column
-zmax = FT(0)
-zmin = FT(-0.35)
-nelems = 5
-soil_domain = Column(; zlim = (zmin, zmax), nelements = nelems)
-z = ClimaCore.Fields.coordinate_field(soil_domain.space.subsurface).z;
-
-# Soil model, and create the prognostic vector Y and cache p:
-soil = Soil.EnergyHydrology{FT}(;
-    parameters = params,
-    domain = soil_domain,
-    boundary_conditions = boundary_fluxes,
-    sources = (),
-)
-
-Y, p, cds = initialize(soil);
-
-# Set initial conditions
+# IC functions
 function hydrostatic_equilibrium(z, z_interface, params)
     (; ν, S_s, hydrology_cm) = params
     (; α, n, m) = hydrology_cm
@@ -167,12 +150,101 @@ function init_soil!(Y, z, params)
     Y.soil.ρe_int =
         Soil.volumetric_internal_energy.(FT(0), ρc_s, T, params.earth_param_set)
 end
+
+# Domain - single column
+zmax = FT(0)
+zmin = FT(-0.35)
+nelems = 28
+soil_domain = Column(; zlim = (zmin, zmax), nelements = nelems)
+z = ClimaCore.Fields.coordinate_field(soil_domain.space.subsurface).z;
+
+# Soil model, and create the prognostic vector Y and cache p:
+soil = Soil.EnergyHydrology{FT}(;
+    parameters = params,
+    domain = soil_domain,
+    boundary_conditions = boundary_fluxes,
+    sources = (),
+)
+
+Y, p, cds = initialize(soil);
+
 init_soil!(Y, z, soil.parameters);
 
 # Timestepping:
 t0 = Float64(0)
 tf = Float64(24 * 3600 * 13)
-dt = Float64(900.0)
+dt = Float64(450.0)
+
+# We also set the initial conditions of the cache here:
+set_initial_cache! = make_set_initial_cache(soil)
+set_initial_cache!(p, Y, t0);
+
+# Define the tendency functions
+exp_tendency! = make_exp_tendency(soil)
+imp_tendency! = make_imp_tendency(soil);
+jacobian! = ClimaLand.make_jacobian(soil);
+jac_kwargs =
+    (; jac_prototype = ClimaLand.FieldMatrixWithSolver(Y), Wfact = jacobian!)
+
+timestepper = CTS.ARS111();
+ode_algo = CTS.IMEXAlgorithm(
+    timestepper,
+    CTS.NewtonsMethod(
+        max_iters = 6,
+        update_j = CTS.UpdateEvery(CTS.NewNewtonIteration),
+    ),
+);
+
+# Define the problem and callbacks:
+prob = SciMLBase.ODEProblem(
+    CTS.ClimaODEFunction(
+        T_exp! = exp_tendency!,
+        T_imp! = SciMLBase.ODEFunction(imp_tendency!; jac_kwargs...),
+        dss! = ClimaLand.dss!,
+    ),
+    Y,
+    (t0, tf),
+    p,
+);
+saveat = Array(t0:3600.0:tf)
+sv_hr = (;
+    t = Array{Float64}(undef, length(saveat)),
+    saveval = Array{NamedTuple}(undef, length(saveat)),
+)
+saving_cb = ClimaLand.NonInterpSavingCallback(sv_hr, saveat)
+updateat = deepcopy(saveat)
+model_drivers = ClimaLand.get_drivers(soil)
+updatefunc = ClimaLand.make_update_drivers(model_drivers)
+driver_cb = ClimaLand.DriverUpdateCallback(updateat, updatefunc)
+cb = SciMLBase.CallbackSet(driver_cb, saving_cb);
+
+# Solve
+sol_hr =
+    SciMLBase.solve(prob, ode_algo; dt = dt, callback = cb, saveat = saveat);
+
+# Repeat at lower resolution
+zmax = FT(0)
+zmin = FT(-0.35)
+nelems = 7
+soil_domain = Column(; zlim = (zmin, zmax), nelements = nelems)
+z = ClimaCore.Fields.coordinate_field(soil_domain.space.subsurface).z;
+
+# Soil model, and create the prognostic vector Y and cache p:
+soil = Soil.EnergyHydrology{FT}(;
+    parameters = params,
+    domain = soil_domain,
+    boundary_conditions = boundary_fluxes,
+    sources = (),
+)
+
+Y, p, cds = initialize(soil);
+
+init_soil!(Y, z, soil.parameters);
+
+# Timestepping:
+t0 = Float64(0)
+tf = Float64(24 * 3600 * 13)
+dt = Float64(1800.0)
 
 # We also set the initial conditions of the cache here:
 set_initial_cache! = make_set_initial_cache(soil)
@@ -189,7 +261,7 @@ timestepper = CTS.ARS111();
 ode_algo = CTS.IMEXAlgorithm(
     timestepper,
     CTS.NewtonsMethod(
-        max_iters = 1,
+        max_iters = 6,
         update_j = CTS.UpdateEvery(CTS.NewNewtonIteration),
     ),
 );
@@ -206,80 +278,116 @@ prob = SciMLBase.ODEProblem(
     p,
 );
 saveat = Array(t0:3600.0:tf)
-sv = (;
+sv_lr = (;
     t = Array{Float64}(undef, length(saveat)),
     saveval = Array{NamedTuple}(undef, length(saveat)),
 )
-saving_cb = ClimaLand.NonInterpSavingCallback(sv, saveat)
+saving_cb = ClimaLand.NonInterpSavingCallback(sv_lr, saveat)
 updateat = deepcopy(saveat)
 model_drivers = ClimaLand.get_drivers(soil)
 updatefunc = ClimaLand.make_update_drivers(model_drivers)
 driver_cb = ClimaLand.DriverUpdateCallback(updateat, updatefunc)
 cb = SciMLBase.CallbackSet(driver_cb, saving_cb);
 
 # Solve
-sol = SciMLBase.solve(prob, ode_algo; dt = dt, callback = cb, saveat = saveat);
+sol_lr =
+    SciMLBase.solve(prob, ode_algo; dt = dt, callback = cb, saveat = saveat);
 
 # Figures
 
 # Extract the evaporation at each saved step
-evap = [
-    parent(sv.saveval[k].soil.turbulent_fluxes.vapor_flux_liq)[1] for
-    k in 1:length(sol.t)
+evap_hr = [
+    parent(sv_hr.saveval[k].soil.turbulent_fluxes.vapor_flux_liq)[1] for
+    k in 1:length(sol_hr.t)
+]
+evap_lr = [
+    parent(sv_lr.saveval[k].soil.turbulent_fluxes.vapor_flux_liq)[1] for
+    k in 1:length(sol_lr.t)
 ]
 evaporation_data = ClimaLand.Artifacts.lehmann2008_evaporation_data();
 ref_soln_E = readdlm(evaporation_data, ',')
 ref_soln_E_350mm = ref_soln_E[2:end, 1:2]
 data_dates = ref_soln_E_350mm[:, 1]
 data_e = ref_soln_E_350mm[:, 2];
 
-fig = Figure(size = (800, 400))
+fig = Figure(size = (800, 400), fontsize = 22)
 ax = Axis(
     fig[1, 1],
     xlabel = "Day",
     ylabel = "Evaporation rate (mm/d)",
-    title = "Bare soil evaporation",
+    xgridvisible = false,
+    ygridvisible = false,
 )
-CairoMakie.xlims!(minimum(data_dates), maximum(data_dates))
+CairoMakie.xlims!(minimum(data_dates), maximum(sol_lr.t ./ 3600 ./ 24))
 CairoMakie.lines!(
     ax,
     FT.(data_dates),
     FT.(data_e),
     label = "Data",
+    color = :orange,
+    linewidth = 3,
+)
+CairoMakie.lines!(
+    ax,
+    sol_lr.t ./ 3600 ./ 24,
+    evap_lr .* (1000 * 3600 * 24),
+    label = "Model, 7 elements",
     color = :blue,
+    linewidth = 3,
 )
 CairoMakie.lines!(
     ax,
-    sol.t ./ 3600 ./ 24,
-    evap .* (1000 * 3600 * 24),
-    label = "Model",
-    color = :black,
+    sol_hr.t ./ 3600 ./ 24,
+    evap_hr .* (1000 * 3600 * 24),
+    label = "Model, 28 elements",
+    color = :blue,
+    linestyle = :dash,
+    linewidth = 3,
 )
-CairoMakie.axislegend(ax)
+CairoMakie.axislegend(ax, framevisible = false)
 
 ax = Axis(
     fig[1, 2],
     xlabel = "Mass (g)",
     yticksvisible = false,
     yticklabelsvisible = false,
+    xgridvisible = false,
+    ygridvisible = false,
 )
 A_col = π * (0.027)^2
-mass_0 = sum(sol.u[1].soil.ϑ_l) * 1e6 * A_col
-mass_loss =
-    [mass_0 - sum(sol.u[k].soil.ϑ_l) * 1e6 * A_col for k in 1:length(sol.t)]
+mass_0_hr = sum(sol_hr.u[1].soil.ϑ_l) * 1e6 * A_col
+mass_loss_hr = [
+    mass_0_hr - sum(sol_hr.u[k].soil.ϑ_l) * 1e6 * A_col for
+    k in 1:length(sol_hr.t)
+]
+
+mass_0_lr = sum(sol_lr.u[1].soil.ϑ_l) * 1e6 * A_col
+mass_loss_lr = [
+    mass_0_lr - sum(sol_lr.u[k].soil.ϑ_l) * 1e6 * A_col for
+    k in 1:length(sol_lr.t)
+]
 CairoMakie.lines!(
     ax,
     cumsum(FT.(data_e)) ./ (1000 * 24) .* A_col .* 1e6,
     FT.(data_e),
-    label = "Data",
+    color = :orange,
+    linewidth = 3,
+)
+CairoMakie.lines!(
+    ax,
+    mass_loss_lr,
+    evap_lr .* (1000 * 3600 * 24),
     color = :blue,
+    linewidth = 3,
 )
 CairoMakie.lines!(
     ax,
-    mass_loss,
-    evap .* (1000 * 3600 * 24),
-    label = "Model",
-    color = :black,
+    mass_loss_hr,
+    evap_hr .* (1000 * 3600 * 24),
+    color = :blue,
+    linewidth = 3,
+    linestyle = :dash,
 )
+
 save("evaporation_lehmann2008_fig8b.png", fig);
 # ![](evaporation_lehmann2008_fig8b.png)

src/standalone/Soil/boundary_conditions.jl

@@ -696,6 +696,9 @@ These variables are updated in place in `soil_boundary_fluxes!`.
 """
 boundary_vars(bc::AtmosDrivenFluxBC, ::ClimaLand.TopBoundary) = (
     :turbulent_fluxes,
+    :q_sfc,
+    :r_sfc,
+    :ice_frac,
     :R_n,
     :top_bc,
     :top_bc_wvec,
@@ -724,6 +727,9 @@ boundary_var_domain_names(bc::AtmosDrivenFluxBC, ::ClimaLand.TopBoundary) = (
     :surface,
     :surface,
     :surface,
+    :surface,
+    :surface,
+    :surface,
     :subsurface,
     Runoff.runoff_var_domain_names(bc.runoff)...,
 )
@@ -747,6 +753,9 @@ boundary_var_types(
         Tuple{FT, FT, FT, FT, FT},
     },
     FT,
+    FT,
+    FT,
+    FT,
     NamedTuple{(:water, :heat), Tuple{FT, FT}},
     ClimaCore.Geometry.WVector{FT},
     FT,