Function Documentations

run_example()

Run example simulation located in "/examples/BEA2016-reset-FALSE-input" for 100 days and return a Dict containing the solution (solution object of DifferentialEquations.jl) and other variables useful for plotting.

Kwargs can be provided and are passed through to loadSPAC().

Example:

using LWFBrook90
example = LWFBrook90.run_example()

###
using Plots, Measures
# A) Use in-built plotting function
pl1 = plotamounts(example, :above_and_belowground, :showRWUcentroid)
pl2 = plotisotopes(example, :d18O)
pl3 = plotforcingandstates(example)
savefig(pl1, "Example_pl1.png")
savefig(pl2, "Example_pl2.png")
savefig(pl3, "Example_pl3.png")

###
# B) Construct plots yourself using the solution object
# Plot scalar solution
# Using simple plot recipe that interpolates, but without dates
plot(example.ODESolution;
    idxs = [1, 2, 3, 4, 5, 6],
    label=["GWAT (mm)" "INTS (mm)" "INTR (mm)" "SNOW (mm)" "CC (MJ/m2)" "SNOWLQ (mm)"])

# Plot vector solution
x = example.ODESolution_datetime
y = cumsum(example.ODESolution.prob.p.p_soil.p_THICK)
z = get_theta(example)
heatmap(x, y, z,
    yflip = true,
    xlabel = "Date",
    ylabel = "Depth [mm]",
    c = cgrad(:blues,  rev = false),
    colorbar_title = "θ [m3/m3]

(of fine soil volume)")#, rightmargin = 10mm)

DiscretizedSPAC

An discretization of a soil-plant-atmopsheric continuum model ready for simulation.

SPAC

An instance of a soil-plant-atmopsheric continuum model with the following fields:

• reference_date: DateTime to relate internal use of numerical days to real-world dates

• tspan: Tuple (0, Int) Time span of available input data in days since reference date

• solver_options: NamedTuple: (computeintermediatequantities, simulate_isotopes, DTIMAX, DSWMAX, DPSIMAX), containing some solver options for the model

• soil_discretization: DataFrame with contents from soil_discretizations.csv, i.e. containing columns: Upper_m,Lower_m,Rootden_,uAux_PSIM_init_kPa,u_delta18O_init_permil,u_delta2H_init_permil

    Either loaded from `soil_discretizations.csv` or then specified as argument to loadSPAC().
    Unused values are NaN.

• forcing: NamedTuple: (:meteo, :meteo_iso, :storm_durations), containing atmospheric forcings

  • forcing.meteo: DataFrame with daily atmospheric variables
  • forcing.meteo_iso: DataFrame with isotopic signatures of precipitation
  • forcing.storm_durations: DataFrame with approximate typical storm duration in hours for each month

• pars: NamedTuple: (:params, )

  • pars.params: NamedTuple: (), containing scalar parameter values for the model
  • pars.root_distribution: either:
    • NamedTuple: (beta, theta, z_rootMax_m = nothing) parametrization of root distribution with depth (f(z_m) with z_m in meters and negative downward)."
    • NamedTuple: (alpha = 0.97, z_rootMax_m = nothing) parametrization of root distribution with depth (f(z_m) with z_m in meters and negative downward)."
    • or String: "soil_discretization.csv" meaning that it must be defined in soil_discretizations.csv
  • pars.IC_scalar: NamedTuple: (), containing initial conditions of the scalar state variables
  • pars.IC_soil: initial conditions of the state variables (scalar or related to the soil). Either:
    • NamedTuple: (PSIM_init, δ18O_init, δ2H_init), containing constant values for the initial conditions
    • [NOT IMPELMENTED:]NamedTuple: (PSIM_init, δ18O_init, δ2H_init), containing functions of the initial values with argument Δz
    • or String: "soil_discretization.csv" meaning that it must be defined in soil_discretizations.csv
  • pars.canopy_evolution: canopy parameters (LAI, SAI, DENSEF, HEIGHT) as relative in percent. Either:
    • NamedTuple: (DENSEFrel = 100, HEIGHTrel = 100, SAIrel = 100, LAIrel = (DOY_Bstart, Bduration, DOY_Cstart, Cduration, LAI_perc_BtoC, LAI_perc_CtoB)), containing constant values and parameters for LAI_relative interpolation
    • or DataFrame: with daily values
  • pars.soil_horizons: DataFrame containing description of soil layers/horizons and soil hydraulic parameters

is_setup(parametrizedSPAC::SPAC)

Internal function that determines whether a SPAC has been discretized with setup().

remakeSPAC(discrSPAC::DiscretizedSPAC;
            requested_tspan = nothing,
            soil_output_depths_m::Vector = zeros(Float64, 0),
            kwargs...)

or

remakeSPAC(parametrizedSPAC::SPAC;
            requested_tspan = nothing,
            soil_output_depths_m::Vector = zeros(Float64, 0),
            kwargs...)

Generates a copy of the provided SPAC or DiscretizedSPAC and modifies all the parameter that are provided as kwargs. This is useful running the same model with a range of different parameter.

Possible kwargs are:

• soil_horizons = (ths_ = 0.4, Ksat_mmday = 3854.9, alpha_per_m = 7.11, gravel_volFrac = 0.1)
• LAI_rel = (DOY_Bstart = 115,)
• root_distribution = (beta = 0.88, z_rootMax_m = -0.6,)
• params = (DRAIN=.33, BYPAR=1, IDEPTH_m=0.67, INFEXP=0.33, ALB=0.15, ALBSN=0.7, RSSA=720., PSICR=-1.6, FXYLEM=0.4, MXKPL=16.5, MAXLAI=9.999, GLMAX=.00801, R5=235., CVPD=1.9, CINTRL=0.18,)
• IC_soil = (PSIM_init_kPa = -3.0, delta18O_init_permil = -15.55, )
• IC_scalar = ICscalar_tochange

When the argument soil_horizons contains scalar values, the parameter of the toplayer are set to this value and parameter of all other layers are modified proportionally. Alternatively, the argument soil_horizons can be provided with a vector or vectors with a value for each soil horizon. A mix of vectors and scalars is also possible.

setup(parametrizedSPAC::SPAC;
    requested_tspan = nothing,
    soil_output_depths_m::Vector = zeros(Float64, 0))

Takes the definition of SPAC model and discretize to a system of ODEs that can be solved by the package DifferentialEquations.jl

If needed, the computational grid of the soil is refined to output values at specific depths, e.g. by doing setup(SPAC; soil_output_depths_m = [-0.35, -0.42, -0.48, -1.05]).

The argument requested_tspan is a tuple defining the duration of the simulation either with a start- and an end-day relative to the reference date: setup(SPAC; requested_tspan = (0,150)). or alternatively as a tuple of 2 DataTime's. Note that the reference date given by the SPAC.reference_date refers to the day 0. Further note that it is possible to provide a tspan not starting at 0, e.g. (120,150). In that case, initial conditions are applied tspan[1], but atmospheric forcing is correctly taken into account.

simulate!(s::DiscretizedSPAC; kwargs...)

Simulates a SPAC model and stores the solution in s.ODESolution.

kwargs... are passed through to solve(SciML::ODEProblem; ...) and are documented under https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/#solver_options

DateTime2RelativeDaysFloat(x,reference_DateTime)

Transforms DateTimes x to simulation time in days since reference date.

RelativeDaysFloat2DateTime(t, reference_DateTime)

Transforms simulation time t in days since reference date to DateTimes.

loadSPAC(folder::String, prefix::String;
        # OPTIONAL ARGUMENTS ARE:
        simulate_isotopes::Bool = true,
        canopy_evolution  = "meteoveg.csv",
        Δz_thickness_m    = "soil_discretization.csv",
        root_distribution = "soil_discretization.csv",
        IC_soil           = "soil_discretization.csv",
        IC_scalar         = "initial_conditions.csv",
        storm_durations_h = "meteo_storm_durations.csv")

Create instance of SPAC model by loading different input files for LWFBrook90 from folder folder` with the names:

• [PREFIX]_meteoveg.csv
• [PREFIX]_meteo_iso.csv (needed for isotopic)
• [PREFIX]_param.csv
• [PREFIX]_meteo_storm_durations.csv
• [PREFIX]_initial_conditions.csv
• [PREFIX]_soil_horizons.csv

These files were created with an R script generate_LWFBrook90jl_Input.R that takes the same arguements as the R function LWFBrook90R::run_LWFB90() and generates the corresponding input files for LWFBrook90.jl.

Soil discretization can be provided as vector with the thickness of each cell, e.g: Δz_thickness_m = fill(0.04, 20).

Root distribution can be provided as arguments to function Rootden_() as (beta = 0.98, ) or (beta = 0.98, z_rootMax_m=-0.5) or for the gamma distribution (root_k = 1.0, root_θ_cm = 20.0, z_rootMax_m=-0.5).

Meteo storm durations can be provided as vector for each month, e.g.: [4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4]

Initial conditions in the soil can be provided as NamedTuple IC_soil = (PSIM_init_kPa = -7.0, delta18O_init_permil = -9.0, delta2H_init_permil = -11.0)

Canopy evolution can be provided as NamedTuple, giving the constant values of DENSEF, HEIGHT, SAI and dynamic evolution of LAI:

canopy_evolution = (DENSEF_rel = 100,
                    HEIGHT_rel = 100,
                    SAI_rel    = 100,
                    LAI_rel = (DOY_Bstart = 120,
                            Bduration  = 21,
                            DOY_Cstart = 270,
                            Cduration  = 60,
                            LAI_perc_BtoC = 100,
                            LAI_perc_CtoB = 20))

The vegetation season in terms of LAI is defined by defining the year into 4 parts: A->B, B->B+, B+->C, C->C+, C+->A where A is the start of the year (January 1st). effectively giving the 4 parts: C+->B, B->B+, B+->C, C->C+. The position and durationof these periods are defined in days by parameters: DOY_Bstart, Bduration, DOY_Cstart, and Cduration. LAI (in percent) is constant from C+->B and B+->C (given by LAI_perc_CtoB and LAI_perc_BtoC) and a simple linear interpolation is done for in_between (i.e. budburst and leaffall).

Initial conditions of states other than soil can be provided as NamedTuple, e.g.:

IC_scalar = (amount = (u_GWAT_init_mm = 0,
                       u_INTS_init_mm = 0,
                       u_INTR_init_mm = 0,
                       u_SNOW_init_mm = 0,
                       u_CC_init_MJ_per_m2 = 0,
                       u_SNOWLQ_init_mm =  0),
            d18O    = (u_GWAT_init_permil = -13.,
                       u_INTS_init_permil = -13.,
                       u_INTR_init_permil = -13.,
                       u_SNOW_init_permil = -13.),
            d2H     = (u_GWAT_init_permil = -95.,
                       u_INTS_init_permil = -95.,
                       u_INTR_init_permil = -95.,
                       u_SNOW_init_permil = -95.))

p_DOY(t::Float64, reference::DateTime)

Get DOY (Day Of Year) from simulation time

p_MONTHN(t::Float64, reference::DateTime)

Get Month from simulation time

read_path_param(path_param; simulate_isotopes::Bool = false)

Reads in the param.csv based on a provided path. The param.csv has a structure shown in the documentation (User Guide -> Input data). Structure of input data

function saveSPAC(sim::DiscretizedSPAC, out_dir; prefix = basename(dirname(out_dir)))

Writes input files for a simulation sim to a specified directory out_dir, with a specified file prefix.

HammelKennel_transient_root_density()

Take root density distribution parameters and generates root density in time (t) and space (z, negative downward).

Rootden_(
    [[beta],[root_k, root_θ_cm]];
    Δz_m,
    z_rootMax_m = maximum(cumsum(Δz_m)),
    z_Upper_m = 0)

Define the relative root density in each discretized soil layer either as either

• exponential distribution with mean b = -1cm/ln(β), with input argument beta (Gale and Grigal, 1987),
• gamma distribution with shape root_k>0 and scale root_θ, resulting in a distribution with mean k*θ

This function returns the instantaneous root fraction (dY/dd) (units of -/cm). The values are normalized so that the area under the curve sums up to 1.0 (when plotted vs cm). The method reduces the amount of roots in a discretization layer in case the effective maximal rooting depth comes to lie within that layer - it does not need to modify the discretization in that case.

(Normalization extends the possible range of β to values larger than 1, and thereby allows
to describe profiles that are denser at the bottom. Although it is recommended to use the gamma distribution for these cases.)

discretizesoil(; Δzm::Vector{T}, Rootden::Function, uAuxPSIMinitkPa::Function, udelta18Oinitpermil::Function = ((Δzm) -> missing), udelta2Hinitpermil::Function = ((Δzm) -> missing))

Manually generate a soil- and root-discretization for LWFBrook90.jl. This function can be used as alternative to loading an input file soil_discretization.csv.

Reqired arguments are a vector with thickness of the discretization cells/layers in meter Δ_m, and functions Rootden_(Δ_m) that generates the root density, a function uAux_PSIM_init_kPa(Δ_m) that generates the initial values of ψ (kPa) for all cells as a function of Δ_m.

Optional arguments are functions initial conditions of the two isotopic signatures based on Δ_m.

extend_lowest_horizon(soil_horizons, soil_discretizationDF)

Extends the physical description of the soil domain further downward, it this is necessary by the requested discretization.

map_soil_horizons_to_discretization(soil_horizons, soil_discretization)

Attribute soil parameters based on the defined horizons to computational cells.

refine_soil_discretization(
    Δz,
    soil_output_depths_m,
    IDEPTH_m, #=modified_SPAC.pars.params[:IDEPTH_m],
    QDEPTH_m)

Densify discretization of soil domain whenever an interface or an additional cell is needed. This is the case for interfaces

• when a new soil horizons begins or
• when a state variable needs to be extracted at a specified output depth (add an additional computational cell (consisting of upper AND lower interfaces))
• when a parameter such as IDEPTHm or QDEPTHm requires an interface (cell boundary)

get_fluxes(simulation::DiscretizedSPAC; days_to_read_out_d = nothing)

Returns a DataFrame with amounts (and isotopoic compositions) of the ecosystem water fluxes*.

By default, the values are returned for each simulation day. The user can optionally define timestep using the input argument days_to_read_out_d as:

• a numeric vector, e.g. days_to_read_out_d = 1:1.0:100 for specific days since simulation.parametrizedSPAC.reference_date

Subsets of scalar state variables or soil states can be mad afterwards with:

simout_fluxes = get_fluxes(simulation)
simout_fluxes[:, Not(r"[0-9]mm$")]    # select only scalar fluxes (by de-selecting columns containing depth information e.g. "_150mm")
simout_fluxes[:, r"(dates)|([0-9]mm)"] # select only vector fluxes (by selecting columns containing depth information e.g. "_150mm")

* fluxes returned as columns in the DataFrame are:

:dates

:cum_d_prec,:cum_d_sfal,:cum_d_sthr,:cum_d_sint,
:cum_d_rfal,:cum_d_rint,:cum_d_rthr,:cum_d_rsno,:cum_d_rnet,:cum_d_smlt,
:cum_d_irvp,:cum_d_isvp,:cum_d_snvp,:cum_d_slvp,
:cum_d_tran,
:cum_d_pint, :cum_d_ptran, :cum_d_pslvp,
:srfl,:slfl,:byfl,:dsfl,:gwfl,:vrfln,

:flow, :seep, :evap,

:TRANI_mmday_50mm, :TRANI_mmday_100mm, ... :TRANI_mmday_1050mm, :TRANI_mmday_1100mm

Returns:

366×50 DataFrame
 Row │ dates                cum_d_prec  cum_d_sfal  cum_d_sthr  cum_d_sint  cum_d_rfal  cum_d_rint  cum_d_rthr  cum_d_rsno  cum_d_rnet  cum_d_smlt   ⋯
     │ DateTime             Float64     Float64     Float64     Float64     Float64     Float64     Float64     Float64     Float64     Float64      ⋯
─────┼────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── ─
   1 │ 2021-01-01T00:00:00         0.0    0.0         0.0        0.0          0.0         0.0         0.0         0.0          0.0         0.0       ⋯
   2 │ 2021-01-02T00:00:00         0.0    0.0         0.0        0.0          0.0         0.0         0.0         0.0          0.0         0.0       ⋯

     ⋯ cum_d_irvp  cum_d_isvp  cum_d_snvp  cum_d_slvp  cum_d_tran   cum_d_pint  cum_d_ptran  cum_d_pslvp  srfl     slfl      byfl     dsfl      ⋯
     ⋯ Float64     Float64     Float64     Float64     Float64      Float64     Float64      Float64      Float64  Float64   Float64  Float64   ⋯
     ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
     ⋯   0.0        0.0         0.0        0.0         0.0            0.0       0.0            0.0            0.0   0.0          0.0      0.0   ⋯
     ⋯   0.0        0.0         0.0        0.0448292   0.0566162      2.64652   0.0566162      0.364727       0.0   0.0          0.0      0.0   ⋯

     ⋯ gwfl      vrfln     flow      seep     evap       TRANI_mmday_50mm  TRANI_mmday_100mm  TRANI_mmday_150mm  TRANI_mmday_200mm  TRANI_mmday_250mm  ⋯ TRANI_mmday_1050mm TRANI_mmday_1100mm
     ⋯ Float64   Float64   Float64   Float64  Float64    Float64           Float64            Float64            Float64            Float64            ⋯ Float64            Float64
     ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
     ⋯ 0.0       0.0       0.0           0.0  0.0             0.0                0.0                0.0                0.0                0.0          ⋯                0.0               0.0
     ⋯ 0.154122  0.154122  0.154122      0.0  0.101445        0.0434988          0.00995972         0.00240976         0.000574825        0.00013463   ⋯                0.0               0.0

get_forcing(simulation::DiscretizedSPAC)

Returns a DataFrame with amounts (and isotopoic compositions) of the input forcing. Note that meteorologic forcing (amounts and isotopic signatures) is same as in input csvs, vegetation evolution (lai, height, sai) have been transformed to absolute units.

 Row │ dates                globrad_MJDayM2  tmax_degC  tmin_degC  vappres_kPa  windspeed_ms  prec_mmDay  precdelta18O_permil  precdelta2H_permil  densef_percent  height_m  lai_m2m2  sai_m2m2
     │ DateTime             Float64          Float64    Float64    Float64      Float64       Float64     Float64              Float64             Float64         Float64   Float64   Float64
─────┼──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
   1 │ 2021-01-01T00:00:00             5.53        0.9      -10.1         0.26           1.5         0.0               -15.04             -111.96           100.0      25.0       2.1       1.0
   2 │ 2021-01-02T00:00:00             5.1        -0.6       -9.3         0.3            1.6         0.0               -15.04             -111.96           100.0      25.0       2.1       1.0

get_soil_(symbols, simulation; depths_to_read_out_mm = nothing, days_to_read_out_d = nothing))

Returns a 2D DataFrame of soil variables with soil layers as columns and time steps as rows. Supports a number of variables: - :θ (= :theta) = volumetric soil moisture values (m3/m3) - :ψ (= :psi) = soil matric potential (kPa) - :W = soil wetness (-) - :SWATI = soil water volumes contained in discretized layers (mm) - :K = soil hydraulic conductivities (mm/day) - :δ18O, :δ2H, :d18O, :d2H = isotopic signatures (delta) - TRANI, RWU = root water uptake flux from each cell (mm/day) The user can define timesteps as days_to_read_out_d or specific depths as depths_to_read_out_mm, that are both optionally provided as numeric vectors, e.g. depths_to_read_out_mm = [100, 150] or days_to_read_out_d = 1:1.0:100

Function to read out soil variables from a simulated simulation.

• symbols: can be a single symbol (:θ) or a vector of symbols [:θ] or [:θ, :ψ, :δ18O, :δ2H, :W, :SWATI, :K]
  • Please note that also non-Unicode symbols are accepted: e.g. [:theta, :psi, :delta18O, :delta2H, :d18O, :d2H]
• simulation: is a DiscretizedSPAC that has been simulated.
• depths_to_read_out_mm: either nothing or vector of Integers.
• days_to_read_out_d: either nothing or vector of Floats representing days.

Examples getsoil(:θ, simulation) getsoil([:θ, :ψ, :K], simulation; depthstoreadoutmm = [100, 200, 500, 1200])

get_states(simulation::DiscretizedSPAC; days_to_read_out_d = nothing)

Returns a DataFrame with amounts (and isotopoic compositions) of the inputs and state variables: PREC, GWAT, INTS, INTR, SNOW, (amount only: CC, SNOWLQ), (signature only: XYLEM). RWU flux can be obtained with get_fluxes(). By default, the values are returned for each simulation day. The user can optionally define timestep using the input argument days_to_read_out_d as:

• :integrator_step to have it each simulation time step
• a numeric vector, e.g. days_to_read_out_d = 1:1.0:100 for specific days since simulation.parametrizedSPAC.reference_date

Subsets of scalar state variables or soil states can be mad afterwards with:

simout_states = get_states(simulation)
simout_states[:, Not(r"[0-9]mm$")]    # select only scalar states (by de-selecting columns containing depth information e.g. "_150mm")
simout_states[:, r"(dates)|([0-9]mm)"] # select only vector states (by selecting columns containing depth information e.g. "_150mm")

Returns:

366×108 DataFrame
 Row │ dates                GWAT_mm  INTS_mm   INTR_mm       SNOW_mm  CC_MJm2     SNOWLQ_mm  SWAT_mm  GWAT_d18O  INTS_d18O  INTR_d18O  SNOW_d18O  XYLEM_d18O  ⋯ δ18O_permil_50mm  δ18O_permil_100mm  δ18O_permil_1 ⋯ ψ_kPa_1100mm
     │ DateTime             Float64  Float64   Float64       Float64  Float64     Float64    Float64  Float64    Float64    Float64    Float64    Float64     ⋯ Float64           Float64            Float64       ⋯ Float 64
─────┼──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
   1 │ 2021-01-01T00:00:00      1.0  0.0        0.0           0.0     0.0          0.0       90.2641  -13.0        -13.0      -13.0      -13.0     -10.1111  ⋯         -10.1111           -10.1111           -10. ⋯ -7
   2 │ 2021-01-02T00:00:00      1.0  0.0        0.0           0.0     0.0          0.0       90.4532  -12.5973     -15.04     -15.04     -15.04    -10.1111            -10.1111           -10.1111           -10. ⋯ -6.43

get_water_partitioning(simulation)

Returns three 2D DataFrame of water fluxes with different fluxes as columns and time steps as rows. The three DataFrames are in daily, monthly and yearly resolution and span the entire simulation.

Examples dfpartdaily, dfpartmonthly, dfpartyearly = getwaterpartitioning(simulation)

plotamounts(simulation::DiscretizedSPAC)
plotamounts(simulation::DiscretizedSPAC, compartments::Symbol)
plotamounts(simulation::DiscretizedSPAC, compartments::Symbol, RWUcentroid::Symbol)

Plots the amount results of a SPAC Simulation. By default both above and belowground. The user can override this with the second argument isotope as one of :aboveground, :belowground, or :above_and_belowground. RWUcentroid can have values of either :dontShowRWUcentroid or :showRWUcentroid.

plotamounts(simulation::DiscretizedSPAC)
plotamounts(simulation::DiscretizedSPAC, compartments::Symbol)
plotamounts(simulation::DiscretizedSPAC, compartments::Symbol, RWUcentroid::Symbol)

Plots the amount results of a SPAC Simulation. By default both above and belowground. The user can override this with the second argument isotope as one of :aboveground, :belowground, or :above_and_belowground. RWUcentroid can have values of either :dontShowRWUcentroid or :showRWUcentroid.

plotamounts(simulation::DiscretizedSPAC)
plotamounts(simulation::DiscretizedSPAC, compartments::Symbol)
plotamounts(simulation::DiscretizedSPAC, compartments::Symbol, RWUcentroid::Symbol)

Plots the amount results of a SPAC Simulation. By default both above and belowground. The user can override this with the second argument isotope as one of :aboveground, :belowground, or :above_and_belowground. RWUcentroid can have values of either :dontShowRWUcentroid or :showRWUcentroid.

plotforcingandstates(simulation::DiscretizedSPAC)

Plots the forcing, states and major fluxes as results of a SPAC Simulation.

plotforcingandstates(simulation::DiscretizedSPAC)

Plots the forcing, states and major fluxes as results of a SPAC Simulation.

plotforcingandstates(simulation::DiscretizedSPAC)

Plots the forcing, states and major fluxes as results of a SPAC Simulation.

plotisotopes(simulation::DiscretizedSPAC)
plotisotopes(simulation::DiscretizedSPAC, isotope::Symbol)
plotisotopes(simulation::DiscretizedSPAC, isotope::Symbol, (d18O = (-6, -16), d2H = (-125, -40)))
plotisotopes(simulation::DiscretizedSPAC, isotope::Symbol, (d18O = (-6, -16), d2H = (-125, -40)), RWUcentroid::Symbol))

Plots the isotope results of a SPAC Simulation. By default both δ18O and δ2H. The user can override this with the second argument isotope as one of :d18O, :d2H, or :d18O_and_d2H. RWUcentroid can have values of either :dontShowRWUcentroid or :showRWUcentroid.

plotisotopes(simulation::DiscretizedSPAC)
plotisotopes(simulation::DiscretizedSPAC, isotope::Symbol)
plotisotopes(simulation::DiscretizedSPAC, isotope::Symbol, (d18O = (-6, -16), d2H = (-125, -40)))
plotisotopes(simulation::DiscretizedSPAC, isotope::Symbol, (d18O = (-6, -16), d2H = (-125, -40)), RWUcentroid::Symbol))

Plots the isotope results of a SPAC Simulation. By default both δ18O and δ2H. The user can override this with the second argument isotope as one of :d18O, :d2H, or :d18O_and_d2H. RWUcentroid can have values of either :dontShowRWUcentroid or :showRWUcentroid.

plotisotopes(simulation::DiscretizedSPAC)
plotisotopes(simulation::DiscretizedSPAC, isotope::Symbol)
plotisotopes(simulation::DiscretizedSPAC, isotope::Symbol, (d18O = (-6, -16), d2H = (-125, -40)))
plotisotopes(simulation::DiscretizedSPAC, isotope::Symbol, (d18O = (-6, -16), d2H = (-125, -40)), RWUcentroid::Symbol))

Plots the isotope results of a SPAC Simulation. By default both δ18O and δ2H. The user can override this with the second argument isotope as one of :d18O, :d2H, or :d18O_and_d2H. RWUcentroid can have values of either :dontShowRWUcentroid or :showRWUcentroid.

prepare_for_LWFBrook90R(s::DiscretizedSPAC; return_value = "inputs", out_dir = ".", file_suffix = "")

Transforms a SPAC model to input for the R package LWFBrook90R. Its return value depends on the argument returnvalue ∈ ["results", "inputs", "csvpaths"] The options are:

• "results": run LWFBrook90R and return the simulation results in Julia without writing to disk
• "inputs": return the input as DataFrames in Julia without writing to disk
• "csv_paths": write the input as .csv files to disk and return the paths ()

the current working directory for running with R.

Note that if return_value = "results" a working installation of R and LWFBrook90R v0.5.3 is required.

define_LWFB90_u0()

Generate vector u0 needed for ODE() problem in DiffEq.jl package.

HammelKennel_lateral_rootgrowth(;final_relden, tini_yrs, INITRDEP_m, INITRLEN_m_per_m2, RGROPER_yrs, age_yrs)

Compute root growth according to LWF Bayern root growth model, (Hammel and Kennel 2000).

Arguments:

- `final_relden[i]`  :  final relative values of root length per unit volume in a specific layer i
- `INITRDEP_m`       :  intial root depth, m
- `INITRLEN_m_per_m2`:  initial water-absorbing root length per unit area, m m-2
- `t_y`              :  current age of vegetation, yrs
- `tstart_y[i]`      :  initial age for root growth in a specific layer i, yrs
- `RGROPER_yrs`      :  period of root growth in layer, yrs

Returns: - RELDEN[] : current, age-dependent relative values of root length per unit volume

define_diff_eq_parameters()

Generate vector p needed for ODE() problem in DiffEq.jl package.

Arguments

• parametrizedSPAC::SPAC: Instance of a definition of a SPAC-model (soil-plant-atmosphere continuum)
• compute_intermediate_quantities::...: TODO argument description.
• simulate_isotopes::...: TODO argument description.
• soil_output_depths_m: vector of depths at which state variables should be extractable (negative numeric values [in meter])

generate_canopy_timeseries_relative(canopy_evolution; days, reference_date)

Take the canopy_evolution argument of a SPAC and generate the realtive time series.

interpolate_aboveground_veg(veg_evolution::DataFrame)

Take canopy evolution defined as DF and generate parameters continuous in time.

interpolate_meteo(
    meteo_forcing::DataFrame,
    meteo_iso_forcing::Union{DataFrame,Nothing})

Take meteorologic parameters in input_meteoveg and input_meteoiso and generate continuous parameters.

make_from_relative_to_absolute(;
    canopy_evolution::DataFrame,
    p_MAXLAI,
    p_SAI_baseline_,
    p_DENSEF_baseline_,
    p_AGE_baseline_yrs,
    p_HEIGHT_baseline_m)

Take canopy evolution defined as DataFrame in relative terms and generate a DataFrame containing the time series in absolute terms.

MSBDAYNIGHT()

Calculate average daily rate of potential and actual interception, evaporation, and transpiration by considering weighted average of rate during day and rate during night: Subroutine MSBDAYNIGHT - day-night loop: Compute day and night rates

Arguments

• many

http://www.ecoshift.net/brook/b90.html: " Subroutine MSBDAYNIGHT contains routines that calculate the five components of evaporation (see Flow Chart):

• evaporation of intercepted rain (IRVP)
• evaporation of intercepted snow (ISVP)
• evaporation from snow (SNVP)
• soil evaporation (SLVP) from the top soil layer
• transpiration (TRANI) from each soil layer that contains roots

These evaporation values are obtained separately for daytime and nightime, then combined into a daily values. Note that evaporation of intercepted storage and snow evaporation (IRVP, ISVP, SNVP) are reduced if their sources disappear. (This happens later in MSBPREINT.)

Potential evaporation rates are obtained using the Shuttleworth and Wallace (1985) modification of the Penman-Monteith approach. Daily solar radiation is corrected for slope, is allocated to the daytime, and is converted to average daytime rate. Subroutine AVAILEN calculates available energy (net radiation minus SHEAT=0) at the top (AA) and at the bottom (ASUBS) of the canopy, using a Beers Law extinction coefficient. The three aerodynamic resistances (RAA, RAC, RAS) needed by the Shuttleworth-Wallace method are obtained in subroutine SWGRA, using algorithms of Shuttleworth and Gurney (1990). These resistances depend on leaf area index (LAI), which can vary seasonally, and on canopy height, which determines stem area index (SAI). The canopy surface resistance to transpiration (RSC) for the daytime is obtained in subroutine SRSC; it depends on maximum leaf conductance, reduced for humidity, temperature, and light penetration. At night RSC is the reciprocal of leaf area index (LAI) times minimum leaf conductance (GLMIN). Soil evaporation resistance (RSS) depends on soil water potential in the top soil layer. Subroutine SWPE uses AA, ASUBS, RSC, RSS, RAA, RAC, RAS and the vapor pressure deficit (VPD) to calculate potential transpiration (PTR) and the associated ground or soil evaporation (GER) as given by Shuttleworth and Wallace (1985). Subroutine SWPE is then called again with RSC = 0 to give the intercepted evaporation rate and its associated soil evaporation (PIR and GIR). Subroutine TBYLAYER obtains actual transpiration by layer (TRANI). Actual transpiration is the lesser of potential transpiration and a soil water supply rate determined by the resistance to liquid water flow in the plants and on root distribution and soil water potential in the soil layers. If the actual transpiration is less than the potential, a new, higher GER is calculated by subroutine SWGE. After the MSBDAYNIGHT day-night loop, these evaporation rates are weighted for daytime and nighttime according to daylength (DAYLEN), and the daily average rates are then used in later calculations.

The "precipitation interval" is equal to one day when daily precipitation is provided along with other daily weather data in the Data File; parameter NPINT is then set to 1 and the precipitation loop is passed through once a day. " http://www.ecoshift.net/brook/pet.html " BROOK90 obtains evaporation rates separately for daytime and nighttime within a day-night evaporation loop. All solar radiation (SOLRAD) is assigned to the daytime. The atmospheric humidity (EA) is assumed constant through the day ("day" refers to 24 hours). The daytime and nighttime values of air temperature and wind speed are obtained in subroutine WEATHER using function WNDADJ. Vapor pressure deficit (VPD) is obtained using subroutine ESAT. Subroutine CANOPY uses function INTERP to obtain canopy structure variables for the day. Subroutine ROUGH gets canopy roughness parameters. Within a day-night loop, the three aerodynamic resistances needed by the Shuttleworth-Wallace method are calculated in subroutine SWGRA. The canopy surface resistance (RSC) for the daytime is obtained from subroutine SRSC, and the soil surface resistance (RSS) in function FRSS. Subroutine SWPE uses function PM along with the various resistances to obtain potential transpiration rate (PTR) and the associated ground or soil evaporation rate (GER) by the Shuttleworth-Wallace equations. Subroutine SWPE is called again with RSC = 0 to give potential interception rate (PIR) and its associated soil evaporation rate (GIR). Subroutine TBYLAYER obtains actual transpiration by layer (ATRANI). If the actual transpiration is less than the potential, a new, higher GER is calculated by subroutine SWGE. BROOK90 then weights the daytime and nighttime rates by the solar daylength (DAYLEN) to obtain average rates for the day, PTRAN, GEVP, PINT, GIVP, and TRANI, which are used in later calculations. "

MSBDAYNIGHT_postprocess()

Calculate average daily rate of potential and actual interception, evaporation, and transpiration by considering weighted average of rate during day and rate during night: Subroutine MSBDAYNIGHT_postprocess - Combine day and night rates to average daily rate

Arguments

• many

http://www.ecoshift.net/brook/pet.html

MSBITERATE()

http://www.ecoshift.net/brook/b90.html

Subsurface water movement is determined in several to many iterations per precipitation time-step. Remaining model calculations are done within subroutine MSBITERATE for each iteration loop.

Net throughfall (RNET) plus snowmelt (SMLT) may:

1. infiltrate into the soil matrix of the surface horizon (INFLI(1)),
2. infiltrate directly to deeper horizons via vertical macropore flow (INFLI),
3. go immediately to streamflow via vertical macropore flow followed by downslope pipe flow (BYFLI),
4. or go immediately to streamflow via impaction on a variable saturated source area (SRFL).

The fraction of area acting as a saturated source area (SAFRAC) is obtained in subroutine SRFLFR. Source area flow (SRFL) is obtained as SAFRAC plus impervious area (IMPERV) times RNET + SMLT. Infiltration rate (SLFL) is RNET + SMLT - SRFL. The fraction of infiltration to each layer that bypasses the layer and becomes an output via bypass flow (BYFLI) is calculated in subroutine BYFLFR. For each layer, the downslope flow rate by matrix flow (DSFLI) is obtained from subroutine DSLOP. In general, one or more of SRFL, BYFL, and DSFL will be set to zero by the user.

If the water potential difference between layers is less than the parameter DPSIMX, vertical flow (VRFLI) is zero; otherwise subroutine VERT obtains VRFLI between layers from a weighted hydraulic conductivity and the water potential difference between the layers. For the bottom layer, outflow to groundwater is the hydraulic conductivity of the layer times a parameter (DRAIN), which can vary from 0 to 1. This assumes a gravity potential gradient.

Subroutine INFLOW is called to get net inflow into each layer (NTFLI) using parameter DTIMAX as a first approximation for iteration time step. The rate of change of matric potential with water content (DPSIDW) from function FDPSIDW is used with NTFLI in subroutine ITER to obtain the maximum iteration time step (DTI) allowed by two parameters. The parameters are DPSIMX and the maximum allowed change in soil water content (DSWMAX). INFLOW is called again with the new DTI to get the final NTFLI, VRFLI, BYFLI, and matric uptake (INFLI).

Groundwater discharge to streamflow (GWFL) and deep seepage (SEEP) are obtained from subroutine GWATER. GWFL is simulated as a fixed fraction of groundwater each day and SEEP is a fixed fraction of GWFL.

Simulated streamflow is the sum of SRFL, BYFL, DSFL, and GWFL. This can be compared with measured streamflow if that is available.

At the end of each iteration time-step, soil water content of each layer (SWATI) is updated by adding NTFLI * DTI. Groundwater storage is also updated. Then new soil water variables are calculated using function FPSIM and subroutine SOILVAR. Available water (AWAT) is calculated for output as water held in the root zone between field capacity and PSICR.

SWCHEK tests that the SWATI remain between zero and saturation; if not, the program stops. At the end of each day the water balance is checked and the program stops if there is an error greater than 0.003 mm. These crashes should only occur if there are parameter or programming errors.

MSBPREINT()

http://www.ecoshift.net/brook/b90.html In subroutine MSBPREINT, precipitation is separated into rain and snow using SNOFRC. If NPINT = 1, subroutine INTER24 is called twice, once for snow interception and once for rain interception; this routine uses the monthly parameter DURATN, which is the average storm duration in hours. If NPINT > 1, subroutine INTER is used instead, and the precipitation is assumed to occur over the whole precipitation interval. Transpiration (TRAN) for the interval is summed over layers and reduced by the fraction of time the canopy is wet; soil evaporation (SLVP) is GIR when the canopy is wet and GER when it is dry. If a snowpack exists, subroutine SNOWPACK is called to use SNOEN and the rain and snow throughfall to calculate snowmelt (SMLT), cold content (CC), and liquid water content of the snow. Net rain to the soil surface (RNET) is rain throughfall (RTHR) minus rain absorbed by the snowpack (RSNO). Water reaching the ground surface is RNET + SMLT.

MSBSETVARS()

Compute state dependent parameters for updating states INTS, INTR, SNOW, CC, SNOWLQ in callback function.

A) get: 1) sunshine durations, 2) SFAL, 3) plant resistance (TODO: could be done before simulation) B) getwindspeedfromcanopyandsnowpack: (pfuUADTM, pfuUANTM) = f(...) C) getpotentialsnowpackandsoilevaporationfromcanopy,snowpack,atmosphere,: (pfuUADTM, pfuUANTM) = f(...)

Arguments

• many

" From ecoshift: Subroutine MSUBSETVARS contains subroutines that calulate derived variables for the day. SUNDS, CANOPY, ROUGH, and PLNTRES are called to get solar, canopy structure, roughness, and plant resistance variables that depend on day of the year. Subsurface heat flux (SHEAT) is always set to zero. Subroutine WEATHER estimates missing values, modifies input data as necessary, corrects weather station wind speed to wind speed above the canopy, and determines daytime and nighttime temperatures and wind speeds. Subroutine SNOFRAC determines the fraction of daily precipitation that is snow (SNOFRC). If there is no snow on the ground, the soil evaporation resistance (RSS) is obtained as function FRSS. When there is snow on the ground, the snowpack temperature (TSNOW) is calculated from cold content (CC). Subroutine SNOVAP then estimates the snow evaporation rate. Subroutine SNOENRGY obtains the energy available for snowmelt (SNOEN) from mean daily temperature. The factor is modified for canopy cover as determined by LAI and SAI. Snow evaporation or condensation depends on the aerodynamic resistances and the vapor gradient; however, an arbitrary reduction factor is required. "

AbstractSoilHydraulicParams

Represents an abstract parametrization of soil hydraulics. Examples of soil hydraulic parametrizations are:

• Mualem-van Genuchten
• Clapp-Hornberger

A summary is available in Shao, Y. and Irannejad, P.: On the Choice of Soil Hydraulic Models in Land-Surface Schemes, Boundary Layer Meterol., 90, 83–115, https://doi.org/10.1023/A:1001786023282, 1999.

ClappHornbergerSHP

Clapp-Hornberger parametrization of soil hydraulics

\[\frac{\psi}{\psi_s} = w^{-1/\lambda} \qquad \textup{if: }\psi \leq \psi_i \\ \psi = -m(w-n)(w-1) \qquad \textup{if: }\psi_i \leq \psi \leq 0 \\ K(\theta) = K_s \left( \frac{\theta}{\theta_s} \right)^{(3+2/\lambda)} \\ \\ \textup{Parameters:} \\ m = \frac{\psi_i}{(1-w_i)^2} - \frac{\psi_i}{w_i(1-w_i)\lambda} \\ n = 2w_i - \frac{\psi_i}{mw_i\lambda} -1\]

source: Shao, Y. and Irannejad, P.: On the Choice of Soil Hydraulic Models in Land-Surface Schemes, Boundary Layer Meterol., 90, 83–115, https://doi.org/10.1023/A:1001786023282, 1999.

Represents a discretized 1D column of soil with Clapp-Hornberger parametrization.

Input fields: pTHICK, pSTONEF, pTHSAT, pPSIF, pTHETAF, pKF, pBEXP, pWETINF Derived fields: NLAYER, pCHM, pCHN, pTHETAF, pPSIG, pSWATMAX, pWETF

Represents a discretized 1D column of soil with Mualem-van Genuchten parametrization.

Input fields: pTHICK, pSTONEF, pTHSAT, pKθfc, pKSAT, pMvGα, pMvGn, pMvGl, pθr Derived fields: pPSIF, pTHETAF, pPSIG, pSWATMAX, pWETF

MualemVanGenuchtenSHP

Mualem-van Genuchten parametrization of soil hydraulics

\[\Theta = \frac{\theta - \theta_r}{\theta_s - \theta_r} = \left( \frac{1}{1+(-\alpha \psi)^n} \right)^m \\ n = \frac{1}{m-1} \\ K(\theta) = K_s \Theta^{1/2}\left[ 1 - (1 - \Theta^{1/m})^m \right ]^2 \\ K(\psi) = K_s\frac{\left[ 1- (-\alpha\psi)^{n-1} (1 + (-\alpha \psi)^n)^{-m} \right ]^2}{\left[ 1 + (-\alpha\psi)^n \right ]^{m/2}}\]

source: Shao, Y. and Irannejad, P.: On the Choice of Soil Hydraulic Models in Land-Surface Schemes, Boundary Layer Meterol., 90, 83–115, https://doi.org/10.1023/A:1001786023282, 1999.

FDPSIDWF!(dψδW, u_aux_WETNES, p_soil)

Compute derivative dψi/dWi for each layer i.

Ecoshift: FDPSIDW returns dψi/dWi for one layer, which is needed for the selection of iteration time-step.

For Clapp-Hornberger: Differentiation of (7) and (4) leads to dψi / dWi = ( -b ψf / Wf ) ( Wi / Wf )-b-1 in the unsaturated range, dψi / dWi= m ( 2 Wi - n - 1 ) in the near saturation range, and dψi / dWi = 0 when the soil is saturated (Wi = 1).

For Mualem-van Genuchten: dψi / dWi = TODO.... write documentation here

FK_MvG(WETNES, KSAT, MvGl, MvGn)

Compute hydraulic conductivity from wetness for the Mualem van Genuchten parametrization.

Compute unsaturated hydraulic conductivity: K(Se) a.k.a. K(W) using MvG equation 8) K = KsW^l[ 1 - (1-W^(1/m))^m ]^2 using m = 1-1/n yields: K = KsW^l[ 1 - (1-W^(n/(n-1)))^(1-1/n) ]^2

FPSIM(u_aux_WETNES, p_soil)

Compute ψ(Se) = h(Se) a.k.a ψ(W) = h(W).

FTheta(u_aux_WETNES, p_soil)

Compute θ based on Se.

FWETNES(u_aux_PSIM, p_soil)

Computes θ(ψ) = θ(h) by computing first Se(ψ)=Se(h) a.k.a W(ψ)=W(h)

SWCHEK(u_SWATI)

Correct u_SWATI and throw error if too far away.

derive_auxiliary_SOILVAR(u_SWATI, p_soil)

Derive alternative representations of soil water status.

Based on the state u_SWATI it returns (u_aux_WETNES, u_aux_PSIM, u_aux_PSITI, p_fu_KK)

• u_aux_WETNES: wetness, fraction of saturation
• u_aux_PSIM: matric soil water potential for layer, kPa
• p_PSIG: gravity potential, kPa
• u_aux_PSITI: total potential ψt = ψm + ψg (sum of matrix potential and gravity potential)
• u_aux_θ: volumetric soil water content, m3/m3
• p_fu_KK: unsaturated hydraulic conductivity: K(Se) a.k.a. K(W)

BYFLFR!(p_fu_BYFRAC, p_QFPAR, p_QFFC, u_aux_WETNES, p_soil)

Compute fraction of bypass flow.

Ecoshift

" Bypass flow (BYFL) and surface or source area flow (SRFL) are the two stormflow or quickflow generating mechanisms in BROOK90. The conceptual difference is that SRFL is "new" water that has not infiltrated but has moved across the surface to a channel, whereas BYFL is "new" water that has moved to a channel below the surface via macropores or pipes. In BROOK90 the amount of BYFL from each layer depends on the wetness of that particular layer and on the amount of infiltration to it, which is controlled by INFEXP. The amount of SRFL depends on the total wetness of all soil layers down to and including input parameter QDEPTHm . In general users should not try to model BYFL and SRFL simultaneously because trying to fit parameters for both at the same time would be too complicated. The parameter BYPAR is set to 1 to allow BYFL and zero to prevent it. SRFL is prevented by setting both QDEPTHm and IMPERV to zero. The same parameters, QFFC and QFPAR, are used for SRFL and for BYFL from all layers.

When BYPAR = 0, there is no bypass flow from deeper layers, but bypass flow is still generated from layer 1 when the layer would otherwise become oversaturated. When INFEXP = 0 or IDEPTH_m = 0, BYFL can only be generated from layer 1 because it is the only layer receiving infiltrated water.

In each iteration loop, subroutine BYFLFR calculates the fraction of the water infiltrating to each layer that becomes bypass flow (BYFRACi) as

BYFRACi = QFFC ^ [ 1 - (1 / QFPAR) * (WETNESi - WETFi) / (1 - WETFi ) ]

where WETNESi is the layer wetness, WETFi is the layer wetness at field capacity, and QFFC and QFPAR are parameters. QFFC is the bypass fraction at field capacity (WETNESi = WETFi). QFPAR represents the fraction of the water content between field capacity and saturation at which BYFRAC reaches 1 (Fig. WAT-2). BYFRAC increases exponentially with layer wetness and is prevented from exceeding 1. Raising QFFC will raise BYFL proportionally at all water contents. Raising QFPAR will increase BYFL for soil drier than field capacity and decrease it for soil above field capacity (Fig. WAT-2).

With BYPAR = 1 and QFPAR = 0, BYFRAC is 0 below field capacity and 1 above it, and QFFC is ignored. With a single soil layer and no SRFL, this produces a classic "bucket" model, with the leakiness of the bucket determined by DRAIN. With multiple layers and BYPAR = 1, the bucket model does not work because a layer at field capacity diverts all excess water to BYFL, preventing wetting of deeper layers. See SRFLFR for a modified bucket model with multiple layers.

When QFPAR = 1, BYFRAC reaches 1 when the layer is saturated (Fig WAT-2). A very large QFPAR produces a constant BYFRAC of QFFC.

Note that BYFRAC is calculated from soil water prior to the input of water for the time step. "

DSLOP()

Compute downslope flow rate from layer.

GWATER(u_GWAT, p_GSC, p_GSP, aux_du_VRFLIN)

Calculate groundwater flow and seepage loss.

INFLOW(NLAYER, DTI, p_INFRAC, p_fu_BYFRAC, p_fu_SLFL,
aux_du_DSFLI, aux_du_TRANI, aux_du_SLVP, p_SWATMAX, u_SWATI, VRFLI_prior)

Compute net inflow to soil layer.

Ecoshift

" In this routine, infiltrating water (SLFL) is allocated to soil water in each layer (INFLIi ) and to bypass flow from each layer (BYFLIi ). The fraction of SLFL going to each layer (INFRACi ) is constant and is obtained in subroutine INFPAR. This fraction is separated into water to bypass flow (BYFLIi ) and water to the soil matrix (INFLIi ) by the bypass flow fraction (BYFRACi ) from subroutine BYFLFR. The routine then calculates net inflow to each layer, including withdrawal by transpiration and soil evaporation. INFLOW is called once each iteration time step and then is called once again if subroutine ITER produces a new, shorter, iteration time step.

The INFLOW routine is passed through once for each layer, working from the bottom up, preventing oversaturation of any layer. The total water allocated to layer i is

INFIL = INFRACi * SLFL.

Bypass flow rate is

BYFLIi = BYFRACi * INFIL

and the infiltration to soil matrix water in the layer is

INFLIi = INFIL - BYFLIi.

INFLOW next determines the maximum inflow rate (MAXINi ), to the layer that can be allowed in the iteration time-step. The vertical flow to the next layer down (VRFLIi ) (which may be negative), the transpiration withdrawal (TRANIi ), and the downslope flow (DSFLIi ) are fixed. So

MAXINi = (SWATMAXi - SWATIi) / DTI + VRFLIi + DSFLIi + TRANIi

where DTI is the iteration time step.

If VRFLIi-1 + INFLIi exceeds MAXIN, then oversaturation will occur. If BYFRACi > 0, INFLIi is first reduced toward zero, then, if necessary, VRFLIi-1 is reduced, or even made negative if VRFLIi is negative. BYFLIi is increased by the amount that INFLIi is reduced. If BYFRACi = 0, VRFLIi-1 is reduced or even made negative. INFLOW finally calculates the net water flux rate, NTFLIi into each soil layer.

NTFLI(I%) = VV(I% - 1) + INFLI(I%) - VV(I%) - DSFLI(I%) - TRANI(I%)

where VV is the final value of VRFLI.

In the top layer, soil evaporation withdrawal is also added to MAXIN. Because there is no VRFLI(0) to reduce, excess water becomes negative INFLI(1) and increases BYFLI(1).

The modified values of VRFLIi are output from the INFLOW routine as variable VV because the original VRFLIi are needed again if the iteration time step (DTI) is reduced. "

INFPAR(p_INFEXP, p_ILAYER, p_soil)

Compute fraction of infiltration to each soil layer.

Arguments:

• p_INFEXP: infiltration exponent: 0 all to top, 1 uniform with depth, >1.0=more at bottom than at top
• p_ILAYER: number of layers over which infiltration is distributed
• p_soil

KKMEAN(KK_i, KK_iplus1, THICK_i, THICK_iplus1)

Compute average hydraulic conductivity.

Note that between version 3.1 (where LWFBrook90 was forked), and 4.8 of Brook90 there were different variants how to compute the average.

VERT()

Compute vertical flow rate.

AVAILEN(SLRAD, p_fu_ALBEDO, p_C1, p_C2, p_C3, p_fT_TA, p_fT_VAPPRES, RATIO, p_fu_SHEAT, p_CR,
p_fu_LAI, p_fu_SAI)

Estimate the available energy above and below the canopy.

Ecoshift

" Estimates of available energy above and below the canopy are made in subroutine AVAILEN. Available energy is net radiation minus subsurface heat flux (SHEAT), and is the energy available for partitioning into heating the air and evaporating water. SHEAT is set to zero in code in MSBSETVARS.

AVAILEN calculates the available energies on an instantaneous basis, that is, in W/m2. MSBDAYNIGHT supplies to AVAILEN the average daytime solar radiation on the given slope (SLRAD, W/m2) as

SLRAD = SLFDAY * SOLRAD / (WTOMJ * DAYLEN)

where SOLRAD is the input daily solar radiation (MJ m-2) on a horizontal surface, and WTOMJ converts W m-2 to MJ m-2 d-1.

Net solar radiation (SOLNET) is

SOLNET = (1 - ALBEDO) * SLRAD

where ALBEDO is the albedo of the surface (above any canopy). When there is snow on the ground, albedo is the parameter ALBSN, otherwise it is the parameter ALB.

Estimation of net longwave radiation has been the subject of much research. BROOK90 uses Brutsaert's (1982) equation for effective clear sky emissivity (EFFEM)

EFFEM = 1.24 * [ EA * 10 / ( TA + 273.15 ] 1 / 7

where EA is the vapor pressure in kPa, and TA is the air temperature in °C. Alternative formulations for clear sky emissivity generally differ by 30 to 50 W/m2 in net longwave radiation under clear sky for the same temperature and humidity. This is equivalent to 1.1 to 1.8 mm/d of evaporated water, a substantial amount! The differences among equations were not systematic but depended on temperature and humidity (see the next section). Consequently the differences are reduced when totalled over a long time, and are further reduced by cloudiness. The Brutsaert equation is the most central.

A cloud cover correction (CLDCOR) to net longwave radiation has been used widely in the form

CLDCOR = C3 + ( 1 - C3 ) * NOVERN

where NOVERN is the sunshine duration for the day, normally written as n/N, or the fraction of possible hours of sunshine. C3 is a fixed parameter in BROOK90. NOVERN is obtained by inverting a common relation of solar radiation to sunshine duration to the form

NOVERN = ( RATIO - C1 ) / C2 , but not < 0 or > 1,

where RATIO is SOLRAD / I0HDAY, or the ratio of solar radiation to potential radiation on a horizontal surface, and C1 and C2 are empirical parameters.

The net longwave radiation (LNGNET) is

LNGNET = ( EFFEM - 1 ) * CLDCOR * SIGMA * ( TA + 273.15 )4

where SIGMA is the Stefan-Boltzmann constant. Then net radiation (RN) is

RN = SOLNET + LNGNET

and available energy above the canopy (AA) is

AA = RN - SHEAT.

Net radiation is assumed to be reduced exponentially down through the canopy according to an extinction coefficient (CR) and the sum of leaf area index (LAI) and stem area index (SAI), so the available energy at the ground (ASUBS) is

ASUBS = RN * exp [-CR * ( LAI + SAI ) ] - SHEAT.

Applying the extinction coefficient for photosynthetically-active radiation (CR) to net radiation as well is theoretically incorrect. The extinction coefficient for net radiation should be less than for PAR. But given the other uncertainties in both uses of CR there is not much point in differentiating the two values.

Clear-sky Longwave Radiation

A number of equations exist for estimating downward longwave radiation from clear sky, EFFEM, based only on weather station temperature and humidity. I have compared these methods by showing the net longwave radiation (to avoid the fourth power temperature response of the downward flux) as a function of temperature for three different relative humidities. I have not seen such a presentation in the literature, but it seems to show the differences among methods quite well.

Fig SUN-1Fig. SUN-1. Net longwave radiation from a clear sky estimated by five different equations, as a function of temperature for three relative humidities. The Y axis should read "Net Longwave Radiation from Clear Sky".

The methods differ only in how they express EFFEM (e):

Brunt (1932) e = ca + cb ea0.5 Satterlund (1979) e = 1.08 {1 - exp[-(10 ea)K / 2016]} Brutsaert (1982) e = 1.24 (10 ea / K)1/7 Idso and Jackson (1969) e = 1 - 0.261 exp(-0.000777 T2) Swinbank (1963) e = 0.0000092 K2 where ea is the vapor pressure in kPa, T is the temperature in °C, and K is the temperature in °K. The (mostly) upper Brunt curve has ca = 0.65, cb = 0.134 (kPa-0.5) (Fitzpatrick and Stern 1966), the middle curve has ca = 0.52, cb = 0.206 (Brunt 1932), and the lower curve has ca = 0.44, cb = 0.253 (Penman 1948). Note that the Swinbank and Idso-Jackson curves have no humidity dependence. Brutsaert (1982) admits that Satterlund's equation matches the data better below freezing, but Satterlund has a very flat temperature response and is the highest method of all for temperatures of 0 to 20°C and humidities above 60%. The Brutsaert method tends to be the most central over the whole range of possible conditions.

At most temperatures, the range of net longwave radiation is about 50 W/m2, which is equivalent to an evaporation of 1.8 mm/d! The choice of method for clear-sky emissivity thus plays a major role in the value of PE estimates when net radiation is estimated. A major effort using worldwide longwave data (not estimates from models!) will be needed to improve this situation. Fortunately, the cloud cover correction, approximate as it is, brings the net longwave closer to zero and helps wash out the emissivity error. "

EQUIVSLP(p_LAT, SLOPE, p_ASPECT)

Correct solar radiation for slope and aspect and compute "equivalent slope" parameters.

Ecoshift

" Correction of solar radiation for slope and aspect requires calculation of an "equivalent slope" in subroutine EQUIVSLP. The equivalent slope is defined as the location on the earth's surface where a horizontal surface is parallel to the given sloping surface. Following Swift (1976), L1 is the latitude of this "equivalent slope" (which is actually a horizontal surface), and L2 is the difference in hour angle (longitude) between the two locations. For any given slope and aspect, L1 and L2 need be found only once; they do not change over time. So they are calculated in EQUIVSLP at the beginning of B90, as:

L1 = ASIN[ COS(SLOPE) * SIN(LAT) + SIN(SLOPE) * COS(LAT)* COS(ASPECT) ] and L2 = ATAN { SIN(SLOPE) * SIN(ASPECT) /[ COS(SLOPE) * COS(LAT) - SIN(SLOPE) * SIN(LAT) * COS(ASPECT) ] }

with fixes for negative or zero denominator in L2, where SLOPE (ESLOPE), ASPECT, and latitude (LAT) are input parameters describing the location. All angles in the subroutine are in radians. "

SUNDS()

Return pfTDAYLEN, pfTI0HDAY, pfTSLFDAY.

Ecoshift

" Several radiation-related variables depend only on day of the year and location. These are calculated in SUNDS, which is called once a day.

SUNDS requires the solar constant (SCD), which is the radiation (W/m2) on a surface normal to the sun outside the atmosphere. It depends on day-of-the-year (DOY) to determine the earth-sun distance and is

SCD = SC / (1 - .0167 * COS(.0172 * (DOY - 3))) ^ 2

where SC is the solar constant at the mean earth-sun distance (Swift 1976). SC is set to 1367 W/m2 (Lean 1991) and can not be changed.

The declination of the sun (DEC) is the angle by which the sun is above or below the plane of the earth's equator. DEC is zero at the equinoxes and +23.5° or -23.5 at the solstices. Swift (1976) gives the solar declination (radians) as:

DEC = ASIN { 0.39785 * SIN [ 4.86961 + 0.017203 * DOY + 0.033446 * SIN ( 6.224111 + 0.017202

• DOY ) ] }

The daily integral of potential insolation on a slope (I0SDAY, MJ/m2) is given by Swift (1976) as:

I0SDAY = WTOMJ * SCD * FUNC3(DEC, L2, L1, T3, T2)

where

FUNC3 = (π/ 2) { SIN(DEC) * SIN(L1) * (T3 - T2) + COS(DEC) * COS(L1) * [ SIN(T3 + L2) - SIN(T2 + L2) ] }

is a program function and T2 and T3 are the hour angles of sunrise and sunset on the slope, which are obtained from function HAFDAY using the latitude of the equivalent slope. FUNC3 has units of d-1 and WTOMJ is the conversion factor 0.0864 (MJ m-2 d-1) / (W m-2). The actual SUNDS algorithm is more complicated than this because it must consider the possibility of two sunrises and sunsets on the slope in one day. The details of the algorithm are given by Swift (1976). Note that this algorithm assumes that the "opposing slope" is horizontal. In reality in mountainous terrain, the potential insolation is further reduced by any distant terrain that obscures the horizon more than the given slope itself does. The calculation of such obscuration is difficult and is outside the scope of BROOK90.

The daily integral of potential insolation on a horizontal surface (I0HDAY, MJ/m2) is found from the I0SDAY equation with L1 = LAT, L2 = 0, and T3 and T2 for a horizontal surface at LAT. The daylength (DAYLEN), which is the fraction of a day that the sun is above a horizontal horizon, is HAFDAY / π where function HAFDAY is used with L = LAT. SLFDAY is the ratio of I0SDAY to I0HDAY. SUNDS outputs DAYLEN, I0HDAY, and SLFDAY. "

" CANOPY_snowCover computes plant "parameters" that can change due to the evolution of a snowpack: The height of the canopy above any snowpack, h (HEIGHT), is

HEIGHT = RELHIT * MAXHT - pfuSNODEP

where MAXHT is the maximum height for the year, which is an input parameter, and RELHIT is the relative height for the day of the year (doy), as obtained with function INTERP from the RELHT parameter array. HEIGHT is not allowed to be less than 0.01 m, which gives an appropriate roughness parameter for "smooth" surfaces. The snowpack depth (pfuSNODEP, m) is the snow water content (SNOW, mm) divided by 1000 times snow density (SNODEN), which is assumed constant. Although snow density can actually vary from 0.05 to 0.5, the constant value is good enough to account for burying of the canopy in BROOK90. The RATIO of uncovered HEIGHT to total height (RELHT * MAXHT) is also calculated. "

CANOPY_timeEvolution()

Compute evolution of plant parameters over the season.

Ecoshift

" Subroutine CANOPY_timeEvolution calculates plant "parameters" that can vary with day of the year ( DOY).

Actual projected leaf area index, Lp (LAI), is

LAI = MAXLAI * RELLAI(DOY) * DENSEF * RATIO

where MAXLAI is the maximum LAI for the year, RELLAI(DOY) is the relative LAI for the doy as obtained with function INTERP from the RELLAI parameter array, DENSEF is a thinning parameter between zero and one (see below). The use of RATIO assumes that LAI is distributed uniformly with height. LAI is prevented from being less than 0.00001 to avoid zero divides; this can cause small amounts of transpiration, which may be ignored.

Actual projected stem area index Sp (SAI), is assumed proportional to HEIGHT following Federer et al. (1996), so

SAI = CS * HEIGHT * DENSEF

where CS is a parameter that is the ratio of SAI to HEIGHT.

Total root length per unit area (RTLEN) is

RTLEN = MXRTLN * RELHT * DENSEF

where MXRTLN is the maximum root length for the year. Correction for seasonal RELHT assumes that root length increases proportionally with height growth.

The total plant resistance to water movement (RPLANT) is

RPLANT = 1 / (MXKPL * RELHT * DENSEF)

where MXKPL is the plant conductivity at maximum height growth. RPLANT is not allowed to be greater than 1E8 MPa d/mm, which is effectively infinite. Correction for seasonal RELHT assumes that canopy conductance increases proportionally with height growth.

DENSEF is normally 1.0 in the above four equations. This parameter was included in the model as a convenient way to "thin" a canopy by removing a fraction of the plants. LAI, SAI, and RTLEN are all reduced proportionally to DENSEF, and RPLANT is increased. However DENSEF does NOT reduce HEIGHT because the remaining canopy still has the same height. Therefore DENSEF should NOT be set to 0 to simulate a clearcut as HEIGHT is unchanged and the aerodynamic resistances will be wrong. Probably DENSEF should not be less than 0.05. "

ESAT(TA_degC)

Calculate saturated vp (kPa) and DELTA=dES/dTA (kPa/K) from temperature based on Murray J Applied Meteorol 6:203 (Magnus-Tetens) using as input TA_degC (air temperature in °C).

PM(AA, VPD, DELTA, RA, RC)

Compute Penman-Monteith latent heat flux density, W/m2.

Ecoshift

" The Penman-Monteith equation is

\[L_v ρ_w E = \frac{Δ(Rn - S) + c_p ρ D_a / r_a}{Δ + γ + γ(r_c/r_a)}\]

or alternatively:

\[L_v ρ_w E = \frac{r_a \cdot Δ \cdot(Rn - S) + c_p ρ D_a}{r_a \cdot (Δ + γ) + γ \cdot r_c }\]

where E is the evaporation rate in volume of water per unit land area per unit time, Lv is the latent heat of vaporization for water, ρw is the density of water, Δ is the rate of change of vapor pressure with temperature, Rn is the net radiation above the surface, S is the subsurface heat flux, cp is the heat capacity of air, ρ is the density of air, Da is the vapor pressure deficit in the air, γ is the psychrometer constant, rc is the "canopy resistance", and ra is the aerodynamic resistance between the canopy and a reference height za at which Da is measured. The vapor pressure deficit, Da, is ea* - ea. The equation assumes that the vapor pressure at the effective evaporating surface, e0, is the saturated vapor pressure at the surface temperature. Then rc and ra are the two "resistances" through which water vapor passes as it moves down the vapor pressure gradient from e0 to ea. The canopy resistance, rc, represents resistance to flow of vapor through the stomates and cuticle of individual leaves and through the air around each leaf to some "effective" source height of water vapor in the plant canopy. The aerodynamic resistance, ra, is a measure of the turbulent transfer capability of the atmosphere between the effective source height and za. The Penman-Monteith equation is derived from the energy balance equation and the mass transfer equations for sensible and latent heat fluxes (e.g. Brutsaert 1982). "

ROUGH()

Compute canopy roughness height.

Ecoshift:

ROUGH obtains the roughness parameter, z0 , and the zero-plane displacement, d, based on canopy height, h, the projected leaf area index, Lp, and the projected stem area index, Sp. The methods used follow Shuttleworth and Gurney (1990) with some modifications. Shuttleworth and Gurney (1990) defined plant canopies as either "closed" or "sparse" based on whether Lp is greater or less than some arbitrary value Lpc, which they take as 4. Following Federer et al. (1996), BROOK90 defines a closed canopy as having Lp + Sp greater than Lpc + Spc, where Spc is taken as cs h, as described in the previous section. Spc is not reduced by DENSEF. RATIO is ( Lp + Sp ) / ( Lpc + Spc ) (this RATIO differs from RATIO in subroutine CANOPY). When RATIO is greater than or equal to 1, the canopy is "closed" and z0 and d are the values for a closed canopy.

The roughness parameter for closed canopies, z0c (Z0C), has been determined experimentally from wind profile data to be about 13% of h for relatively short canopies, but only 5% or less of h for some forests (Brutsaert 1982; Shuttleworth 1989). Following Federer et al. (1996) BROOK90 parameterizes the ratio z0c / h as follows: czs is the ratio for smooth surfaces with h less than hs; czr is the ratio for rough surfaces with h greater than hr; and for h between hs and hr, z0c is interpolated linearly between z0 = czs hs at h = hs and z0 = czr hr at h = hr. The effect of stems and branches is assumed to be included in z0c.

The zero-plane displacement of a closed canopy, dc (DISPC), is generally related to canopy height, h, and to z0c by

(22) dc = h - z0c / 0.3

(Monteith 1973; Shuttleworth and Gurney 1990).

When h (HEIGHT) is small, it is possible for z0c to be less than z0g (Z0GS), which is the parameter giving roughness of the ground surface under the canopy. This impossible situation would cause problems later on. It is prevented by reducing z0g to the value of z0c, which works nicely because z0c is always added to dc.

The values of d and z0 for sparse canopies (RATIO < 1) are obtained by a modification of the method that Choudhury and Monteith (1988) developed from curves of Shaw and Pereira (1982). A drag coefficient per unit leaf area and stem area, cd (CDRAG), is calculated from z0c and dc as

\[(23)\qquad c_d = 4 \frac{(-1 + \exp(0.909 - 3.03 \frac{z_{0c}}{h} )}{(L_{pc} + S_{pc})} .\]

This assumes that a unit of Spc produces the same drag as a unit of Lpc. Then

\[(24)\qquad d = 1.1 h \ln(1 + (c_d(L_p+S_p))^{0.25}).\]

Shuttleworth and Gurney (1990) also used this approach, but assumed a constant cd of 0.07, which is not appropriate when z0c / h is not 0.13.

For z0, Shuttleworth and Gurney (1990) chose between a sparse canopy equation and the closed canopy value of 0.3 (h - d), depending on a fixed value of cd Lp. However, that procedure leads to a step change at the transition. The step change is avoided by choosing the minimum of the two z0 values, as

\[(25)\qquad z_0 = \min(0.3(h-3), z_{0g} + 0.3 h (c_d(L_p+S_p))^{0.5})\]

where z0g is the roughness parameter of the ground surface. When a snowpack is present, the value of z0g is replaced by the parameter Z0S before ROUGH is entered.

The reference height at which weather variables are known, za (ZA), is set to a fixed amount (ZMINH) above h. It thus is assumed to move up and down if h changes through the year.

SRSC()

Compute canopy surface resistance.

Ecoshift

" This routine obtains the canopy surface resistance, rsc, which is the classic canopy resistance in the Penman-Monteith equation, using the Jarvis (1976) expression for the factors that control the individual leaf resistance, r, and its reciprocal the leaf conductance, g.

Confusion arises because r and g may be given on the basis of both total leaf surface area and projected leaf area. To straighten this out, consider that evaporation is the product of leaf conductance, g, in units of m3H2O m-2leaf s-1, and leaf area index, L, in units of m-2leaf m-2ground to give a canopy conductance in m3H2O m-2ground s-1. But the leaf area index can be given for projected leaf area, Lp, or for total leaf area, Lt. Similiarly, the conductance can be for projected leaf area, gp, or for total leaf area, gt. But the total evaporation must be the same, so gp Lp = gt Lt is required. In the literature, gt Lt is commonly used for needle-leaved plants and gp Lp for broadleaved plants. Following Körner et al. (1979) BROOK90 uses gp Lp for all plants. For needles, where gt, or its reciprocal rt, is known, then gp = ρ tp gt = ρtp / rt, where ρtp = Lt / Lp. Note that ρtp can vary from 2 for flat needles to π for cylindrical needles. For broad leaves, the leaf resistance is usually given separately as rb for the abaxial (lower) surface, and rd for the adaxial (upper) surface. Each of these is based on the one-sided area of the leaf or the projected leaf area index, so gp Lp = gb Lp + gd Lp = Lp ( 1 / rb + 1 / rd ). For amphistomatous leaves (stomates on both sides) the simplification gb = gd is often assumed, giving gp Lp = 2 gb. The factor of 2 can easily be confused with ρtp = 2, but they are not the same. For hypostomatous leaves (stomates only on the abaxial side), gd = 0 is often assumed, giving gp Lp = gd Lp. BROOK90 uses only gp, which hereafter is called gl, and is independent of assumptions about how gp was obtained.

Stomates open and close in response to several external and internal variables. Jarvis (1976) proposed that the effects could be considered as multiplicative such that

\[(26)\qquad g_l = g_{lmin} + (g_{lmax} - g_{lmin}) f_T f_D f_R f_W f_C\]

where gl is the leaf conductance, glmin is its minimum value (closed stomates), glmax is its maximum value, and fT, fD, fR, fW, and fC are reduction factors, varying between 0 and 1, that account for effects of temperature, vapor pressure deficit, radiation (light), leaf water stress, and atmospheric carbon dioxide concentration, respectively, on stomatal opening. This multiplicative expression has been widely used, though more by default and for simplicity than because of any extensive empirical testing. BROOK90 is not concerned with effects of changing atmospheric CO2 effects, so fC = 1 always.

BROOK90 also always uses fw = 1. The resulting rsc is therefore appropriate to the definition of potential transpiration as occurring when soil water is at field capacity. Reduction of fw below 1 in response to drying soil would allow a direct calculation of actual transpiration using the Penman-Monteith or Shuttleworth-Wallace equations. However, BROOK90 does not do this; it simulates actual transpiration as the lesser of potential transpiration using fw = 1 and a supply rate controlled by the water potential gradient and plant resistance.

Of the remaining three dependencies, the temperature response is the least well-documented. Jarvis (1976) proposed a skewed parabolic response, such that stomatal conductance is reduced both by low and high temperatures. BROOK90 uses a slightly simpler response consisting of an optimum temperature range with fT = 1 and inverted half parabolas to reduce fT to 0 at each end of the range, so

\[(27)\\ \begin{aligned} f_T &= 0 &\text{if } &T_a <T_L \\ f_T &= 1-(\frac{T_1-T_a}{T_1-T_L})^2 &\text{if }T_L < &T_a < T_1\\ f_T &= 1 &\text{if }T_1 < &T_a < T_2 \\ f_T &= 1-(\frac{T_a-T_2}{T_H-T_2})^2 &\text{if }T_2 < &T_a < T_H\\ f_T &= 0 &\text{if }T_H < &T_a \\ \end{aligned}\]

T1, T2, and TH are all parameters. If T1 = T2, there is no optimum range. If TL = T1 and T2 = TH, fT is a square wave function. The use of mean daily temperature, Ta, here instead of some other temperature like minimum daily temperature is arbitrary. The actual temperature response is probably much more complex than this and involves acclimation.

Stomates are known to partially close in response to greater dryness of the air surrounding the leaf, though the mechanism is unclear. The dryness is expressed best by the vapor pressure difference between the leaf and its surrounding air, but as this is not generally known, the atmospheric vapor pressure deficit above the canopy, Da (VPD), is usually used. The error induced by using Da is certainly not larger than the uncertainty in the magnitude of the vapor deficit response for most species. Jarvis (1976) assumed a linear conductivity response, but Lohammar et al. (1980) suggest a linear resistance response such that

\[(28)\qquad f_D = c_D / ( c_D + D_a )\]

where cD (CVPD) is a constant, which is the Da at which fD = 0.5. This has the advantage of approaching 0 as an asymptote when Da is large.

For fR also, several functional forms have been used. BROOK90 uses the form given by Stewart (1988) in which fR = 1 when solar radiation incident on the leaf, RL, is equal to its nominal maximum value, Rm (1000 W/m2), so

\[(29)\qquad f_R = \frac{(R_m+R_0)R_L}{R_m(R+R_0)}\]

and R0 is a second parameter. A more convenient parameter than R0 is R.5, defined as the radiation level at which fR = 0.5. Then from (29)

\[(30)\qquad R_0 = \frac{R_m}{\frac{R_m}{R_{.5}} - 2}\]

For most plants, the value of R.5 is relatively low, around 50 to 100 W/m2. In these expressions, the solar radiation at the leaf is used, whereas the stomatal opening is actually influenced by the photosynthetically active radiation. Use of solar radiation, RL, therefore assumes that the spectral distribution of the radiation does not change with depth into the canopy, which is incorrect (Federer and Tanner 1966, can't resist the opportunity to quote my ancient Ph.D. work here!).

The canopy conductance, gsc, which is the reciprocal of rsc, is the integral of gl over each increment of total leaf area in the canopy , dL' , so

\[(31)\qquad g_{sc} = \frac{1}{r_{sc}} = \int_0^{L_p} g_l dL'\]

With (26), where fT, fD, and glmax are all assumed constant through the canopy, (31) becomes

\[(32)\qquad \frac{1}{r_c} = g_c = L_p g_{lmin} + f_T f_D(g_{lmax} - g_{lmin}) \int_0^{L_p} fR dL'\]

where fR is the only variable that depends on location in the canopy.

Integration of fR requires an expression for the penetration of solar radiation into the canopy. A Beer's Law or exponential extinction is commonly used, such that the average radiation flux density on the leaf surface, RL, at any level in the canopy depends on the projected leaf area, L', and projected stem area, S', above that level, as follows:

\[(33)\qquad R_L = C_R R exp(-C_R (L' + S'))\]

where CR is the extinction coefficient and R is the solar radiation at the top of the canopy. The first CR accounts for leaf inclination assuming random azimuthal distribution (Monteith 1973, Campbell 1977, Shuttleworth and Gurney 1990). The S' term was added by Federer et al. (1996) and assumes that a unit of projected leaf area and a unit of projected stem area have the same absorption effect. Assuming that only half of Sp is distributed proportionally to L' (and the other half is below the leaves in a "stem space"), then L' + S' = L' + (L'/Lp)(Sp/2) = fs L' where fs = (Lp + Sp/2) / Lp. The integral in (32) using (29) and (33) is then

\[(34)\qquad \int_0^{L_p} f_R dL' = \frac{R_m+R_0}{f_S R_m C_R} \ln(\frac{R_0 + C_R R}{R_0 + C_R R exp(-C_R f_SL_p)})\]

which corresponds to Shuttleworth and Gurney (1990) and Saugier and Katerji (1991) when Sp =

1. The combination of (32) and (34) provides the value of rsc (RSC)

"

SWGE()

Compute ground evaporation rate (mm/d) using Shuttleworth-Wallace with known transpiration.

Ecoshift

" The Shuttleworth-Wallace approach incorporates the energy tradeoff between transpiration and soil evaporation. When transpiration is reduced by low availability of soil water or is zero, BROOK90 uses the new value of transpiration, Ec (ARATE), in subroutine SWGE to get a new value of soil evaporation, Es (ERATE). With Ec known, substituting (5) into (4) and (4) into (2) and solving for Lv ρw E gives

\[(14)\qquad L_v ρ_w E = \frac{R_s}{R_s + R_a}L_v ρ_w E_c + \frac{c_p ρ D_a + Δr_a^s A_s + Δr_a^a A}{R_s + R_a}\]

then

\[(15)\qquad L_v ρ_w E_s = L_v ρ_w E - L_v ρ_w E_c\]

"

SWGRA()

Compute Shuttleworth - Wallace - Gurney Aerodynamic Resistances.

The function computes and returns the aerodynamic resistances $r_a^a$ (source height to reference height), $r_a^c$ (leaf to zero plane), and $r_a^s$ (ground to source height), see Shuttleworth and Gurney (1990).

Ecoshift

" Shuttleworth - Wallace - Gurney Aerodynamic Resistances The three SW aerodynamic resistances, raa, ras, and rac are obtained in subroutine SWGRA by the methods of Shuttleworth and Gurney (1990). The derivations of their equations are not given here.

The friction velocity, u*, is first obtained from the classic logarithmic wind profile as

\[(16)\qquad u^* = \frac{ k u_a }{ \log((z_a-d)/z_{0}) }\]

where ua is the wind speed at the reference height, za, z0 is the surface roughness parameter, d is the zero-plane displacement, and k is the von Karman constant. The roughness parameter is a measure of the turbulence-inducing properties of the surface. The zero-plane displacement, d, arises because the height of the effective canopy surface is above the ground surface that is taken as zero height. Both z0 and d are obtained in subroutine ROUGH . This u* equation strictly only applies for neutral atmospheric stability. Corrections for non-neutral stability are well-known (Brutsaert 1982), but are not usually considered where the objective is to evaluate PE for periods of a day and are not used in BROOK90.

Shuttleworth and Gurney (1990) assume that the classic logarithmic wind profile applies above the canopy and that an exponential profile applies within the canopy. As the canopy becomes sparser, they further assume that the effective source height of the energy fluxes remains at the same height as for a closed canopy, $D_c = z_{0c} + d_c$; these values are obtained in subroutine ROUGH. The ground to source height resistance, ras, is obtained by integrating the exponential eddy diffusivity from 0 to $D_c$ to give

\[(17)\qquad r_a^s = \frac{h \exp(n)}{ n K_h} \left( \exp\left( \frac{-n z_{0g} }{h} \right) - \exp\left( \frac{-n D_c }{h} \right) \right)\]

and the source height to reference height resistance, raa, is obtained by integrating from $D_c$ to $z_a$ as

\[(18)\qquad r_a^a = \frac{ \log( \frac{z_a-d}{h-d} ) }{k u^*} + \frac{h}{n K_h} \left( -1 + \exp\left( n \frac{h-D_c}{h} \right) \right),\]

where n is the extinction coefficient for eddy diffusivity, $z_{0g}$ is the roughness parameter of the underlying ground surface, and the eddy diffusivity at the canopy height, h, is

\[(19)\qquad K_h = k u^* (h-d).\]

The exponential wind profile and the derived K(z) are known to be incorrect in canopies of intermediate leaf area index, but a better model is not yet available (Choudhury and Monteith 1988).

The leaf to zero plane resistance, rac, is assumed by Shuttleworth and Gurney (1990) to be a sum of the individual leaf laminar boundary layer conductances, so

\[(20)\qquad r_a^c = \frac{r_b}{ρt_p L_p} = \frac{1}{ρ t_p L_p} \frac{ \frac{n}{ab} \left( \frac{w}{u_h} \right)^{1/2} }{ \left( 1 - \exp( -\frac{n}{2}) \right) },\]

where n is the extinction coefficient for eddy diffusivity, ab is a constant = 0.01 m/s0.5 (Campbell 1977), ρtp is the ratio of projected leaf area to total leaf surface area, Lp is the projected leaf area index, w is the representative leaf width in m, and uh is the wind speed in m/s at the top of the canopy. (Note that Shuttleworth and Gurney (1990) should have n in the numerator, not the denominator.) The $r_a^c$ equation strictly applies only to flat leaves and needles, not to cylindrical needles; but $r_a^c$ is small when w is small so this inaccuracy is negligible. The canopy-top wind speed, uh, is

\[(21)\qquad u_h = (u^* / k) \log\left( \frac{h-d}{z_0} \right).\]

In the Shuttleworth and Gurney equation, rac goes to infinity as Lp goes to zero, but this neglects the contribution of stems and branches to interception loss, especially when Lp = 0. Although the $r_a^c$ equation is specifically for leaves with width w, subroutine SWGRA uses (ρtp Lp + π Sp) for r tp Lp in equation (20), where Sp is the projected stem area index, defined analagous to Lp as the ratio of projected surface area of stems and branches to ground area (Federer et al. 1996). This assumes that a unit of stem surface has the same influence on rac as a unit of leaf surface. To avoid zero divides when Sp = 0, Lp is prevented from being less than 0.00001.

The Shuttleworth - Gurney separation of rac and rsc at the canopy level rather than at the leaf level is theoretically incorrect. In reality, the leaf boundary layer resistance and the leaf diffusion resistance are in series on each side of flat leaves (Jarvis and McNaughton 1986). The two resistances should be summed over each side of the leaf before integrating over the canopy, but this is too complicated for practical application (Choudhury and Monteith 1988). "

SWPE(AA, ASUBS, VPD, RAA, RAC, RAS, RSC, p_fu_RSS, DELTA)

Compute Shuttleworth and Wallace (1985) transpiration and ground evaporation.

Ecoshift

" Shuttleworth and Wallace (1985) (SW) modified the Penman-Monteith method to account separately for the different water vapor and sensible heat pathways from the soil and from the leaves. Instead of the two resistances of equation (1), rc and ra, SW define five: rsc, raa, rac, ras, and rss. Resistances rsc and rac are in the transpiration pathway while rss and ras are in the soil evaporation pathway and raa is common to both. The canopy surface resistance, rsc, is the resistance to movement of water vapor out of the leaves. The resistance rac restricts vapor movement from the leaf surfaces to the effective source height for water vapor in the canopy. The resistance between the source height and a reference height above the canopy is raa, which corresponds to ra in equation (1). The reference height is that at which air temperature, humidity, and wind speed are known. The resistance to movement of water vapor from inside the soil to the soil surface is rss. The resistance to vapor movement from the soil surface to the source height is ras. The resistances rac, ras, and raa are assumed also to apply to sensible heat transfer.

Shuttleworth and Wallace (1985) start with

\[(2)\qquad L_v ρ_w E = L_v ρ_w E_c + L_v ρ_w E_s\]

where Ec is transpiration and Es is soil evaporation, then write an equation similar to (1) for each term

\[(3)\qquad L_v ρ_w E_c = \frac{Δ(A-A_s) + c_p ρ D_0/r_a^c}{Δ + γ + γ \frac{r_s^c}{r_a^c}}\\ (4)\qquad L_v ρ_w E_s = \frac{Δ(A-A_s) + c_p ρ D_0/r_a^s}{Δ + γ + γ \frac{r_s^s}{r_a^s}}\]

where D0 is the vapor pressure deficit at the effective source height, A is Rn - S or the available energy above the canopy, and As is Rns - S or the available energy at the ground. From the relationships of sensible and latent fluxes to the gradients and resistances, and using the definition of Δ, Shuttleworth and Wallace obtain

(5) D0 = Da + raa [Δ A - (Δ + γ) Lv ρw E ] / cpρ

Algebraic manipulation of (2), (3), and (4) to eliminate D0 leads to:

(6) Lv ρw E = Cc Mc + Cs Ms

where

\[(7)\qquad M_C = \frac{ΔA + \frac{c_p ρ D_a - Δ r_a^cA_s}{r_a^a + r_a^c}}{Δ + γ + \frac{γ r_s^c}{r_a^a + r_a^c}}\\ (8)\qquad M_S = \frac{ΔA + \frac{c_p ρ D_a - Δ r_a^s(A-A_s)}{r_a^a + r_a^s}}{Δ + γ + \frac{γ r_s^s}{r_a^a + r_a^s}}\]

(9) Cc = 1 / { 1 + Rc Ra / [ Rs ( Rc + Ra ) ] }

(10) Cs = 1 / { 1 + Rs Ra / [ Rc ( Rs + Ra ) ] }

with

(11) Ra = (Δ + g) raa

(12) Rs = (Δ + γ) ras + γ rss

(13) Rc = (Δ + γ) rac + γ rsc

Although the algebra appears complicated, these equations for the first time provide a PM-type theory that includes both transpiration and soil evaporation. The total evaporation, E, must be obtained from (6) first, then D0 from (5), before the two components can be calculated from (3) and (4). Subroutine SWPE includes equations (3) through (13), more or less in reverse order.

The outputs Ec (PRATE) and Es (ERATE) from SWPE are in units of mm/d whereas Lvρw E in (1) is output as W m-2 from function PM. The conversion is ETOM * WTOMJ . "

WEATHER()

Compute solar radiation, temperature and wind speed.

Ecoshift

WEATHER includes all adjustments of input weather data, including separation into daytime and nighttime values.

If daily solar radiation (SOLRAD) is input as zero, it is estimated as 0.55 * I0HDAY, or 55% of the potential solar radiation for the doy and location. The 0.55 value is an overall generalization for the United States, where values range from 0.50 in the east to 0.60 in the west (U.S. Department of Commerce 1968). As of Version 4.8 this value can be changed on the BROOK90 main window.

If vapor pressure (EA) is input as zero, it is estimated as the saturated vapor pressure at the minimum daily temperature (TMIN) using subroutine ESAT.

If daily average wind speed at a weather station (UW) is input as zero, it is estimated as 3 m s-1. This is a surprisingly good approximation for most weather stations in the United States at all seasons of the year (U.S. Department of Commerce 1968). For other default values, enter the value as UA for each day in the data file. Note: measured values of zero wind speed should be entered as 0.1.

The average temperature for the day (TA) is taken as the average of the input maximum and minimum temperatures (TMAX and TMIN). Daytime (TADTM) and nighttime (TANTM) average temperatures are calculated by assuming a sine wave variation between TMAX and TMIN. Integration leads to

TADTM = TA + [ (TMAX - TMIN) / (2 π DAYLEN) ] sin( π DAYLEN )

TANTM = TA - { (TMAX - TMIN) / [2 π (1 - DAYLEN)] } sin( π DAYLEN )

where DAYLEN is the solar daylength determined in subroutine SUNDS.

For wind speed, a parameter, WNDRAT, defines the average ratio of nighttime wind speed (UANTM) to daytime wind speed (UADTM). The default value for WNDRAT is 0.3 based on Hubbard Brook data.

UADTM = UA / [ DAYLEN + (1 - DAYLEN) WNDRAT ]

UANTM = WNDRAT * UADTM

where UA is the daily average wind speed at height za.

TA, UA, and the vapor pressure, EA, which is assumed constant over the day, are all theoretically the values at the reference height above the canopy (ZA) . In practice, these values are rarely measured above the canopy of interest, but are usually from some relatively nearby weather station. Any attempt to theoretically adjust TA and EA would require some information on their profiles, such as surface temperature and vapor pressure, which are not known. So BROOK90 assumes that TA and EA are the same at the weather station and at za. However, UA can be estimated from wind speed at the weather station (UW) because wind speed extrapolates to zero at height z0 + d over both surfaces. This adjustment is done in function WNDADJ. UW is prevented from being less than 0.2 m s-1.

Subroutine WEATHER could be modified to do further adjustments to temperature, vapor pressure, and wind. For instance, an elevation adjustment to temperature could be added, using the average environmental lapse rate of -0.6°C / 100m.

WNDADJ(p_fu_ZA, p_fu_DISP, p_fu_Z0, p_FETCH, p_ZW, p_Z0W)

Compute ratio of wind speed at reference height (above canopy) to wind speed at weather station.

Ecoshift

" This function estimates the wind speed (UA) at reference height ZA above the canopy from input wind speed at a remote weather station (UW). Assume that the weather station represents a new surface downwind that has a roughness of z0w (Z0W) and a fetch of F (FETCH). Brutsaert (1982) gives the height of the internal boundary layer, zb, as

\[z_b = 0.334 F^{0.875} z_{0w}^{0.125}\]

For logarithmic wind profiles over both surfaces to have the same wind speed at zb,

\[u_a = u_w \left( \frac{ \log(z_b/z_{0w}) \log((z_a-d)/z_{0}) }{ \log(z_b/z_{0}) \log(z_{w}/z_{0w}) } \right)\]

where zw (ZW) is the height of wind measurement at the weather station (Federer et al. 1996) and d (DISP) is the zero-plane displacement of the canopy. This assumes that the weather station is over a smooth surface so its zero plane displacement can be ignored. If the parameter Z0W is set to zero, then no adjustment is made and ua = uw. "

SNOENRGY(p_fu_TSNOW, p_fT_TA, p_fT_DAYLEN, p_CCFAC, p_MELFAC, p_fT_SLFDAY, p_fu_LAI,
p_fu_SAI, p_LAIMLT, p_SAIMLT)

Compute energy flux density to snow surface, MJ m-2 d-1.

Ecoshift

" The energy flux density to the snow surface (SNOEN, MJ m-2 d-1) is calculated in subroutine SNOENRGY independently of precipitation for the day.

When TA is <= 0°C, SNOEN is the energy used to heat or cool the snowpack

SNOEN = CCFAC * 2 * DAYLEN * (TA - TSNOW)

where CCFAC is an input parameter, and DAYLEN is the daytime fraction of a day. CCFAC is the below-zero degree-day factor for a day with a daylength of 0.5 d. Anderson (1973) pointed out that this degree-day factor appears to vary seasonally. In BROOK90, this seasonality is incorporated by using 2 * DAYLEN as a multiplier in the above equation. BROOK90 is not very sensitive to CCFAC unless it is close to 0; larger values of CCFAC make the snow melt later because the snowpack cools more. When CCFAC = 0, snow temperature is always 0°C and there is never any cold content.

When TA is greater than 0°C, energy is provided that can melt snow, and SNOEN is calculated differently. The energy supply rate, SNOEN, is then

SNOEN = MELFAC * 2 * DAYLEN * SLFDAY * TA * exp(-LAIMLT * LAI - SAIMLT * SAI)

where MELFAC is the melting degree-day factor for a day with a daylength of 0.5 d and no canopy, SLFDAY is the ratio of potential insolation on the slope to that on a horizontal surface, and the input parameters LAIMLT and SAIMLT express the dependence of SNOEN on leaf area index (LAI) and stem area index (SAI). MELFAC uses 0°C as a base and is zero for TA below this. Inclusion of SLFDAY in the MELT equation arises from an assumption that radiation melt plays an important role. If this is not so, then SLFDAY multiplier could be omitted, and snowmelt would not depend on slope-aspect. The functional forms of the SAI and LAI dependencies are based on the somewhat arbitrary curves used by Federer and Lash (1978). "

SNOFRAC(p_fT_TMAX, p_fT_TMIN, p_RSTEMP)

Separate RFAL from SFAL.

Ecoshift

" Separation of daily precipitation into snow or rain is a major problem in hydrologic modeling. For instance, if the wrong precipitation form is chosen in December, simulated streamflow from that day's precipitation could be shifted from December to April or vice versa, a rather significant effect! BROOK90 uses both daily maximum (TMAX) and daily minimum (TMIN) temperatures to allow days on which mixed rain and snow falls. This reduces the potential error from making the wrong choice. The algorithm seems to have been stated first by Willen and Shumway (1971). When TMAX for the day is greater than the parameter RSTEMP and TMIN is less than RSTEMP, the fraction of precipitation as snow, SNOFRC, is

SNOFRC = (RSTEMP - TMIN) / (TMAX - TMIN)

where RSTEMP is the "base" temperature for the rain-snow transition. When TMAX < RSTEMP, SNOFRC = 1; when TMIN > RSTEMP, SNOFRC = 0. The default value of RSTEMP is -0.5°C because that seems to work best at Hubbard Brook. If precipitation is input more than once a day, the same SNOFRC is used for all precipitation intervals. "

SNOVAP(p_fu_TSNOW, p_fT_TA, p_fT_VAPPRES, UA, p_fu_ZA, p_fu_HEIGHT, p_fu_Z0, p_fu_DISP,
p_fu_Z0C,p_fu_DISPC, p_fu_Z0GS, p_LWIDTH, p_RHOTP, p_NN, p_fu_LAI, p_fu_SAI, p_KSNVP)

Compute potential snow evaporation (mm/d).

Ecoshift

" Evaporation rate from the snowpack, and its negative, condensation, are evaluated using the aerodynamic flux equation

E = (cp ρ / g Ls rw ) (e0 - ea) / (raa + ras)

where ea is the vapor pressure of the air, e0 is the surface vapor pressure, and raa and ras are the Shuttleworth-Wallace resistances described in section PET. The constants cp ρ (CPRHO), γ (GAMMA), and the latent heat of sublimation Ls ρw (LS) are constant. BROOK90 assumes that the snowpack is always isothermal and that its temperature does not change diurnally. When the snowpack temperature (TSNOW) is less than 0°C, the surface vapor pressure, e0, is the saturated vapor pressure at TSNOW and is obtained by calling subroutine ESAT. When TSNOW is 0°C, e0 is 0.61 kPa; use of Ls instead of the latent heat of vaporization, Lv, in this case is slightly wrong. The vapor pressure ea (EA) is the input vapor pressure at reference height za . The resistances raa and ras are obtained from subroutine SWGRA using the daily average wind speed (UA). The value of E in mm/d returned from SNOVAP is called PSNVP because it can be reduced in SNOWPACK if snow disappears.

For evaporation from snow in the open, U.S. Army Corps of Engineers (1956) gives

E = 1.9 ua (e0 - ea)

for evaporation in the open. This yields E = 0.6 mm when ua = 3 m/s and the vapor pressure difference is 0.1 kPa. The two E equations are the same when ras = 0, za - d = 2 m, and z0 = 1 mm. Subroutine SWGRA thus gives the appropriate raa for snow in the open.

Colbeck et al. (1979) state "Evaporation from the snow in a forest has received a great deal of attention, with many investigators concluding that it is small.... Although there are many reports of high evaporative losses from forests, these have not been verified from heat balance considerations." Generally, literature values are around 0.5 mm/d in the open, with monthly values around 10 mm and annual values of 20 mm or more. Anderson (1976) gives 15-20 mm annually for Sleepers River, VT.

The modified Shuttleworth and Gurney resistance formulations in subroutine SWGRA for a leafless tall canopy give raa about 3 s/m and ras about 40 s/m for a weather station wind speed of 3 m/s. A common vapor pressure difference of 0.1 kPa then gives a very high evaporation of 1.6 mm/d or 50 mm/ month. The problem may be that the resistances calculated by SWGRA are too small, either because of too large a roughness parameter for leafless deciduous forests, or because stability effects are ignored. To fix the problem BROOK90 includes KSNVP, which is an arbitrary constant multiplier of PSNVP. Use KSNVP = 1.0 for open open areas or short canopies normally covered by snow; but for forest, KSNVP of 0.3 gives more reasonable values of SNVP. More work is obviously needed on amount and theory of snow evaporation under forests.

Note that although evaporation and condensation of water are simulated in SNOVAP, the accompanying latent transfer is not simulated. The snow energy balance in subroutine SNOENRGY is (unfortunately) decoupled from the snow evaporation-condensation process. "

SNOWPACK(p_fu_RTHR, p_fu_STHR, p_fu_PSNVP, p_fu_SNOEN, u_CC, u_SNOW,
u_SNOWLQ, p_DTP, p_fT_TA, p_MAXLQF, p_GRDMLT,
p_CVICE, p_LF, p_CVLQ)

Update status of snowpack. Generally it computes in terms of mass balance: uSNOWnew = uSNOW + pDTP * (pfuSTHR + auxduRSNO - auxduSMLT - auxduSNVP)

However, there are numerous checks involved regarding availability of snow and energy that modify the above mentioned mass fluxes. The subroutine returns the magnitude of these fluxes in [mm/day] enabling the separate use of them for the isotope mass balance.

Ecoshift

" In each precipitation interval, throughfall of rain (RTHR) and throughfall of snow (STHR) are calculated and subroutine SNOWPACK is entered if there is STHR or if there is SNOW. This subroutine adds throughfall to the snowpack, subtracts groundmelt, snow evaporation, and drainage of liquid water, and calculates the new cold content (CC) and liquid water content (SNOWLQ) of the snowpack. There are a number of different ways to program all this adding and subtracting of energy and water. The program flow of SNOWPACK has been selected to make the many realizable situations as clear as possible. Much shorter algorithms could have been used, but at the expense of clarity. Any alterations to this routine must be made carefully, keeping all the possibilities in mind. However, unlike routine SNOENRGY, the content of the SNOWPACK routine is essentially standard for all snow models and alterations generally should not be necessary.

In SNOWPACK, snow throughfall (STHR) is first added to SNOW. If TA is <0°C, the cold content of the new snow is increased by

CC = CC - CVICE * TA * STHR * DTP

where DTP is the precipitation interval time step and CVICE is the volumetric heat capacity of ice. If this addition of cold snow causes both cold content (CC) and liquid water (SNOWLQ) to coexist, then liquid water is refrozen; CC is used to refreeze part or all of the liquid.

Groundmelt (GRDMLT) and snow evaporation-condensation are dealt with next. In the precipitation time interval (DTP), the fraction of the snowpack that melts by groundmelt and evaporates is

FRAC = (GRDMLT + PSNVP) * DTP / SNOW.

If FRAC is > 1 then all the snow melts and evaporates at rates proportional to GRDMLT and PSNVP. If FRAC is < 1, SNVP is equal to potential snow evaporation (PSNVP) and GRDMLT drains from the snowpack as snowmelt (SMLT). An assumption is made that evaporation and groundmelt remove any liquid water and cold content that is associated with the snow removed, so SNOW, SNOWLQ, and CC are all reduced by 1 - FRAC. If PSNVP is negative, condensation is occurring; SNOW, SNOWLQ, and CC are correspondingly increased. This is simple, but not quite accurate. If no snow is left the routine ends.

The amount of snowpack warming or cooling is calculated next. The equivalent amount of ice melted by the energy input from SNOEN and the heat included in any warm rain is

EQEN = DTP * (SNOEN + RTHR * RMAX(TA, 0) * CVLQ) / LF

where CVLQ is the specific heat of water, and LF is the latent heat of fusion of water. Both CVLQ and LF are constant. When EQEN is less than 0, the snow is cooling; first any SNOWLQ is refrozen, then CC is increased. CC is not allowed to be reduced so that TSNOW is below the mean daily air temperature, TA, although it may remain colder, i.e. if TA < 0°, CC can only be reduced to TSNOW = TA. When EQEN is greater than 0, the snow is warming; first CC is reduced, then SNOWLQ is produced, and finally melt (SMLT) occurs.

Finally, any rain throughfall (RTHR) is added to the snowpack. If any CC exists it refreezes rain until the CC is "used up". Any additional rain then increases SNOWLQ; when the maximum SNOWLQ is reached, the input of rain to the snow (RSNO) has also reached its maximum. In all cases the final results are a new SNOW, new SNOLQ and CC, and a value of RSNO.

MSBPREINT then calculates the rain passing through the snow (RNET) as

RNET = RTHR - RSNO.

When SNOW exists at the beginning of the day, soil evaporation (SLVP) is zero. "

INTER(p_fT_RFAL, p_fu_PINT, p_fu_LAI, p_fu_SAI, p_FRINTL, p_FRINTS, p_CINTRL, p_CINTRS, p_DTP, u_INTR)

Compute rain catch rate (interception) and evaporation rate of intercepted rain in mm/d.

Rain interception, used when p_NPINT > 1.

Arguments

• p_fT_RFAL: rainfall rate, mm/d
• p_fu_PINT: potential interception rate, mm/d
• p_fu_LAI: projected leaf area index, m2/m2
• p_fu_SAI: projected stem area index, m2/m2
• p_FRINTL: intercepted fraction of pfTRFAL per unit pfuLAI
• p_FRINTS: intercepted fraction of pfTRFAL per unit pfuSAI
• p_CINTRL: maximum interception storage of rain per unit pfuLAI, mm
• p_CINTRS: maximum interception storage of rain per unit pfuSAI, mm
• p_DTP: precipitation interval time step, d
• u_INTR: intercepted rain storage, mm,

Ecoshift:

Older studies of rain and snow interception regressed throughfall on precipitation, but such interpretation ignored the fact that energy supply rather than water supply may limit interception and also ignores storm duration/intensity and interstorm interval. Detailed simulation models of rain interception over time through a storm have been developed (Rutter et al. 1972) but these are too complex to include in a hydrologic model like BROOK90. Much less is known about the complicated process of evaporation of intercepted snow.

Subroutine INTER accounts in the simplest way for the concepts of catch rate, evaporation rate, and canopy capacity. The same algorithm is applied to both rain and snow, which are considered to behave independently with respect to their interception. INTER is used when precipitation data are input more than once a day in a precip. interval file (PINT > 1). If only daily precipitation is input, then the modified procedure of INTER24 is used.

The conservation of mass equation for rain interception can be written as

dS/dt = C - I - D

where S is the amount of water stored on the canopy (mm), C is the catch rate, or rate of water input to the canopy, I is the rate of evaporation of intercepted water, and D is the drip rate, or rate of transfer of liquid water to the ground. The same equation applies to snow or mixed snow and rain, when any solid-liquid phase change is ignored and D includes all rain or snow blowing or falling from canopy to ground.

BROOK90 ignores D by defining C as a net catch rate (C - D), or only the portion of the catch that will sooner or later evaporate, so, from the Flow Chart,

d INTR / dt = RINT - IRVP

for rain, and

d INTS / dt = SINT - ISVP

for snow, where INTR and INTS are the canopy storages, RINT and SINT are the net catch rates, and IRVP and ISVP are the evaporation rates, for rain and snow respectively.

BROOK90 assumes that interception catch rates, RINT and SINT, are a constant fraction of rainfall or snowfall until the canopy reaches a storage capacity. Until the capacity is reached, RINT and SINT are assumed to be linear functions of LAI and SAI, so that

RINT = (FRINTLAI * LAI + FRINTSAI * SAI) * RFAL

and

SINT = (FSINTLAI * LAI + FSINTSAI * SAI) * SFAL

where RFAL and SFAL are rainfall rate and snowfall rate as determined from subroutine SNOFRAC, FRINTLAI and FSINTLAI are the catch fraction per unit LAI for rain and snow, respectively, and FRINTSAI and FSINTSAI are the catch fraction per unit SAI for rain and snow, respectively.

The canopy has capacities or maximum values of INTR and INTS that depend on LAI and SAI. In BROOK90 these dependencies are assumed linear. The parameters CINTRL and CINTRS are the capacities for intercepted rain per unit LAI and SAI respectively, so that INTRMX, the capacity for rain, is

INTRMX = CINTRL * LAI + CINTRS * SAI.

For snow,

INTSMX = CINTSL * LAI + CINTSS * SAI,

and the capacity parameters are generally larger than for rain. The eight interception parameters, FRINTLAI, FRINTSAI, FSINTLAI, FSINTSAI, CINTRL, CINTRS, CINTSL, and CINTSS, only control interception loss in small storms; interception loss in large storms is controlled by the evaporation rate of intercepted water (PINT) and the storm intensity and duration.

The rate at which intercepted water evaporates (PINT) is calculated from the Shuttleworth-Wallace equations by calling subroutine SWPE (Section PET) with the canopy resistance rc = 0. The soil surface resistance (RSS) is not reduced for the PINT calculation. The MSBDAYNIGHT routine does this separately for daytime and for nighttime weather variables and the results are weighted by daylength (DAYLEN) to produce PINT. PINT is considered to be constant throughout the daily time step; its actual diurnal variation is ignored.

The canopy is considered to be either completely wetted or completely dry. Partial canopy wetting and drying is not treated in BROOK90, though it is a key component of specific models of the interception process (Rutter et al. 1972). Subroutine INTER determines the actual catch rate (RINT or SINT) and the actual evaporation rate (IRVP or ISVP) for the precipitation time step in the three cases that the canopy dries during the timestep, the canopy wets but does not reach capacity, and the canopy reaches capacity. The routine appropriately handles the case of a wet canopy with decreasing capacity because of decreasing LAI or SAI by allowing RINT or SINT to be negative.

INTER(p_fT_RFAL, p_fu_PINT, p_fu_LAI, p_fu_SAI, p_FRINTL, p_FRINTS, p_CINTRL, p_CINTRS, p_DTP, u_INTR)

Compute rain catch rate (interception) and evaporation rate of intercepted rain in mm/d.

Rain interception with duration in hours, used when p_NPINT = 1. Same routine is used for snow interception, with different calling variables.

Arguments:

• p_fT_RFAL: 24-hour average rainfall rate, mm/d
• p_fu_PINT: potential interception rate, mm/d
• p_fu_LAI: projected leaf area index, m2/m2
• p_fu_SAI: projected stem area index, m2/m2
• p_FRINTL: intercepted fraction of pfTRFAL per unit pfuLAI
• p_FRINTS: intercepted fraction of pfTRFAL per unit pfuSAI
• p_CINTRL: maximum interception storage of rain per unit pfuLAI, mm
• p_CINTRS: maximum interception storage of rain per unit pfuSAI, mm
• p_DURATN: average storm duration, hr
• u_INTR: intercepted rain storage, mm,
• MONTHN: Month of the year

Ecoshift:

" Subroutine INTER24 - daily interception Proper representation and integration of the interception process is a problem for hydrologic models that use a daily interval for precipitation input (p_NPINT = 1), because the storm duration is not known. For a brief, intense storm, the canopy wets once and the interception loss is limited primarily by canopy capacity. For a low intensity, all day storm, the canopy stays wet and the interception loss is limited primarily by the potential interception, PINT. This problem is worst when only daily precipitation is known, and decreases as precipitation is given at shorter intervals.

Subroutine INTER24 was developed because the use of subroutine INTER for daily precipitation consistently produced too much interception. INTER24 is a modification of INTER that loops through the procedure every hour,using the PINT rate for each hour. DURATN is a parameter that specifies the average hourly duration of precipitation for each month of the year. INTER24 truncates DURATN to the next lower even integer, and then centers the "storm" on noon. Thus if DURATN is input as 7.5, the daily precipitation is assumed to occur at a constant rate from time 0900 to 1500. Centering on noon is only used to see how much interception carries over into the next day. The algorithm for each hourly loop is the same as for INTER, except that rates are in mm/hr and amounts are summed over the day. The interception catch rate (RINT or SINT), and the evaporation rate (IRVP or ISVP) are returned to MSBPREINT as average rates over the day.

To determine appropriate values of DURATN I examined hourly precipitation data for 4 years (one year at Hubbard Brook) from the SAMSON data set. Averaging the number of hours per day of precipitation of 0.02 inch (0.5 mm) or greater over days with such precipitation gave the following results after a little smoothing

                    J   F   M   A   M   J   J   A   S   O   N   D
San Juan PR         3   2   2   2   2   2   2   3   3   3   3   3
Atlanta GA          5   5   5   5   4   4   3   3   4   4   5   6
Caribou ME          4   4   5   5   4   4   4   4   4   6   6   5
Madison WI          4   4   5   3   3   2   3   3   4   4   5   5
Lake Charles LA     5   4   3   3   3   3   2   2   3   3   4   5
Phoenix AZ          4   4   4   4   4   2   2   2   2   2   4   4
Rapid City SD       3   3   3   4   4   3   2   2   2   2   4   4
Tacoma WA           6   6   5   4   4   4   4   4   4   4   6   6
Fairbanks AK        3   3   4   4   4   4   3   3   4   4   4   3
Hubbard Brook NH    5   5   5   4   4   4   4   4   4   5   5   5

Apparently a default DURATN of 4 hours is appropriate for all months anywhere in the U.S. "

PLNTRES()

Allocates total plant resistance to xylem and root layers.

Ecoshift: Subroutine PLNTRES is called at the beginning of each day to obtain resistivities to liquid water flow: rhizosphere resistivity for each soil layer, root resistivity in each soil layer, and xylem resistivity. These parameters, together with soil water potential in each layer (PSITI) and critical plant water potential (PSICR) control the supply of water to transpiring leaves and thus the reduction of actual transpiration below potential transpiration. As defined by Hunt et al. (1991) the resistances used here are "potential difference resistivities", because the transpiration flux rate is in units of mm/d and the potential gradient is in MPa. The resistivities have units of MPa d mm-1.

In most of this subsection and the following subsection (TBYLAYER), algebraic notation is used instead of variable names. The correspondence is:

rri RROOTI() ψt PSIT P PTR fi RTFRAC rx RXYLEM ψti PSITI() T ATR Di D() rp RPLANT ψc PSICR Ti ATRANI() Li RTDENI ri RI S SUPPLY δi DELT rt RT αi ALPHA() The following additional parameters or constants are input to the routines

fx FXYLEM di RELDEN ρwg RHOWG R1 RTRAD d DISPC π PI Lr RTLEN Ki KK()

Several other algebraic variables occur in the derivations below, but are not needed in the program.

If soil water potential is uniform through the root zone and rhizosphere resistivity is negligible, a bulk plant resistivity, rp, can be defined by

rp = ( ψt - ψ - ρw gd ) / T'

where T' is the transpiration rate, ψt is the soil water potential, and ψ is the leaf water potential. The ρwgd term is the gravity potential difference above the ground surface, where ρwg is the density of water times the acceleration of gravity, and d is the effective height of canopy evaporation, taken here as the closed canopy zero-plane displacement (DISPC). Change in water storage within the plant is ignored in BROOK90, as Hunt et al. (1991) have shown that it does not matter to total daily transpiration.

In BROOK90, rp is determined primarily by the maximum bulk plant conductivity (MXKPL), which is an input parameter. MXKPL is the water uptake rate for a closed canopy per unit of soil to leaf potential difference. When the soil is wet so ψt is close to 0, many plants have a ψ of around -1.5 MPa when transpiration rate is about 0.5 mm/hr, a normal midday rate for a sunny day in temperate regions. This gives a MXKPL of 8 mm d-1 MPa-1. Abdul-Jabbar et al. (1984) found literature values to range only from 7 to 30 mm d-1 MPa-1, which seems a surprisingly narrow range for plant canopies of any species. Actual plant resistivity, rp, is determined in subroutine CANOPY.

Fig. EVP-1. Resistances and potentials in the liquid flow path for transpiration, for 3 layers with roots. ψ is the leaf potential, ψx is potential at gorund level, and ψti are total soil water potentials. rₓ is xylem resistance, rᵣᵢ are root resistances, and rₛᵢ are rhizosphere resistances.

Figure EVP-1 shows the resistance network used in BROOK90. Each soil layer is considered to have a rhizosphere resistivity, rsi, and a root resistivity, rri, in series. Each layer is considered in parallel with the others. An additional resistivity in series accounts for resistance to flow through the xylem above ground level, rx. The total soil-water potential, ψti, differs among layers. The leaf water potential, ψ, is assumed constant through the canopy. The potential at the ground surface is ψx. This system and the parameterization of it was developed by Federer (1979) and Hunt and others(1991).

The xylem resistance, rx, is

rx = fx rp ;

where fx (FXYLEM) is an input parameter, which is probably close to zero for short canopies, but may be 0.5 or higher for forests (Hunt et al. 1991). The use of fx here specifies the amount of rp (RPLANT) that should not be divided among soil layers, so it is effectively the fraction of the plant resistance that is above ground level.

The plant resistance that is below ground level, rp - rx, is the parallel combination of the individual layer resistances, rri (Fig. EVP-1). BROOK90 assumes that the root resistivity per unit length of root is constant, so rri (RROOTI) is inversely proportional to the fraction of total root length that is in each layer, fi

rri = ( rp - rx ) / fi .

BROOK90 parameterizes root distribution by the relative density of roots in each soil layer di (RELDENi). RELDENi is obtained from ROOTDEN in subroutine RTDEN . Then fi is

fi = di Di / S ( di Di )

where Di is the stone-free layer thickness, THICK(i) * (1 - STONEF(i)).

When the soil is at uniform water potential and rhizosphere resistances, rsi , are negligible and fx = 0, then the relative transpiration withdrawal from each layer is proportional to fi (see subroutine TBYLAYER). Increasing fx makes the withdrawal less dependent on fi .

Usually the rhizosphere resistance only becomes significant when the soil is dry. Following Federer (1979) and Cowan (1965), the rhizosphere resistance, rsi , calculated in subroutine TBYLAYER is rsi = αi / Ki , where Ki is the hydraulic conductivity of the rhizosphere. The value of αi (ALPHAi) for a layer is

αi = Ai / ρwg Di

where Ai from Cowan (1965) is

Ai = ( 1 / 8 π Li ) [ δi - 3 - 2 ( ln δi ) / ( 1 - δi ) ]

where Li is the root density in the layer (mm/mm3) and δi is the root volume fraction in the layer, obtained as

δi = π R12 Li

where R1 is the average radius of the absorbing roots (RTRAD) , which is an input parameter. (Note that δi has no relation to the daylength, δ ). The root density is

Li = fi Lr / Di

where Lr is the total length of absorbing roots per unit ground area in mm/mm2. The value of Lr in m/m2 is found in subroutine CANOPY as the parameter MXRTLN reduced by DENSEF and RELHT. The dependence on the seasonal RELHT assumes that root length increases proportionally with height growth. MXRTLN is an input parameter expressing the length of absorbing roots per unit ground area in m/m2. Lr and R1 only affect rhizosphere resistance and thus are only important for dry soil or when Lr is small. "

TBYLAYER()

Allocate transporation among soil layers.

Ecoshift: TBYLAYER determines the rate at which liquid water can be supplied to transpiring leaves, compares this rate with the potential transpiration rate, sets actual transpiration equal to the lesser of the two, and then allocates the transpiration among soil layers. This routine is based on the model of Federer (1979), which has been widely used, e.g. Wetzel and Chang (1987), Levine and Salvucci (1999).

The routine requires iteration when outflow from roots is prevented for layers in which xylem potential ψx is greater than soil water potential ψti(Fig. EVP-1).

The resistance to uptake from each layer, ri, is

ri = rri + rsi = rri + αi / Ki

(Fig. EVP-1) where rri (RROOTI) and αi (ALPHAi) are variables from subroutine PLNTRES and Ki is the rhizosphere conductivity of the layer (Cowan 1965, Federer 1979). Estimating rhizosphere conductivity requires an iterative solution as described in Federer (1979), so BROOK90 uses the Ki (KKi) of the bulk soil instead. This tends to underestimate the rhizosphere resistance, but the error is probably unimportant unless rp is quite small. The total root resistance to uptake is

rt = 1 / Σ ( 1 / ri ) .

The transpiration rate for each layer, Ti , is

Ti = ( ψti - ψx ) / ri

where ψti is the total water potential for the layer (PSITI) , and ψx is the xylem potential. The total transpiration from all layers is

T = ΣTi = Σ( ψti / ri ) - ψx Σ ( 1 / ri ) = ( ψt - ψx ) / rt

where ψt is a weighted mean soil water potential defined as

ψt = rt Σ (ψti / ri ).

Eliminating ψx from the Ti and T equations gives

Ti = ( ψti - ψt + rt T ) / ri

which is the equation used to distribute T among layers.

In BROOK90 actual transpiration rate, T, is the lesser of potential transpiration rate, P, and the maximum possible rate of water uptake from the soil, S. Following Federer (1979) and Lynn and Carlson (1990), BROOK90 assumes that the S obtains when the leaf potential, ψ, is the critical potential at which stomates close, ψc , which is an input parameter (PSICR). The assumption that the T is the lesser of the P and S is equivalent to an assumption that stomatal opening is not limited by leaf water potential, ψ, until the critical potential, ψc, is reached. The stomata are then assumed to close as much as necessary to maintain ψ = ψc. This abrupt switch from a demand-limited to supply-limited regime is oversimplified, but is convenient for modeling. This type of behavior induces a flat-topped diurnal transpiration curve (Lynn and Carlson 1990). The critical potential, ψc, represents the maximum suction that the plant can exert to get water from the soil, the minimum soil-water potential that the plant can induce, and the lower limit of soil water availability.

The water supply rate, S, is the potential difference between leaf water potential, ψ, and the xylem potential, ψx, divided by the xylem resistance, rx , when ψ is equal to the critical potential, ψc (see Fig EVP-1). With allowance for the gravity potential between the ground and the canopy zero-plane displacement height, d (DISP), this is

S = ( ψx - ψc - ρw g d ) / rx

S below ground must be T = ( ψt - ψx ) / rt from above, and equating the two S equations to eliminating ψx, gives

S = ( ψt - ψc - ρw g d) / ( rt + rx )

which is assumed to be constant throughout a day.

Fig. EVP-2. Transpiration for the day, T' (shaded area), as the lesser of a constant water supply rate, S, and a half-sine potential transpiration, P'.

Following Federer (1982), BROOK90 assumes that the daytime potential transpiration rate, P', varies as a half sine wave. The actual transpiration rate, T' , is the lesser of S and P' (Fig. EVP-2). The average value of T' over the daytime, T, is

T = P ( 1 + R cos-1 R - sin (cos-1 R) ) R < 1

T = P R >= 1

where

R = 2 S / π P

where P is the average of P' over the daytime. (Notation here differs slightly from Federer (1982) who used S' for S, D/d for P, and T/d for T where d is daylength.) At night, P is assumed constant and T is the lesser of S and P.

Normally, all layers with roots are included in the above calculations. However, when some layers are wet and others are dry, it is possible for ψti (PSITI) in one or more layers to be lower than ψx so that the roots in those layers are releasing water to the soil. The question of whether such outflow occurs or is somehow prevented by the roots is controversial (Hunt et al. 1991). More recent work shows that outflow does occur (Dawson 1993). Richards and Caldwell (1987) name the process "hydraulic lift" because it moves water upward from wetter deeper soil layers through the plant to shallow drier layers. In BROOK90 outflow from roots is prevented when the parameter NOOUTF is set to 1, and is allowed if NOOUTF is 0. In general, BROOK90 usually will work better when NOOUTF = 1; this is the default. When NOOUTF = 1 and any Ti is negative, the layer with the most negative Ti is eliminated and new values of rt , ψt , T, and Ti are obtained. If any Ti are still negative, the elimination process is repeated. This elimination procedure causes transpiration from a layer to cease when its potential is still greater than PSICR.

compute_isotope_U_of_INTS_INTR_SNOW_and_SLFL(
    p_δ2H_PREC, p_δ18O_PREC, p_fT_TADTM, p_fT_VAPPRES,
    # for INTS (in: SINT; out: ISVP):
    u_INTS, aux_du_SINT, aux_du_ISVP, p_DTP, u_δ2H_INTS, u_δ18O_INTS,
    # for INTR (in: RINT; out: IRVP):
    u_INTR, aux_du_RINT, aux_du_IRVP, u_δ2H_INTR, u_δ18O_INTR,
    # for SNOW (in: STHR, RSNO (both δ_PREC); out: SMLT, SNVP (δ_SNOW and fractionated)):
    p_fu_STHR, aux_du_RSNO, aux_du_SMLT, aux_du_SNVP, u_δ2H_SNOW, u_δ18O_SNOW,
    # to compute isotopic signature of soil infiltration: SLFL
    p_fu_RNET)

Computes updated values of states INTS, INTR, and SNOW as well as their isotopic composition. Compute mixing and evaporative fractionation, and also compute the isotopic composotion of the resulting flux that infiltrates into the soil: pfuδ18OSLFL, pfuδ2HSLFL

The function is called in the daily callback.

The function returns: pfuδ18OSLFL, pfuδ2HSLFL, as well as: uINTS uδ18OINTS uδ2HINTS uINTR uδ18OINTR uδ2HINTR uSNOW uδ18OSNOW uδ2H_SNOW