SPA handbook v 1.2

The Soil-Plant-Atmosphere model

Manual, version 1.2, April 2005

Dr. Mathew Williams

School of Geosciences

University of Edinburgh

Darwin Building

Kings Buildings

Edinburgh EH9 3JU, UK

changes to version 1.1.

New note on the iota parameter in “Running the model”.

Index

The Soil-Plant-Atmosphere model

Manual, version 3.0, 4 December 2002

Introduction

Ecological and Physical Fundamentals

Radiative Transfer

Boundary layer conductance

Plant and leaf ecophysiology

Soil energy and water balance

Model Subroutines (SR) and Functions

All_declarations.f90

Scale declarations.f90

Canopy.f90

Leaf.f90

Light.f90

Soil functions.f90

Soil and air.f90

RUNGE_KUTTA.f90

ZBRENT.f90

IO.f90

Main.f90

Running the model

Model Applications

History

References

Introduction

The Soil-Plant-Atmosphere model (SPA, Williams et. al 1996) is a process-based model that simulates ecosystem photosynthesis and water balance at fine temporal and spatial scales (30 minute time-step, ten canopy and twenty soil layers). The scale of parametrization (leaf-level) and prediction (canopy level) have been designed to allow the model to diagnose eddy flux data, and to provide a tool for scaling up leaf level processes to canopy and landscape scales (Williams et al. 2001c).

The model is written in FORTRAN 90. The code is divided among several files, each associated with a particular component of the overall model. This means that the core files can be shared among users, who should only have to alter the input/output file and the main program file in order to customise the model for their own uses.

In this document we will describe the underlying ecological and physical principles that govern the construction of SPA. We will also work through and describe the subroutines and functions that make up the code. We will describe the input and output files – how to create drivers, alter parameters, and interpret model products.

Ecological and Physical Fundamentals

Radiative Transfer

The model employs a detailed radiative transfer scheme that determines the time-varying transmittance, reflectance and absorption of long wave, near infra-red and direct and diffuse photosynthetically active radiation (PAR) by canopy layers and the soil surface. Absorption of PAR in each canopy layer is determined by Beer-Lambert’s Law, and is partitioned between sunlit and shaded foliage fractions.

The model requires, as inputs, measured PAR above the canopy at each time step. Short wave radiation can be estimated from PAR or vice versa. If the data are available, then the model can be driven with measured proportions of direct and diffuse PAR for each time step. If these data are not available, then the ratio of diffuse to direct radiation is estimated using an empirical formula (Erbs et al. 1982). Solar geometry (i.e. solar declination) is determined from the day, the hour and the latitude of the simulation.

Canopy interaction with radiation is determined by the area and distribution of foliage. At its simplest the canopy is described by the vertical patterns of leaf area density (the amount of leaf area in each canopy layer, the sum of which gives the leaf area index, or LAI). The model assumes no foliage clumping, a spherical leaf angle distribution, and no interaction of radiation with non-photosynthetic tissues. There is a clumping parameter (‘clump’) which can be adjusted if relevant data are available, and the leaf angle parameter (‘spherical’, G in the equation below) can also be adjusted manually.

We assume a spherical leaf angle distribution (Russell et al. 1989), so the proportion (P) of radiation incident on the top a layer that passes through without striking a leaf is given by

P = exp(-GLi/sin)

where G is 0.5 for a spherical leaf angle distribution, Li is the leaf area (m2 m-2) in canopy layer i, and  is the elevation of the radiation source (which varies with latitude and time of day for beam radiation (Jones 1992); for diffuse radiation it is held constant at 30).

Radiation not intercepted in a layer is incident on the layer below (for downward radiation) or the layer above (for upward radiation). Radiation that strikes leaves in a layer is either absorbed, reflected or transmitted. Reflectivity and transmissivity of leaves and soils are estimated from empirical data (Baldocchi & Harley 1995). After the first downward pass, there is an upward pass through the canopy to determine the fate of the reflected radiation. Any beam radiation that is reflected or transmitted is converted to diffuse radiation. This whole process is repeated two or three times until all radiation has been absorbed by leaves or soil, or reflected back to the sky. Both PAR and NIR regimes are calculated in a similar fashion, but separately, because leaf optical properties in relation to these wavebands differ markedly.

In 1998 the model was modified to introduce a revised radiation subroutine that calculates sun-lit and shaded fractions of the foliage in each canopy layer (Norman 1981) for the PAR routines only. Sun-lit foliage fractions are lit by both incident direct and diffuse PAR, while the shaded fractions receive only diffuse PAR. For each canopy layer the model calculates absorbed radiation separately for the sun-lit and shaded fractions.

Sunlit leaf area (Ls) is determined by

Ls = (1-e(-k*L))/k

Where k = G/(sin). The sunlit fraction of each layer receives the full amount of the direct beam radiation incident at the top of the canopy, plus the attenuated diffuse radiation that penetrates to that layer.

Radiative transfer of NIR follows the same principles as for PAR, but there are different reflectance/transmittance constants, and no separation of sunlit and shaded leaf fractions.

For long wave energy balance, down-welling long wave radiation is either provided directly from data, or we estimate it under clear skies (Ld) from an empirical relation (Jones 1992)

Ld (Ta - 20)4

The upward emittance of long wave radiation of the canopy layers (l, W m-2) is estimated by

lu Ta4 (1- exp(-kLi))

where Ta is air temperature (K) of the canopy layer, which approximates leaf temperature, k is the extinction coefficient (0.5, as for diffuse radiation), and Li is leaf area in the canopy layer (m2 m2). The final term scales emittance according to reabsorptance within the same layer(Amthor et al. 1994). The downward emittance of longwave radiation is set equal to the upward emittance. The fate of longwave radiation is determined in a similar manner to short-wave radiation. The sum of PAR (derived from PPFD), NIR and longwave radiation absorbed (la) determines the net isothermal radiation (ni; W m-2) for each canopy layer. Long wave emissions from foliage are estimated by assuming leaves are at ambient temperatures. Thus the radiation model is not currently coupled directly to the leaf energy balance model in the leaf level process routines. Such a coupling would increase the model accuracy, but at a large computational cost, and likely would only make a small change to current predictions.

Boundary layer conductance

The multilayer approach requires that we estimate the variation of boundary layer characteristics of leaves within the canopy. Windspeeds, and thus boundary layer conductances, decline within the canopy (Roberts et al. 1990). We use an exponential wind relationship (Cionco 1985) to approximate the windspeed at different heights within the canopy;

u(h) = uz exp( (h/hz - 1))

where uZis measured windspeed above the canopy at height hZ (m), u(h) is windspeed at height h (m) within the canopy, and  is the canopy flow index. We use a value of 4.0 for  from Cionco (1978) for maple-fir forest.

Within each canopy layer, we approximate the thickness of the boundary layer of a leaf (m) using a relationship from Nobel (1983);

where llis the characteristic dimension of the leaf (m). For a default value, we use 0.08m, a representative measure of the size of Quercus/Acer leaves.

One-sided leaf boundary layer conductance to water vapour (gb; m s-1) varies with the diffusion coefficient of water vapour in air (Dwv; m2 s-1) and the thickness of the boundary layer (),

gb = Dwv/

Boundary layer conductance to heat (gH; m s-1) is calculated in a similar manner, using thermal diffusivity (DH; m2 s-1 ), and the factor 2 accounting for transfer from both leaf surfaces,

gH = 2DH/

Plant and leaf ecophysiology

The SPA model employs some well tested theoretical representations of eco-physiological processes, such as the Farquhar model of leaf-level photosynthesis (Farquhar & von Caemmerer 1982), and the Penman-Monteith equation to determine leaf-level transpiration. These two processes are linked by a novel model of stomatal conductance that optimizes daily carbon (C) gain per unit leaf nitrogen (N), within the limitations of canopy water storage and soil-to-canopy water transport. The maximum flux rate of water through vegetation is determined by the difference between soil water potential and the minimum sustainable leaf water potential, and by the hydraulic resistance of the soil-root-leaf pathway. Stomata adjust to equalise evaporative losses with the maximum hydraulic supply, minimising the risk of cavitation.

The principle behind stomatal dynamics is that C assimilation is maximised within the limitations of the hydraulic system, so stomatal resistance is adjusted to balance atmospheric demand for water with rates of water uptake and supply from soils. Atmospheric demand is governed by the vapour pressure difference between leaf internal air spaces and the atmosphere, and vapour phase exchange (E) is determined using the Penman-Monteith equation. Water loss (= E) is linked to changes in leaf water potential (l), according to the water potential gradient between leaf and soil, liquid-phase hydraulic resistances (both in the rhizosphere, Rs, plant stems, Rp, and roots, Rr) and the capacitance (C) of the pathway that links soil to leaf. Stomatal resistance is varied to maintain evapotranspiration (E) at the rate that keeps Ψl from falling below a critical threshold value (Ψlmin), below which potentially dangerous cavitation of the hydraulic system may occur. Thus, once Ψl = Ψlmin, E is set so that dΨl /dt= 0, where

The gravitational component of leaf water potential is determined by the density of water (w), acceleration due to gravity (g), and the height above the reference plane (h), set to the soil surface.

Plant hydraulics

We assume that the plant resistance per unit leaf area to flow of water along the xylem to a layer in the canopy (Rp; MPa s m2 mmol-1) increases with the height (h; m) of each layer above the ground;

where Gp is the canopy hydraulic conductivity (mmol s-1 m-1 MPa-1). The layers higher in the canopy are subject to the greatest resistance to xylem water supply (Hellkvist et al. 1974). In some plants there is evidence that resistance is not strongly related to path length (it may be that the resistance is concentrated at nodes in the conductive network, for example). In such a case, it may be better to specify a constant plant resistance for each canopy layer. There is a switch in the model (‘CONDUCTIVITY’) to select either option. More details are given in the operating instructions.

The following equations can be used to describe the relationship between steady-state flow and water potential between the soil and the roots, between the roots and the leaves, and for the whole plant system;

where Jp is water flux (mmol m-2 s-1), l is leaf water potential (MPa), d is gravitational head of pressure (MPa m-1), and the denominators contain belowground and aboveground resistances to water transport (MPa s m2 mmol-1).

The water in the xylem can rupture under the extreme tensions that occur naturally - there is often a threshold water potential for such cavitation (Jones 1992). The onset of xylem cavitation can lead to a rapid and catastrophic decline in stem hydraulic conductivity by inducing further cavitations (Tyree & Sperry 1989). Jones and Sutherland (1991) have shown that the maintenance of a maximally efficient conducting system requires that stomata close as evaporative demand rises to prevent shoot water potentials falling below the threshold value (lmin).

The relationship between Jp and the water potential drop is not unique; initially water is drawn from stores within plant tissues, so that liquid flow lags behind the evaporative demand (Landsberg et al. 1976; Schulze et al. 1985). This hysteresis can be modelled by incorporating capacitors into the circuit analogue (Jones 1978). The capacitance (C, mmol MPa-1 m2) of any part of the system is defined as the ratio of the change in tissue water content (W, mmol m-2) to the change in water potential ();

As shown above, the flow through the plant to the leaves (Jp; mmol m-2 s-1) can be written;

Jp = (s - l - hd)/(Ra +Rb)

The rate of change of layer water content (dWl/dt) is given by the difference between the flow of water into the layer and that lost by evaporation;

dWl/dt = Jp - E

Thus,

dWl/dt = (s - l - hd)/(Ra +Rb) - E

Assuming constant capacitance, the first order differential equation describing leaf water potential is;

Calibration of the belowground resistances is described below (Root distribution and water uptake).

Leaf level processes

For each canopy layer, once every half-hour an iterative procedure is used to determine the maximum stomatal conductance (gs) and the assimilation rate associated with this conductance. We use a bisection routine to calculate the stomatal conductance very rapidly, but in essence the procedure is hypothesized to work in the following manner:

Starting from a very low gs:

  1. Increment gs.
  2. Determine leaf temperature (Tl, C) resulting from the leaf energy balance at this gs using a steady-state approximation (Jones 1992).

where Ta (C) is air temperature in the canopy layer, ni is net isothermal radiation (KW m-2),r are leaf resistances (m s-1; subscripts: rHR parallel heat and radiative transfer (Jones 1992); rW water; ra boundary layer; rl leaf),  is the psychrometer constant (KPa K-1), S is the rate of change of saturation vapour pressure with temperature (KPa K-1, assumed constant between Tl and Ta), e is vapour pressure deficit (Pa),  is air density (g m-3), and cp is the specific heat capacity of air (KJ g-1 K-1).

  1. Determine biochemical parameters for Farquhar & Von Caemmerer (Farquhar & von Caemmerer 1982) equations. These photosynthetic parameters are dependent on foliar N concentrations (Field & Mooney 1986) and leaf temperature. Photosynthesis is limited by the minimum of the ribulose bisphosphate carboxylation rate and the rate of ribulose bisphosphate (RuBP) regeneration. We determine the maximum rates of these two processes.

Maximum Carboxylation Capacity (Vcmax)

Vcmax = Ntcc where N is leaf nitrogen content (g m-2), and tc is a temperature coefficient (described below), and c is the catalytic rate coefficient, determined from A/Ci curve data.

Maximum electron transport rate (Jmax)

Jmax = Ntjj

where tj is a temperature coefficient.

The temperature dependencies of Jmax (tj) and Vcmax (tc) are described using response curves (Rastetter et al. 1991) fitted to the polynomial relationships of McMurtrie et al. (1992), but with temperature optima set to 30C to reflect those of temperate deciduous vegetation.

Biochemical Constants

The values and temperature dependencies of Kc and Ko, the Michaelis-Menten constants for enzyme catalytic activity for CO2 and O2 respectively are taken from Kirschbaum and Farquhar (1984) and McMurtrie et al. (1992). *, the CO2 compensation point in the absence of non-photorespiratory respiration, is derived from measurements of Vcmax, amongst other parameters. We approximatewoody plant * by 36.5 mol mol-1 at 25C (Epron et al. 1995), using an Arrhenius relationship for temperature sensitivity (McMurtrie et al. 1992).

  1. Determine the mesophyll CO2 concentration (Cc) that satisfies both diffusion and metabolic uptake;

a) Diffusion model;

A = gt (Ca - Cc)(23)

where Ca is ambient air CO2 concentration, A is assimilation rate (mol m-2 s-1) and gt is total CO2 conductance (mol m-2 s-1), determined from stomatal (gs), leaf boundary layer (gb) and internal mesophyll conductances (gi), all m s1, to CO2 ;

(24)

where P is atmospheric pressure (Pa), R is the gas constant (Pa m3 mol-1 K-1) and Tl is leaf temperature (K).

The existence of a limitation to CO2 assimilation by internal resistances to CO2 transfer has been shown in several woody plant species. However, for simplicity, we assume that gi is very large (=1.0) and thus insignificant. To assume otherwise (as we did in early publications) means that calibrating Vcmax and Jmax from A/Ci curves becomes far more complex. This simplification of gi is not really satisfactory, and the issue of finding a way to deal with a realistic mesophyll conductance is important.

b) Metabolic model;

We determine the Rubisco limited carboxylation rate (Wc,mol m2 s-1) as (Farquhar et al. 1982)

(25)

where Cc is CO2 concentration at the catalytic site of Rubisco, Oc is the oxygen concentration at the same site (210 mmol mol-1), and the other parameters are described in (3), above.

The RuBP regeneration-limited carboxylation rate is described by

(26)

where J is the potential rate of electron transport, which is a function of Jmax and absorbed photosynthetic photon flux density (Q;mol m2 s-1), satisfying the relationship (McMurtrie et al. 1992),

J2 - (jQ + Jmax)J +jQJjmax = 0(27)

in which j and  are the initial slope and the curvature of the quantum response of the electron transport slope (Farquhar & Wong 1984).

We determine the actual rate of carboxylation by

Vc = min (Wc, Wj)(28)

and the photosynthetic rate (A, mol m2 s-1) by

A = Vc(1 - */Cc) - Rd(29)

where Rd is respiration rate (mol m2 s-1). Rd is parametrised from leaf level data, and is a function of leaf temperature and foliar N.

  1. Calculate evapotranspiration rate (mmol m-2 s-1) at this gs with the Penman-Monteith equation,

(31)

where gb is leaf boundary layer conductance (m s-1),  is the latent heat of vaporisation (KJ g-1), cw is atmospheric water deficit (g m-3), c converts from g m-2 s-1 to mmol m-2 s-1,  is the ratio of latent heat content increase to sensible heat content increase (unitless),and n is net radiation (KWm-2), approximated by (Jones 1992),

nni -4Ta3(Tl - Ta)(32)

where  is emissivity,  is the Stefan-Boltzmann constant and the temperatures are in K.

  1. Calculate the change in l after one time step of evapotranspiration at the specified gs;

(33)

  1. Return to step 1 (for a further increment of gs), unless either;

a)previous gs increment failed to raise assimilation appreciably. After an initial rapid increase, the response of A to gs becomes asymptotic, and the carbon gain per unit water loss declines. We argue that plants use their water store conservatively; it is more efficient to limit gs when water deficits are low (morning), so that stored water can be utilised later to buffer the impacts of high afternoon water deficits (Cowan 1977). For this reason, optimal gs is fixed when an increment in gs (~1 mmol m-2 s-1) fails to increase A by more than 1- (the threshold parameter iota was initially selected to give maximum conductances of 300 mmol m-2 s1 (Weber & Gates 1990).