GitHub

Evaporation from Gilat Loess soil

This sets up the simulation that mimicks the lab experiment presented in Gardner [15] and modeled also by Lehmann and Or [16].

For further details on how to setup a simulation, please see our other Soil tutorials. This one is very terse and does not provide complete explanations

The same experiment is carried out 3 times

  1. No evaporation (zero flux boundary conditions)
  2. With evaporation but no drainage (Ksat = 0)
  3. With evaporation and drainage
using CairoMakie
import ClimaTimeSteppers as CTS
using Thermodynamics

using ClimaCore
import ClimaParams as CP
using SurfaceFluxes
using StaticArrays
using Dates
using DelimitedFiles: readdlm

using ClimaLand
using ClimaLand.Domains: Column
using ClimaLand.Soil
import ClimaLand
import ClimaLand.Parameters as LP
import ClimaLand.Simulations: LandSimulation, solve!
import SurfaceFluxes.Parameters as SFP

FT = Float64;
toml_dict = LP.create_toml_dict(FT)
earth_param_set = LP.LandParameters(toml_dict)
thermo_params = LP.thermodynamic_parameters(earth_param_set);
┌ Warning: Error requiring `IntervalArithmetic` from `TaylorSeries`
│   exception =
│    LoadError: UndefVarError: `IntervalBox` not defined in `TaylorSeries`
│    Suggestion: check for spelling errors or missing imports.
│    Stacktrace:
│      [1] top-level scope
│        @ ~/.julia/packages/TaylorSeries/ZBRjU/src/intervals.jl:171
│      [2] include(mapexpr::Function, mod::Module, _path::String)
│        @ Base ./Base.jl:307
│      [3] top-level scope
│        @ ~/.julia/packages/Requires/1eCOK/src/Requires.jl:40
│      [4] eval(m::Module, e::Any)
│        @ Core ./boot.jl:489
│      [5] EvalInto
│        @ ./boot.jl:494 [inlined]
│      [6] (::TaylorSeries.var"#__init__##4#__init__##5")()
│        @ TaylorSeries ~/.julia/packages/Requires/1eCOK/src/require.jl:101
│      [7] macro expansion
│        @ ./timing.jl:461 [inlined]
│      [8] err(f::Any, listener::Module, modname::String, file::String, line::Any)
│        @ Requires ~/.julia/packages/Requires/1eCOK/src/require.jl:47
│      [9] (::TaylorSeries.var"#__init__##2#__init__##3")()
│        @ TaylorSeries ~/.julia/packages/Requires/1eCOK/src/require.jl:100
│     [10] withpath(f::Any, path::String)
│        @ Requires ~/.julia/packages/Requires/1eCOK/src/require.jl:37
│     [11] (::TaylorSeries.var"#__init__##0#__init__##1")()
│        @ TaylorSeries ~/.julia/packages/Requires/1eCOK/src/require.jl:99
│     [12] foreach(f::typeof(invokelatest), itr::Vector{Function})
│        @ Base ./abstractarray.jl:3191
│     [13] loadpkg(pkg::Base.PkgId)
│        @ Requires ~/.julia/packages/Requires/1eCOK/src/require.jl:27
│     [14] run_package_callbacks(modkey::Base.PkgId)
│        @ Base ./loading.jl:1467
│     [15] _require_search_from_serialized(pkg::Base.PkgId, sourcepath::String, build_id::UInt128, stalecheck::Bool; reasons::Dict{String, Int64}, DEPOT_PATH::Vector{String})
│        @ Base ./loading.jl:2201
│     [16] _require_search_from_serialized
│        @ ./loading.jl:2076 [inlined]
│     [17] __require_prelocked(pkg::Base.PkgId, env::String)
│        @ Base ./loading.jl:2695
│     [18] _require_prelocked(uuidkey::Base.PkgId, env::String)
│        @ Base ./loading.jl:2560
│     [19] macro expansion
│        @ ./loading.jl:2488 [inlined]
│     [20] macro expansion
│        @ ./lock.jl:376 [inlined]
│     [21] __require(into::Module, mod::Symbol)
│        @ Base ./loading.jl:2453
│     [22] require(into::Module, mod::Symbol)
│        @ Base ./loading.jl:2429
│     [23] eval(m::Module, e::Any)
│        @ Core ./boot.jl:489
│     [24] include_string(mapexpr::typeof(identity), mod::Module, code::String, filename::String)
│        @ Base ./loading.jl:2952
│     [25] include_string
│        @ ./loading.jl:2962 [inlined]
│     [26] #61
│        @ ~/.julia/packages/Literate/On3EJ/src/Literate.jl:959 [inlined]
│     [27] task_local_storage(body::Literate.var"#61#62"{String, Bool, Module, String}, key::Symbol, val::String)
│        @ Base ./task.jl:298
│     [28] #59
│        @ ~/.julia/packages/Literate/On3EJ/src/Literate.jl:955 [inlined]
│     [29] (::IOCapture.var"#12#13"{Type{Union{}}, Literate.var"#59#60"{String, Bool, Module, String}, IOContext{Base.PipeEndpoint}, IOContext{Base.PipeEndpoint}, IOContext{Base.PipeEndpoint}, IOContext{Base.PipeEndpoint}})()
│        @ IOCapture ~/.julia/packages/IOCapture/MR051/src/IOCapture.jl:170
│     [30] with_logstate(f::IOCapture.var"#12#13"{Type{Union{}}, Literate.var"#59#60"{String, Bool, Module, String}, IOContext{Base.PipeEndpoint}, IOContext{Base.PipeEndpoint}, IOContext{Base.PipeEndpoint}, IOContext{Base.PipeEndpoint}}, logstate::Base.CoreLogging.LogState)
│        @ Base.CoreLogging ./logging/logging.jl:542
│     [31] with_logger(f::Function, logger::Base.CoreLogging.ConsoleLogger)
│        @ Base.CoreLogging ./logging/logging.jl:653
│     [32] capture(f::Literate.var"#59#60"{String, Bool, Module, String}; rethrow::Type, color::Bool, passthrough::Bool, capture_buffer::IOBuffer, io_context::Vector{Any})
│        @ IOCapture ~/.julia/packages/IOCapture/MR051/src/IOCapture.jl:167
│     [33] capture
│        @ ~/.julia/packages/IOCapture/MR051/src/IOCapture.jl:100 [inlined]
│     [34] execute_block(sb::Module, block::String; inputfile::String, fake_source::String, softscope::Bool, continue_on_error::Bool)
│        @ Literate ~/.julia/packages/Literate/On3EJ/src/Literate.jl:953
│     [35] execute_block
│        @ ~/.julia/packages/Literate/On3EJ/src/Literate.jl:935 [inlined]
│     [36] execute_markdown!(io::IOBuffer, sb::Module, block::String, outputdir::String; inputfile::String, fake_source::String, flavor::Literate.DocumenterFlavor, image_formats::Vector{Tuple{MIME, String}}, file_prefix::String, softscope::Bool, continue_on_error::Bool)
│        @ Literate ~/.julia/packages/Literate/On3EJ/src/Literate.jl:683
│     [37] (::Literate.var"#42#43"{Dict{String, Any}, IOBuffer, Module, Literate.CodeChunk, Int64})()
│        @ Literate ~/.julia/packages/Literate/On3EJ/src/Literate.jl:651
│     [38] cd(f::Literate.var"#42#43"{Dict{String, Any}, IOBuffer, Module, Literate.CodeChunk, Int64}, dir::String)
│        @ Base.Filesystem ./file.jl:112
│     [39] markdown(inputfile::String, outputdir::String; config::Dict{Any, Any}, kwargs::@Kwargs{execute::Bool, documenter::Bool, preprocess::var"#mdpre#generate_tutorial##2"{String}})
│        @ Literate ~/.julia/packages/Literate/On3EJ/src/Literate.jl:650
│     [40] markdown (repeats 2 times)
│        @ ~/.julia/packages/Literate/On3EJ/src/Literate.jl:614 [inlined]
│     [41] (::var"#generate_tutorial##0#generate_tutorial##1"{String})()
│        @ Main ~/work/ClimaLand.jl/ClimaLand.jl/docs/make.jl:47
│     [42] cd(f::var"#generate_tutorial##0#generate_tutorial##1"{String}, dir::String)
│        @ Base.Filesystem ./file.jl:112
│     [43] generate_tutorial(tutorials_dir::String, tutorial::String)
│        @ Main ~/work/ClimaLand.jl/ClimaLand.jl/docs/make.jl:37
│     [44] #9
│        @ ./none:-1 [inlined]
│     [45] (::Serialization.__deserialized_types__.var"#9#10")(t::String)
│        @ Main ~/work/ClimaLand.jl/ClimaLand.jl/docs/make.jl:0
│     [46] exec_from_cache(rr::RemoteChannel{Channel{Any}}, args::String; kwargs::@Kwargs{})
│        @ Distributed /opt/hostedtoolcache/julia/1.12.7/x64/share/julia/stdlib/v1.12/Distributed/src/workerpool.jl:387
│     [47] (::Distributed.var"#handle_msg##4#handle_msg##5"{Distributed.CallMsg{:call_fetch}})()
│        @ Distributed /opt/hostedtoolcache/julia/1.12.7/x64/share/julia/stdlib/v1.12/Distributed/src/process_messages.jl:287
│     [48] run_work_thunk(thunk::Distributed.var"#handle_msg##4#handle_msg##5"{Distributed.CallMsg{:call_fetch}}, print_error::Bool)
│        @ Distributed /opt/hostedtoolcache/julia/1.12.7/x64/share/julia/stdlib/v1.12/Distributed/src/process_messages.jl:70
│     [49] (::Distributed.var"#handle_msg##2#handle_msg##3"{Distributed.CallMsg{:call_fetch}, Distributed.MsgHeader, Sockets.TCPSocket})()
│        @ Distributed /opt/hostedtoolcache/julia/1.12.7/x64/share/julia/stdlib/v1.12/Distributed/src/process_messages.jl:287
│    in expression starting at /home/runner/.julia/packages/TaylorSeries/ZBRjU/src/intervals.jl:171
└ @ Requires ~/.julia/packages/Requires/1eCOK/src/require.jl:51

Parameters

K_sat = FT(0.01 / 3600 / 24)
vg_n = FT(1.55)
vg_α = FT(1.5)
hcm = vanGenuchten{FT}(; α = vg_α, n = vg_n)
ν = FT(0.4)
θ_r = FT(0.04)
S_s = FT(1e-3)
ν_ss_om = FT(0.0)
ν_ss_quartz = FT(0.3)
ν_ss_gravel = FT(0.0)
emissivity = FT(1.0)
z_0m = FT(1e-3)
z_0b = FT(1e-4)
d_ds = FT(0.01)# 10mm
params = ClimaLand.Soil.EnergyHydrologyParameters(
    toml_dict;
    ν,
    ν_ss_om,
    ν_ss_quartz,
    ν_ss_gravel,
    hydrology_cm = hcm,
    K_sat,
    S_s,
    θ_r,
    emissivity,
    z_0m,
    z_0b,
    d_ds,
);

start_date = DateTime(2005)
stop_date = start_date + Day(15)
dt = Float64(900.0)
SW_d = (t) -> 0
LW_d = (t) -> 294.15^4 * 5.67e-8
radiation = PrescribedRadiativeFluxes(
    FT,
    TimeVaryingInput(SW_d),
    TimeVaryingInput(LW_d),
    start_date,
)
PrescribedRadiativeFluxes{Float64}
  start_date: 2005-01-01T00:00:00
  drivers (SW_d, LW_d): AnalyticTimeVaryingInput
  cosθs: computed from date and location
  frac_diff: computed empirically

Atmos

T_air = FT(301.15)
rh = FT(0.38)
esat = Thermodynamics.saturation_vapor_pressure(
    thermo_params,
    T_air,
    Thermodynamics.Liquid(),
)
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
q_atmos = (t) -> q
h_atmos = FT(0.1)
P_atmos = (t) -> 101325
gustiness = FT(1e-2)
atmos = PrescribedAtmosphere(
    TimeVaryingInput(precip),
    TimeVaryingInput(precip),
    TimeVaryingInput(T_atmos),
    TimeVaryingInput(u_atmos),
    TimeVaryingInput(q_atmos),
    TimeVaryingInput(P_atmos),
    start_date,
    h_atmos,
    toml_dict;
    gustiness = gustiness,
)
PrescribedAtmosphere{Float64}
  start_date: 2005-01-01T00:00:00
  reference height: 0.1 m, gustiness: 0.01 m/s
  drivers: AnalyticTimeVaryingInput (liquid/snow precip, T, u, q, P)
  CO2: AnalyticTimeVaryingInput

Simulation setup - no evaporation Boundary conditions

zero_water_flux = WaterFluxBC((p, t) -> 0)
zero_heat_flux = HeatFluxBC((p, t) -> 0)
no_flux_boundary_fluxes = (;
    top = WaterHeatBC(; water = zero_water_flux, heat = zero_heat_flux),
    bottom = WaterHeatBC(; water = zero_water_flux, heat = zero_heat_flux),
);

Δz = 0.01
zmax = FT(0)
zmin = FT(-1.6)
nelems = Int((zmax - zmin) / Δz)
soil_domain = Column(; zlim = (zmin, zmax), nelements = nelems)
z = ClimaCore.Fields.coordinate_field(soil_domain.space.subsurface).z

soil = Soil.EnergyHydrology{FT}(;
    parameters = params,
    domain = soil_domain,
    boundary_conditions = no_flux_boundary_fluxes,
    sources = (),
);

Initial conditions

function estimated_ic(z)
    0.34 / (1 + exp(-(z + 0.165) / 0.005)) + 0.05
end
function set_ic!(Y, p, t0, model)
    params = model.parameters
    z = model.domain.fields.z
    FT = eltype(Y.soil.ϑ_l)
    Y.soil.ϑ_l .= estimated_ic.(z)
    Y.soil.θ_i .= 0
    T = FT(294.15)
    ρc_s = @. Soil.volumetric_heat_capacity(
        Y.soil.ϑ_l,
        Y.soil.θ_i,
        params.ρc_ds,
        params.earth_param_set,
    )
    Y.soil.ρe_int = Soil.volumetric_internal_energy.(
        Y.soil.θ_i,
        ρc_s,
        T,
        params.earth_param_set,
    )
end
timestepper = CTS.ARS111();
ode_algo = CTS.IMEXAlgorithm(
    timestepper,
    CTS.NewtonsMethod(
        max_iters = 1,
        update_j = CTS.UpdateEvery(CTS.NewNewtonIteration),
    ),
);

Saving callback

saveat = Hour(1)
saving_cb = ClimaLand.NonInterpSavingCallback(start_date, stop_date, saveat)
sv = saving_cb.affect!.saved_values
simulation = LandSimulation(
    start_date,
    stop_date,
    dt,
    soil;
    set_ic! = set_ic!,
    updateat = nothing, # disables driver updates
    solver_kwargs = (; saveat),
    timestepper = ode_algo,
    user_callbacks = (saving_cb,),
    diagnostics = (),
);
sol_no_evap = solve!(simulation);

Repeat with evaporation and drainage This requires different initial conditions

top_bc = ClimaLand.Soil.AtmosDrivenFluxBC(atmos, radiation)
evap_boundary_fluxes = (;
    top = top_bc,
    bottom = WaterHeatBC(; water = zero_water_flux, heat = zero_heat_flux),
)

soil = Soil.EnergyHydrology{FT}(;
    parameters = params,
    domain = soil_domain,
    boundary_conditions = evap_boundary_fluxes,
    sources = (),
)

timestepper = CTS.ARS111();
saveat = Hour(1)
saving_cb = ClimaLand.NonInterpSavingCallback(start_date, stop_date, saveat)
sv = saving_cb.affect!.saved_values

simulation = LandSimulation(
    start_date,
    stop_date,
    dt,
    soil;
    set_ic! = set_ic!,
    updateat = Hour(1),
    solver_kwargs = (; saveat),
    timestepper = ode_algo,
    user_callbacks = (saving_cb,),
    diagnostics = (),
);
sol = solve!(simulation);
evap = [
    parent(sv.saveval[k].soil.turbulent_fluxes.vapor_flux_liq)[1] for
    k in 1:length(sol.t)
];

# Repeat with no drainage (Ksat = 0, different BC), and with evaporation, in shorter domain
[ Info: Warning: No runoff model was provided; zero runoff generated.

This requires different boundary conditions yet again: Wet boundary at bottom, zero heat flux at bottom, the previously defined atmos driven evaporation at the top.

bottom_water_bc = MoistureStateBC((p, t) -> 0.35)
no_drainage_boundary_fluxes = (;
    top = top_bc,
    bottom = WaterHeatBC(; water = bottom_water_bc, heat = zero_heat_flux),
)
zmax = FT(0)
zmin = FT(-0.16)
nelems = Int((zmax - zmin) / Δz)
dt = Float64(10.0)
soil_domain = Column(; zlim = (zmin, zmax), nelements = nelems)
z_no_evap = ClimaCore.Fields.coordinate_field(soil_domain.space.subsurface).z
soil = Soil.EnergyHydrology{FT}(;
    parameters = params,
    domain = soil_domain,
    boundary_conditions = no_drainage_boundary_fluxes,
    sources = (),
)
saveat = Hour(1)
saving_cb = ClimaLand.NonInterpSavingCallback(start_date, stop_date, saveat);
sv = saving_cb.affect!.saved_values;
simulation = LandSimulation(
    start_date,
    stop_date,
    dt,
    soil;
    set_ic! = set_ic!,
    updateat = Hour(1),
    solver_kwargs = (; saveat),
    timestepper = ode_algo,
    user_callbacks = (saving_cb,),
    diagnostics = (),
);
sol_no_drainage = solve!(simulation);
evap_no_drainage = [
    parent(sv.saveval[k].soil.turbulent_fluxes.vapor_flux_liq)[1] for
    k in 1:length(sol.t)
];

Figures

fig = Figure(size = (800, 400))
ax = Axis(fig[1, 1], xlabel = "Day", ylabel = "Evaporation rate (mm/d)")
CairoMakie.lines!(
    ax,
    FT.(sol.t) ./ 3600 ./ 24,
    evap .* (1000 * 3600 * 24),
    label = "With drainage",
    color = :red,
)
CairoMakie.lines!(
    ax,
    FT.(sol_no_drainage.t) ./ 3600 ./ 24,
    evap_no_drainage .* (1000 * 3600 * 24),
    label = "No drainage",
    color = :blue,
)

CairoMakie.axislegend(ax)
ax2 = Axis(fig[1, 2], xlabel = "Day", ylabel = "Cumulative evaporation (mm)")
CairoMakie.lines!(
    ax2,
    FT.(sol.t) ./ 3600 ./ 24,
    cumsum(evap) .* (1000 * 3600),
    color = :red,
)
CairoMakie.lines!(
    ax2,
    FT.(sol_no_drainage.t) ./ 3600 ./ 24,
    cumsum(evap_no_drainage) .* (1000 * 3600),
    color = :blue,
)
save("evaporation_lehmann2024_figS6.png", fig);

fig2 = Figure(size = (800, 1200))
ax1 = Axis(fig2[1, 1], title = "Drainage only")
CairoMakie.ylims!(-0.35, 0)
CairoMakie.xlims!(0.0, 0.4)
linestyles = [:solid, :dash, :dashdot, :dashdotdot, :dot]
days = [0, 1, 2, 10]
for i in 1:1:4
    CairoMakie.lines!(
        ax1,
        parent(sol_no_evap.u[days[i] * 24 + 1].soil.ϑ_l)[:],
        parent(z)[:],
        label = "$(days[i]) days",
        color = :black,
        linestyle = linestyles[i],
    )
end
ax2 = Axis(fig2[2, 1], title = "Evap+Drainage", ylabel = "Depth(cm)")

CairoMakie.ylims!(-0.3, 0)
CairoMakie.xlims!(0.0, 0.4)
days = [0, 1, 2, 5, 13]
for i in 1:1:5
    CairoMakie.lines!(
        ax2,
        parent(sol.u[days[i] * 24 + 1].soil.ϑ_l)[:],
        parent(z)[:],
        label = "$(days[i]) days",
        color = :black,
        linestyle = linestyles[i],
    )
end
ax3 = Axis(fig2[3, 1], title = "Evap only", xlabel = "Volumetric Water Content")
CairoMakie.ylims!(-0.15, 0)
CairoMakie.xlims!(0.0, 0.4)
days = [0, 2, 9, 14]
for i in 1:1:4
    CairoMakie.lines!(
        ax3,
        parent(sol_no_drainage.u[days[i] * 24 + 1].soil.ϑ_l)[:],
        label = "$(days[i]) days",
        parent(z_no_evap)[:],
        color = :black,
        linestyle = linestyles[i],
    )
end

CairoMakie.axislegend(ax3, position = :lt)
CairoMakie.axislegend(ax2, position = :lt)
CairoMakie.axislegend(ax1, position = :lt)
save("evaporation_gardner_fig1.png", fig2);

This page was generated using Literate.jl.