Code Listing

Isotopic mixing and fractionation

This module defines constants and functions associated with mixing and isotopic fractionation processes simulated by LWFBrook90.jl.

States INTS, INTR, SNOW are updated once per day. The isotopic composition at the end of the day is computed based on an operator splitting approach. First, non-fractionating in- and outflows are considered for each compartment and a resulting isotopic composition is computed. Second, the fractionation associated with the evaporation process is computed, assuming a Craig-Gordon fractionation process applied to a reservoir with the evaporation process as single sink. The reservoirs inital condition is the isotopic composition resulting from the mixing computed in the first operator step.

Turbulence conditions are fixed to different values for each compartment (XINTR, XINTS, X_SNOW).

Evaporation from the soil is modelled from the uppermost layer only.TODO(bernhard): Should this be improved? Turbulence conditions for evaporation from the soil layer are conditional on the soil moisture status $θ$ following (Zhou-2021-EnvironModelSoftw (X is called nk there)) as a weighted average between Xs = 1$ (molecular diffusion only) and $X_s = 0.5$ (both molecular and turbulent diffusion):

\[X_{soil} = \frac{(θ-θ_{res})*X_a + (θ_{sat}-θ)*X_s}{(θ_{sat}-θ_{res})}\]

Enrichment of xylem water due to evaporation from the stomatas is not modelled.

Soil water properties

Text copied from Ecoshift on module KPT:

" This section describes the soil water properties in BROOK90, and the four routines used to calculate their values. The section also discusses the concept of field capacity, and how it is or is not used. Finally the section comments on the selection of soil water parameters and the use of the menu item EditParameters-KPTGraph. This section uses both the standard algebraic notation for soil water variables as in Clapp and Hornberger(1978) and the BROOK90 variable names. The correspondence is θ = THETA W = WETNES ψ = PSIM K = KK b = BEXP θf = THETAF Wf = WETF ψf = PSIF Kf = KF n = CHN θs = THSAT Wi = WETINF ψi = PSIINF Ks = KSAT m = CHM T = THICK S = STONEF

Ecoshift-KPT: Functional relationships among soil-water content, θ, matric potential, ψ, and hydraulic conductivity, K, are required for any simulation that moves water through the soil. There is a rather vast literature on this subject, though much of it emphasizes agricultural soils rather than natural soils, which are less disturbed and higher in organic matter.

Clapp Hornberger (FLAG_MualVanGen=0) see: http://www.ecoshift.net/brook/kpt.html

BROOK90 uses a modification of the Campbell (1974) expressions with the near-saturation interpolation of Clapp and Hornberger (1978).

1. W = (θ - θr) / θs - θr)
2. Ψ = Ψs*W^(-b) with b=1/λ
3. K = Ks*W^(2b+3) with b=1/λ

near saturated:

1. Ψ = -m(W - n)(1 - W) = (-mW + mn) +mW^2 -Wmn = mW^2 -m(1+n)W + m*n

2.  m = -Ψi/(1 - Wi) * [1/(1-Wi) - b/Wi]

3.  n = 2 Wi - 1 + (b*Ψi / (m*Wi))  ### where Ψi is Ψ at Wi obtained from 2)

clearly unsaturated:

1. Ψ = Ψf(W/Wf)^(-b) = ΨfWf^(b)*W^(-b)

2. K = Kf*(W/Wf)^(2b+3) ### where Wf = θf/θs is the wetness at θf (with θr=0)

   ==> dψdW = ΨfWf^(b) (-bW^(-b-1)) # in clearly unsaturated range ==> dψdW = 2mW - m(1+n)

Mualem van Genuchten (FLAG_MualVanGen=1)

1. W = 1/(1 + (αh)^n )^m (Radcliffe eq.2.47)
2. ψ = -1/α * [ (1-W^(1/m))*(W^(-1/m)) ]^(1/n) (Radcliffe p.122)

...

1. K = KsW^l[ 1 - (1-W^(1/m))^m ]^2

   ==> dψdW = .....

"

Water movement in soil

Text copied from Ecoshift on module WAT:

" BROOK90 does not try to account for all possible paths of water movement through soil, which is a very complicated subject (McDonnell 1990). As a lumped parameter model it does not move water laterally between subareas. However, it does simulate several different pathways of water movement over or through the soil to streamflow and groundwater (see Flow Chart):

snowmelt and throughfall on impervious areas goes directly to streamflow (SRFL) snowmelt and precipitation on variable saturated source areas goes directly to streamflow (SRFL) remaining snowmelt and streamflow (SLFL) infiltrates either all to the surface layer or to several layers via vertical pipes or macropores (INFL) some infiltrated water can go directly to streamflow as "bypass" flow in pipes (BYFL) classic vertical matric flow between soil layers (VRFLI) lateral or downslope movement of matric water to streamflow (DSFL) vertical drainage of matric water to groundwater (VRFLN) discharge of groundwater to streamflow (GWFL) deep seepage loss from groundwater (SEEP)

SLFL, INFL, VRFLI, and VRFLN are internal flows. SRFL and BYFL produce streamflow only on the day of precipitation, simulating a streamflow response of less than 24 hours duration. Users generally should choose one or the other of these flows. DSFL produces a response over several days only if VRFL is limited by a low conductivity layer within the profile or VRFLN is limited by setting DRAIN < < 1. VRFL from the bottom of the profile produces a response of several days only if there is no groundwater. GWFL reponse can vary from several to many days. SEEP produces no streamflow at all.

If the surface horizon becomes saturated, infiltration-excess overland flow is simulated as BYFL from the top layer. The input rate (SLFL) is constant over the precipitation interval, which may be a full day for most users. So with saturated hydraulic conductivities usually 200 mm/d or more, such overland flow will be rare.

Vertical water movement out of a layer is a combination of matrix flow, VRFL(I), and the macropore infiltration, SLFL(I). SLFL(I) either becomes BYFL from deeper layers or becomes soil water in deeper layers. VRFL(I) will generally increase with depth as SLFL(I) decreases. Lysimeter users may want to assume that suction lysimeters collect VRFL(I) whereas zero-tension lysimeters collect both VRFL(I) and SLFLI(I).

In this section a subscript i is used to indicate individual layers. "

Radiation

Text copied from Ecoshift on module SUN:

" Daily solar radiation on a horizontal surface (SOLRAD in MJ/m2) is an input variable to BROOK90. This value is sometimes called global radiation to emphasize that it includes both direct or beam radiation from the sun and diffuse radiation from the sky hemisphere. It is directly measured by various types of pyranometers at most research sites and some locations of the National Weather Service. SOLRAD is internally prevented from being larger than 0.99 times the potential insolation (in the absence of an atmosphere), I0HDAY.

If no data are available and input SOLRAD is zero, BROOK90 uses a fixed fraction of the potential solar radiation on a horizontal surface. This fraction had been set at 0.55, but as of Version 4.8 it can be changed on the BROOK90 main window. Use of this generalized fraction loses the effects of day-to-day variation of solar radiation in the model.

Alternatively, there are various simple to complicated methods for estimating SOLRAD from cloud cover measurements, which are more widely made (Brutsaert 1982, Dingman 1994).

Another approach to estimating daily solar radiation uses the fact that the daily temperature range is usually larger on clear days than on overcast ones. The equation of Bristow and Campbell (1984) has been used worldwide in various modifications. For Hubbard Brook I use

R / I0 = a (1 - exp(-b Δ ^ c)) + d

where R is the daily solar radiation, I0 is the potential insolation (I0HDAY), and Δ is the daily temperature range (Tmax - Tmin), and a, b, c, and d are constants that depend more or less on location. I added the d value so that R is greater than zero when Δ is zero. For Hubbard Brook I use a = 0.78, c = 1.5, d = 0.05, and b = 0.245 + 0.0035 sin(2π (doy + 80) / 366) to account for some residual seasonality, where doy is day of the year.

In BROOK90 subroutine EQUIVSLP is called once to obtain parameters that depend on latitude, slope, and aspect of the surface being simulated. Subroutine SUNDS uses functions FUNC3 and HAFDAY once a day to calculate daylength, potential radiation on a horizontal surface, and the ratio of potential radiation on the slope to that on a horizontal surface. AVAILEN calculates net radiation minus soil heat flux at the top and the bottom of the canopy, separately for daytime and for nighttime. "

Potential Evaporation

Text copied from Ecoshift on module PET:

" In this section, most equations are given in standard algebraic notation, as extended by Shuttleworth and Wallace (1985). The correspondence between algebraic notation and variable names is: found on http://www.ecoshift.net/brook/b90doc.html

A common approach to estimating evaporation calculates a potential evaporation (PE) from weather variables, and then reduces actual evaporation below PE in response to soil drying. PE quantifies what the evaporation rate would be in the absence of any limitation of liquid water supply to the evaporating surfaces. It is thus an upper limit to the evaporation rate. Actual evaporation falls below the potential rate whenever liquid water supply to the plant leaves or to the soil surface cannot maintain the PE rate. Changing knowledge of limitations on evaporation has produced a variety of definitions of PE, which are related to the methods chosen to calculate it (Shuttleworth, 1991).

Federer et al. (1996) use the term "potential evaporation" (PE) as a generic term to include the general concept and all definitions and methods. Then they distinguish three fundamentally different definitions of PE by subscripts. Reference-surface PE (PEr) is defined as the evaporation that would occur from a specified or reference land surface in given weather conditions with soil water at field capacity (Shuttleworth, 1991); the reference surface is normally defined as a short, complete, green plant cover. Surface-dependent PE (PEs) is defined as the evaporation that would occur from any given land surface in given weather conditions if plant surfaces were externally dry and soil water was at field capacity. Potential interception (PEi) is defined as the evaporation that would occur from any given land surface in given weather conditions if all plant and soil surfaces were externally wetted, as by rain. PEi and PEs depend on surface characteristics, most notably on canopy height and stomatal resistance respectively.

Shuttleworth and Wallace (1985) applied the well-known Penman-Monteith equation separately for the canopy and for the soil surface to give separate estimates of transpiration and soil evaporation. The Shuttleworth-Wallace approach was designed to be applicable to canopies of any leaf area index or "sparseness". Federer and others (1996) developed the Shuttleworth-Wallace method to estimate PEs and PEi for any land surface. To obtain the potential soil evaporation component, they used a soil surface at field capacity for PEs and a saturated soil surface for PEi.

However, BROOK90 separates evaporation into two pathways and five processes controlled by five resistances to vapor transfer:

canopy evaporation
    IRVP    raa + rac         evaporation of intercepted rain
    ISVP    raa + rac         evaporation of intercepted snow
    TRAN    raa + rac + rsc   transpiration
ground evaporation
    SNVP    raa + ras         evaporation from snow on the ground
    SLVP    raa + ras + rss   soil evaporation

where raa is the above canopy resistance, rac and ras are the within canopy resistances from the canopy and from the ground (soil/snow), rsc is the canopy surface resistance (predominantly stomatal resistance), and rss is the resistance to vapor movement within the soil. These are the five resistances in the Shuttleworth-Wallace equations.

In BROOK90 only one canopy process and one ground process can occur at any given time. For IRVP, ISVP, and SNVP, the evaporating surfaces are always assumed saturated when the processes are occurring. In other words, they always evaporate at their own potential rate. The output value PINT is the interception evaporation that would occur if the canopy were continually wetted by snow or rain. IRVP and ISVP are equal to PINT as long as intercepted rain or snow remain on the canopy. Snow interception is treated identically with rain interception, with no allowance for melting the intercepted snow. When there is no interception, BROOK90 calculates and outputs potential transpiration, PTRAN, which is the transpiration that would occur if stomatal opening were not restricted by water supply to the leaves. It then reduces transpiration below PTRAN if low water supply to the leaves causes stomatal closure. For SNVP, BROOK90 avoids or neglects numerous problems of snow energy balance and any interaction with canopy evaporation by using only the vapor gradient and raa + ras, effectively a potential rate (see SNO-SNOVAP).

Because the ground evaporation and the canopy evaporation are calculated simultaneously in the Shuttleworth-Wallace equations, the potential interception or transpiration depend on the value of rss (which is zero for snow or saturated soil). Potential conditions for the canopy processes often do not occur at the same time as potential conditions for the surface processes. When there is no snow on the ground, BROOK90 calculates potential canopy evaporation using the ambient soil resisistance, rss. When there is snow, BROOK90 calculates surface snow evaporation using raa + ras and the vapor gradient, and canopy evaporation using rss = 0. The equations then provide the soil (ground) evaporation (GER) corresponding to PTR and the soil evaporation (GIR) corresponding to PIR. If actual transpiration is less than PTR, soil evaporation is recalculated.

Theoretically the Penman-Monteith and Shuttleworth-Wallace equations are only valid for a short time period over which the weather variables are effectively constant. However, in spite of their non-linearities, the equations also give reasonable values when used with weather variables averaged over daily and even monthly periods (Federer et al. 1996). Although BROOK90 uses only daily weather data, PE estimates could be made for shorter time periods such as hourly by assuming diurnal distributions of weather variables. However, this would be costly of computer time. Tanner and Pelton (1960) suggested that the Penman equation would be more accurate when used on a daily basis if daytime and nighttime were separated. Federer et al. (1996) show that separation of daytime and nighttime also works well with the Shuttleworth-Wallace method.

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.

Potential Evapotranspiration

In BROOK90, PEs (surface-dependent potential evapotranspiration) is not directly calculated, because that would involve changing the value of rss to its potential value (at field capacity), which would then change the values of IRVP, ISVP, and PTRAN. Can users obtain a value of PEs for comparison with other studies? In July 2013, in response to a query by Stefan Plötner, I have figured out how PEs can be obtained. BROOK90 must be run with two parameter changes:

In Fixed Parameters set RSSB to zero. This forces rss = RSSA, which is defined in BROOK90 as the value at "field capacity", and forces SNVP to be at the potential rate. Shuttleworth and Gurney (1990) point out that appropriate values of rss are poorly known (as is the dependence of rss on soil wetness). BROOK90 suggests using their value of 500 s/m for RSSA. In fact, it is perfectly feasible to avoid the "field capacity" concept by DEFINING potential soil evaporation as that which would occur at a fixed rss such as 500 s/m. Note that the potential SLVP is nearly inversely proportional to this rather arbitrary value. In Soil Parameters set THICK(1) = 9999 (max allowed) and NLAYER = 1. This forces a very thick top soil layer that should not dry out. If THICK is too small, SLVP at the potential rate may make the top layer water content go negative, which crashes the model. These parameter changes may increase SLVP a lot, and will reduce PTRAN and ISVP a little, but will not change ISVP and SNVP. An estimate of PEs can then be obtained as: PEs = PTRAN + SLVP + IRVP

• ISVP + SNVP.

"

Snow accumulation and melt

Text copied from Ecoshift on module SNO:

" Simulation of snow accumulation and melt is a complex subject (U.S. Army Corps Eng. 1956, Anderson and Crawford 1964, Colbeck and Ray 1979). Detailed expressions for the energy balance of snow surfaces have been developed (Anderson 1976, Dingman 1994) but they are not generalized to all cover types and are too complex for BROOK90. Leavesley and Striffler (1979) give an energy balance model that includes radiation melt and the effects of canopy cover on it, but ignores convection-condensation melt. More complex algorithms could be developed, but Colbeck et al. (1979) say "the energy exchange processes between snow and a forest cover are not well enough understood to allow detailed modelling of the melt process through use of the energy equation." The energy balance is made complicated because of the heat of fusion in the water-ice phase change, and because the surface vapor pressure for melting snow is fixed. Application of the Shuttleworth-Wallace two layer approach to snow under sparse canopies remains in the future.

BROOK90 therefore falls back on the classic degree-day method for estimating snow energy balance. Anderson (1979, p. 342) states that "Air temperature (ambient) is an adequate index to surface energy exchange in most cases." This is not a perfect solution as Anderson points out three cases in which air temperature fails: 1) warm temperatures with little wind causes overestimates of melt, 2) high dewpoint with high wind causes underestimates, and 3) low temperatures with clear sky and ripe snow causes underestimates. A modified degree-day method that incorporates solar radiation and wind speed could be added but would require development for sparse canopies. Another improvement would separate the day into daytime and nighttime as is done in BROOK90 for evaporation. But BROOK90 currently uses only the mean daily temperature (TA) for the snow energy balance.

The "water equivalent" of snow (SNOW, mm) is the depth of water a snowpack would produce if it were all melted; this is the BROOK90 variable that represents the snowpack. The actual depth of snow, assuming a constant snow density (SNODEN), is used only to calculate the amount of the canopy above the snow in subroutine CANOPY. Variable snow density (mass per unit volume) is not simulated in BROOK90. When the snow is colder than 0°C, it has a "cold content" (CC, MJ/m2), which is the energy needed to warm the snow to 0°C. When the snow is at 0°C, part of SNOW can be liquid water (SNOWLQ, mm). The maximum liquid water that can be retained by the snowpack without draining is a constant fraction (MAXLQF) of SNOW; CC and SNOWLQ are always initialized as zero in BROOK90, so any initial SNOW is considered to be at 0°C with no liquid water.

Groundmelt is snowmelt at the bottom of a snowpack; it occurs because of heat conduction from the soil whenever the soil is unfrozen. A constant groundmelt rate (GRDMLT, mm/d) is an constant parameter in BROOK90 and is applied whenever there is snow on the ground. The possibilities of frozen soil or variable groundmelt are not considered.

Snowmelt (SMLT, mm/d) is the sum of groundmelt and drainage of excess liquid water from the snowpack. Drainage occurs only after the snowpack is both isothermal at 0°C and is holding the maximum possible liquid water; the snowpack is then "ripe". The gains and losses of liquid water by the snowpack, including the refreezing of rain on cold snow, are handled in the somewhat complicated subroutine SNOWPACK.

BROOK90 assumes that the snowpack is always isothermal. In reality, large and variable temperature gradients can exist in thick snowpacks; simulating these is beyond the scope of BROOK90. The snowpack temperature (TSNOW) at the beginning of the day is calculated in MSBSETVARS from the cold content

TSNOW = -CC / (CVICE * SNOW)

where CVICE is the heat capacity of ice (0.00192 MJ m-2 mm-1 K-1) (Leavesley and Striffler, 1979). TSNOW is used both in calculating snow evaporation (SNVP) and snow energy flux (SNOEN). "

Interception and transpiration

Text copied from Ecoshift on module EVP:

" BROOK90 routines in this module relate to interception and actual transpiration. Subroutines INTER and INTER24, which handle the interception of rain or snow by the plant canopy are equivalent routines. INTER24 is used when parameter p_NPINT is 1 and precipitation is input once a day; it assumes that the daily precipitation all occurs within DURATN hours in the middle of the day. INTER is used when p_NPINT > 1 and precipitation is input more than once a day; it assumes that precipitation rate is constant through the precipitation time step. INTER and INTER24 are used both for rain and snow, with different calling parameters and variables. PLNTRES calculates parameters related to rhizosphere, root, and xylem resistance; it is called once at the beginning of each day. These parameters affect only soil water supply rate, not potential transpiration. TBYLAYER calculates the daily transpiration from each soil layer from the potential transpiration (PTRAN) and the total soil water potential in each layer (PSITI). "