BROOK90 - PET - POTENTIAL EVAPORATION

Potential ET - PM - ESAT - SWPE - SWGE - SWGRA - CANOPY - ROUGH - SRSC - FRSS - WEATHER - WNDADJ

July 10, 2015

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:

and for the parameters and constants is:

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 (PE_r) 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 (PE_s) 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 (PE_i) 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. PE_i and PE_s 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 PE_s and PE_i for any land surface. To obtain the potential soil evaporation component, they used a soil surface at field capacity for PE_s and a saturated soil surface for PE_i.

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 r_aa is the above canopy resistance, r_ac and r_as are the within canopy resistances from the canopy and from the ground (soil/snow), r_sc is the canopy surface resistance (predominantly stomatal resistance), and r_ss 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 r_aa + r_as, 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 r_ss (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, r_ss. When there is snow, BROOK90 calculates surface snow evaporation using r_aa + r_as and the vapor gradient, and canopy evaporation using r_ss = 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, PE_s (surface-dependent potential evapotranspiration) is not directly calculated, because that would involve changing the value of r_ss to its potential value (at field capacity), which would then change the values of IRVP, ISVP, and PTRAN. Can users obtain a value of PE_s for comparison with other studies? In July 2013, in response to a query by Stefan Plötner, I have figured out how PE_s can be obtained. BROOK90 must be run with two parameter changes:

1. In Fixed Parameters set RSSB to zero. This forces r_ss = 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 r_ss are poorly known (as is the dependence of r_ss 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 r_ss such as 500 s/m. Note that the potential SLVP is nearly inversely proportional to this rather arbitrary value.
2. 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.

Function PM - Penman-Monteith equation

The Penman-Monteith equation is

where E is the evaporation rate in volume of water per unit land area per unit time, L_v is the latent heat of vaporization for water, ρ_w is the density of water, Δ is the rate of change of vapor pressure with temperature, R_n is the net radiation above the surface, S is the subsurface heat flux, c_p is the heat capacity of air, ρr is the density of air, D_a is the vapor pressure deficit in the air, γ is the psychrometer constant, r_c is the "canopy resistance", and r_a is the aerodynamic resistance between the canopy and a reference height z_a at which D_a is measured. The vapor pressure deficit, D_a, is e_a* - e_a. The equation assumes that the vapor pressure at the effective evaporating surface, e₀, is the saturated vapor pressure at the surface temperature. Then r_c and r_a are the two "resistances" through which water vapor passes as it moves down the vapor pressure gradient from e₀ to e_a. The canopy resistance, r_c, 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, r_a, is a measure of the turbulent transfer capability of the atmosphere between the effective source height and z_a. 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).

Subroutine ESAT - Saturated vapor pressure

The subroutine ESAT calculates the saturated vapor pressure (kPa), e*, from temperature (°C), T , and calculates Δ. The equations of Murray (1967) are used. If T >= 0°C then

e* = 0.61078 exp[ 17.26939 T / ( T + 237.3 ) ]

Δ = 4098 e* / ( T + 237.3 )²

or if T < 0°C

e* = 0.61078 exp[ 21.87456 T / ( T + 265.5 ) ]

Δ = 5808 e* / ( T + 265.5 )².

Subroutine SWPE - Shuttleworth-Wallace Equations

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), r_c and r_a, SW define five: r_sc, r_aa, r_ac, r_as, and r_ss. Resistances r_sc and r_ac are in the transpiration pathway while r_ss and r_as are in the soil evaporation pathway and r_aa is common to both. The canopy surface resistance, r_sc, is the resistance to movement of water vapor out of the leaves. The resistance r_ac 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 r_aa, which corresponds to r_a 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 r_ss. The resistance to vapor movement from the soil surface to the source height is r_as. The resistances r_ac, r_as, and r_aa are assumed also to apply to sensible heat transfer.

Shuttleworth and Wallace (1985) start with

(2) L_v ρ_w E = L_v ρ_w E_c + L_v ρ_w E_s

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

where D₀ is the vapor pressure deficit at the effective source height, A is R_n - S or the available energy above the canopy, and A_s is R_ns - 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) D₀ = D_a + r_aa [Δ A - (Δ + γ) L_v ρ_w E ] / c_pρ

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

(6) L_v ρ_w E = C_c M_c + C_s M_s

where

(9) C_c = 1 / { 1 + R_c R_a / [ R_s ( R_c + R_a ) ] }

(10) C_s = 1 / { 1 + R_s R_a / [ R_c ( R_s + R_a ) ] }

with

(11) R_a = (Δ + g) r_aa

(12) R_s = (Δ + γ) r_as + γ r_ss

(13) R_c = (Δ + γ) r_ac + γ r_sc

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 D₀ 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 E_c (PRATE) and E_s (ERATE) from SWPE are in units of mm/d whereas L_vρ_w E in (1) is output as W m^-2 from function PM. The conversion is ETOM * WTOMJ .

Subroutine SWGE - Shuttleworth-Wallace with known transpiration

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, E_c (ARATE), in subroutine SWGE to get a new value of soil evaporation, E_s (ERATE). With E_c known, substituting(5) into (4) and (4) into (2) and solving for L_v ρ_w E gives

then

(15) L_v ρ_w E_s = L_v ρ_w E - L_v ρ_w E_c

Subroutine SWGRA - Shuttleworth - Wallace - Gurney Aerodynamic Resistances

The three SW aerodynamic resistances, r_aa, r_as, and r_ac 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) u* = k u_a / ln [ ( z_a - d ) / z₀ ]

where u_a is the wind speed at the reference height, z_a, z₀ 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 z₀ 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, r_as, is obtained by integrating the exponential eddy diffusivity from 0 to D_c to give

(17) r_as = [ h exp(n) / n K_h] [exp( -n z_0g / h ) - exp( -n D_c / h ) ]

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

(18) r_aa = ln [ ( z_a - d ) / ( h - d ) ] / k u* + ( h / n K_h ) { -1 + exp[ n (h - D_c) / h ] }

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) 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, r_ac, is assumed by Shuttleworth and Gurney (1990) to be a sum of the individual leaf laminar boundary layer conductances, so

(20) r_ac = ( n / a_b ) ( w / u_h )^1/2 / { ρ t_p L_p [ 1 - exp( -n / 2 ) ] }

where n is the extinction coefficient for eddy diffusivity, a_b is a constant = 0.01 m/s^0.5 (Campbell 1977), ρt_p is the ratio of projected leaf area to total leaf surface area, L_p is the projected leaf area index, w is the representative leaf width in m, and u_h 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_ac equation strictly applies only to flat leaves and needles, not to cylindrical needles; but r_ac is small when w is small so this inaccuracy is negligible. The canopy-top wind speed, u_h, is

(21) u_h = (u* / k) ln [ ( h - d ) / z₀ ].

In the Shuttleworth and Gurney equation, r_ac goes to infinity as L_p goes to zero, but this neglects the contribution of stems and branches to interception loss, especially when L_p = 0. Although the r_ac equation is specifically for leaves with width w, subroutine SWGRA uses (ρt_p L_p + π S_p) for r t_p L_p in equation (20), where S_p is the projected stem area index, defined analagous to L_p 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 r_ac as a unit of leaf surface. To avoid zero divides when S_p = 0, L_p is prevented from being less than 0.00001.

The Shuttleworth - Gurney separation of r_ac and r_sc 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).

Subroutine CANOPY - Canopy structure parameters

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

The height of the canopy above any snowpack, h (HEIGHT), is

HEIGHT = RELHIT * MAXHT - SNODEP

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 (SNODEP, 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.

Actual projected leaf area index, L_p (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 S_p (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.

Subroutine ROUGH

ROUGH obtains the roughness parameter, z₀ , and the zero-plane displacement, d, based on canopy height, h, the projected leaf area index, L_p, and the projected stem area index, S_p. 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 L_p is greater or less than some arbitrary value L_pc, which they take as 4. Following Federer et al. (1996), BROOK90 defines a closed canopy as having L_p + S_p greater than L_pc + S_pc, where S_pc is taken as c_s h, as described in the previous section. S_pc is not reduced by DENSEF. RATIO is ( L_p + S_p ) / ( L_pc + S_pc ) (this RATIO differs from RATIO in subroutine CANOPY). When RATIO is greater than or equal to 1, the canopy is "closed" and z₀ and d are the values for a closed canopy.

The roughness parameter for closed canopies, z_0c (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 z_0c / h as follows: c_zs is the ratio for smooth surfaces with h less than h_s; c_zr is the ratio for rough surfaces with h greater than h_r; and for h between h_s and h_r, z_0c is interpolated linearly between z₀ = c_zs h_s at h = h_s and z₀ = c_zr h_r at h = h_r. The effect of stems and branches is assumed to be included in z_0c.

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

(22) d_c = h - z_0c / 0.3

(Monteith 1973; Shuttleworth and Gurney 1990).

When h (HEIGHT) is small, it is possible for z_0c to be less than z_0g (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 z_0g to the value of z_0c, which works nicely because z_0c is always added to d_c.

The values of d and z₀ 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, c_d (CDRAG), is calculated from z_0c and d_c as

(23) c_d = { -1 + exp [ 0.909 - 3.03 ( z_0c / h ) ] } 4 / (L_pc + S_pc) .

This assumes that a unit of S_pc produces the same drag as a unit of L_pc. Then

(24) 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 c_d of 0.07, which is not appropriate when z_0c / h is not 0.13.

For z₀, 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 c_d L_p. However, that procedure leads to a step change at the transition. The step change is avoided by choosing the minimum of the two z₀ values, as

(25) z₀ = min { 0.3 ( h - d ), z_0g + 0.3 h [ c_d ( L_p+ S_p) ] ^0.5 }

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

The reference height at which weather variables are known, z_a (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.

Subroutine SRSC - Canopy surface resistance

This routine obtains the canopy surface resistance, r_sc, 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 m³_H2O m^-2_leaf s^-1, and leaf area index, L, in units of m^-2_leaf m^-2_ground to give a canopy conductance in m³_H2O m^-2_ground s^-1. But the leaf area index can be given for projected leaf area, L_p, or for total leaf area, L_t. Similiarly, the conductance can be for projected leaf area, g_p, or for total leaf area, g_t. But the total evaporation must be the same, so g_p L_p = g_t L_t is required. In the literature, g_t L_t is commonly used for needle-leaved plants and g_p L_p for broadleaved plants. Following Körner et al. (1979) BROOK90 uses g_p L_p for all plants. For needles, where g_t, or its reciprocal r_t, is known, then g_p = ρ t_p g_t = ρt_p / r_t, where ρt_p = L_t / L_p. Note that ρt_p can vary from 2 for flat needles to π for cylindrical needles. For broad leaves, the leaf resistance is usually given separately as r_b for the abaxial (lower) surface, and r_d 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 g_p L_p = g_b L_p + g_d L_p = L_p ( 1 / r_b + 1 / r_d ). For amphistomatous leaves (stomates on both sides) the simplification g_b = g_d is often assumed, giving g_p L_p = 2 g_b. The factor of 2 can easily be confused with ρt_p = 2, but they are not the same. For hypostomatous leaves (stomates only on the abaxial side), g_d = 0 is often assumed, giving g_p L_p = g_d L_p. BROOK90 uses only g_p, which hereafter is called g_l, and is independent of assumptions about how g_p 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) g_l = g_lmin + f_T f_D f_R f_W f_C ( g_lmax - g_lmin )

where g_l is the leaf conductance, g_lmin is its minimum value (closed stomates), g_lmax is its maximum value, and f_T, f_D, f_R, f_W, and f_C 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 f_C = 1 always.

BROOK90 also always uses f_w = 1. The resulting r_sc is therefore appropriate to the definition of potential transpiration as occurring when soil water is at field capacity. Reduction of f_w 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 f_w = 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 f_T = 1 and inverted half parabolas to reduce f_T to 0 at each end of the range, so

T_L, T₁, T₂, and T_H are all parameters. If T₁ = T₂, there is no optimum range. If T_L = T₁ and T₂ = T_H, f_T is a square wave function. The use of mean daily temperature, T_a, 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, D_a (VPD), is usually used. The error induced by using D_a 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) f_D = c_D / ( c_D + D_a )

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

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

(29) f_R = (R_m + R₀) R_L / [ R_m (R + R₀) ]

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

(30) R₀ = R_m / [ ( R_m / R_.5 ) - 2 ]

For most plants, the value of R_.5 is relatively low, around 50 to 100 W/m². 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, R_L, 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, g_sc, which is the reciprocal of r_sc, is the integral of g_l over each increment of total leaf area in the canopy , dL' , so

With (26), where f_T, f_D, and g_lmax are all assumed constant through the canopy, (31) becomes

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

Integration of f_R 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, R_L, at any level in the canopy depends on the projected leaf area, L', and projected stem area, S', above that level, as follows:

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

where C_R is the extinction coefficient and R is the solar radiation at the top of the canopy. The first C_R 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 S_p is distributed proportionally to L' (and the other half is below the leaves in a "stem space"), then L' + S' = L' + (L'/L_p)(S_p/2) = f_s L' where f_s = (L_p + S_p/2) / L_p. The integral in (32) using (29) and (33) is then

which corresponds to Shuttleworth and Gurney (1990) and Saugier and Katerji (1991) when S_p = 0. The combination of (32) and (34) provides the value of r_sc (RSC)

Function FRSS

FRSS returns the Shuttleworth-Wallace soil surface resistance, r_ss (RSS), which must increase with drying of the surface. BROOK90 currently uses only the top soil layer to calculate RSS, no matter what its thickness. The model assumes that the ratio of matric potential (PSIM) to matric potential at field capacity (PSIF) in the surface layer is the controlling variable. Function FRSS is

(36) RSS = RSSA * (PSIM / PSIF) ^ RSSB

where the parameter RSSA is RSS at field capacity, and RSSB is an exponent. The value of RSSA is poorly known; BROOK90 suggests using the Shuttleworth-Wallace value of 500 s/m; soil evaporation (SLVP) will be inversely proportional to this value. Using RSSB = 1 makes SLVP inversely proportional to PSIM in the top layer. Certainly more investigation is required before changing these values for any soil surfaces. A valid relationship is yet to be developed even for bare soil (van de Griend and Owe 1994).

Currently in BROOK90 the presence of intercepted rain (INTR) does not affect RSS. That is, the soil surface is not assumed to be saturated during and after storms. An alternative would be to make RSS = 0 when INTR > 0.

When there is SNOW, RSS is set to zero, which affects calculated transpiration and interception rates, but there is no soil evaporation and snow evaporation is calculated using subroutine SNOVAP. This leaves something to be desired, but the whole situation with melting and evaporating snow as well as transpiration is a complicated mess.

When the parameter RSSA is set to 0 there is no soil evaporation (SLVP).

Subroutine WEATHER

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 temperatire (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 z_a.

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 z_a. However, UA can be estimated from wind speed at the weather station (UW) because wind speed extrapolates to zero at height z₀ + 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.

Function WNDADJ

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 z_0w (Z0W) and a fetch of F (FETCH). Brutsaert (1982) gives the height of the internal boundary layer, z_b, 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 z_b,

where z_w (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 u_a = u_w.

COMPASS BROOK