MAESPA: a model to study interactions between water limitation, environmental drivers and vegetation function at tree and stand levels, with an example application to [CO2] × drought interactions

Process-based models (PBMs) of vegetation function can be used to interpret and integrate experimental results. Water limitation to plant carbon uptake is a highly uncertain process in the context of environmental change, and many experiments have been carried out that study drought limitations to vegetation function at spatial scales from seedlings to entire canopies. What is lacking in the synthesis of these experiments is a quantitative tool incorporating a detailed mechanistic representation of the water balance that can be used to integrate and analyse experimental results at scales of both the whole-plant and the forest canopy. To fill this gap, we developed an individual treebased model (MAESPA), largely based on combining the well-known MAESTRA and SPA ecosystem models. The model includes a hydraulically-based model of stomatal conductance, root water uptake routines, drainage, infiltration, runoff and canopy interception, as well as detailed radiation interception and leaf physiology routines from the MAESTRA model. The model can be applied both to single plants of arbitrary size and shape, as well as stands of trees. The utility of this model is demonstrated by studying the interaction between elevated [CO 2] (eCa) and drought. Based on theory, this interaction is generally expected to be positive, so that plants growing in e Ca should be less susceptible to drought. Experimental results, however, are varied. We apply the model to a previously published experiment on droughted cherry, and show that changes in plant parameters due to long-term growth at e Ca (acclimation) may strongly affect the outcome ofCa× drought experiments. We discuss potential applications of MAESPA and some of the key uncertainties in process representation.


Introduction
The response of plant carbon uptake and water use to environmental change is complex because there are many interactions and feedbacks that modify the response to single environmental drivers.This complexity is highlighted by the remarkable diversity of experimental outcomes, for example in the response of plant water use and carbon uptake to elevated atmospheric [CO 2 ] (Ainsworth and Long, 2005;Ainsworth and Rogers, 2007;Norby et al., 1999), warming (Rustad et al., 2001), and soil water deficit (Manzoni et al., 2011;Wu et al., 2011).Soil water availability frequently limits plant productivity, but because of interactions with plant properties, and with other environmental drivers, it remains difficult to predict the effect on vegetation water use and carbon uptake (Hanson et al., 2004).Because soil drought is already becoming more frequent as a result of climate change (Huntington, 2006), it is crucial that its effect on vegetation function is more readily quantifiable.
Process-based models (PBMs) can be used as a research tool to clarify interactions among environmental drivers, plant and canopy structure, leaf physiology and soil water availability and their combined effects on water use and carbon uptake (Luo et al., 2008;Norby and Luo, 2004;Williams et al., 2001b).Because PBMs summarize the state-of-art theory of plant functioning in a coherent quantitative framework, they provide a way forward for testing our understanding of how plants respond to environmental change.In this way, they might be used to improve on empirical metaanalyses of experiments.These meta-analyses have been crucial in determining overall responses of vegetation to manipulation of the environment, but the variability among experiments remains incompletely understood (Poorter and Navas, 2003).Meta-analyses usually focus on the effect of a single variable on vegetation function, although it has been recognized that interactions, feedbacks and acclimation are important in determining overall experimental outcomes (Norby and Luo, 2004;Norby et al., 1999;Wullschleger et al., 2002).
Typical experiments or long-term observations are conducted at spatial scales from whole plants (potted seedlings, whole-tree chambers) to entire canopies (eddy-covariance sites, free-air CO 2 enrichment; FACE).A useful PBM to analyse and integrate experiments should therefore operate at both whole-plant and canopy scales.Another requirement is that the PBM incorporate detailed water balance routines and hydraulics of the soil-plant pathway.
Currently, no PBM exists that meets both requirements, and is sufficiently detailed to allow a connection with typical measurements of plant physiological function and plant canopy structure.Here, we introduce a new model, MAESPA, based on a combination of the well-known MAESTRA and SPA models.
The MAESTRA model is a tree array model that uses detailed radiative transfer calculations and leaf physiology to calculate radiation absorption, photosynthesis and transpiration of individual trees, growing singly or in a population (Medlyn, 2004;Wang and Jarvis, 1990).The MAESTRA model has been applied in a wide range of contexts from horticulture (e.g.Bauerle and Bowden, 2011) to regional climate change (e.g.Luxmoore et al., 2000), resulting in well over fifty publications (Medlyn, 2004; updated bibliography available at www.bio.mq.edu.au/maestra/bibliog.html).A major limitation of this model, however, has been the lack of a dynamic water balance.Medlyn et al. (2005) added the capacity to input soil moisture and thereby estimate the effect of a given soil moisture stress on photosynthesis and transpiration (e.g. as used by Reynolds et al., 2009), but the model does not calculate the soil moisture dynamically.This lack severely limits model applications in dry conditions because feedbacks to plant performance via soil moisture cannot be simulated.For example, although the model has often been used to estimate effects of elevated atmospheric [CO 2 ] (C a ) on canopy performance (e.g.Kirschbaum et al., 1994;Medlyn, 1996;Kruijt et al., 1999;Luo et al., 2001), the model cannot be used to investigate elevated [CO 2 ] impacts in drought conditions because it is not possible to calculate how elevated [CO 2 ] modifies soil moisture.
The SPA model (Williams et al., 2001a, b) is another well-known process-based model, with a focus on the impacts of water availability on forest canopies.However, this model is also limited in that it assumes a horizontally homogenous canopy and thus does not allow for applications to individual trees.We added to MAESTRA the detailed soil water balance routines from the SPA model (Williams et al., 2001a, b), thus combining the strengths of the two models into a new soil-plant-atmosphere model that accounts for non-homogenous stand structure.In addition, we made a number of improvements to the combined model that allow a detailed physiologically-based analysis of the complete soilplant-atmosphere pathway.We believe that the new, combined model is a significant advance over either MAESTRA or SPA and will prove a useful tool for researchers in forest ecophysiology.Our goals in this paper are to fully document the new model and to demonstrate its use by applying it to a previously intractable question, namely the interaction between atmospheric [CO 2 ] (C a ) and drought.
Numerous experiments have been carried out on the interaction of C a and drought on plant functioning (see reviews by Rogers et al., 1994 andWullschleger et al., 2002).A simple prediction for this interaction is that the C a response will be higher under soil drought for two reasons: firstly, lower leaflevel water use under eC a (Eamus, 1991;Medlyn et al., 2001) leads to higher soil water content, which helps plants avoid soil water deficits; and secondly, lower intercellular CO 2 concentration (C i ) during drought enhances the C a response because of the non-linear response of photosynthetic rate to C i (Grossman-Clarke et al., 2001;McMurtrie et al., 2008).While both mechanisms have been confirmed in a number of crop and grassland studies (e.g.Morgan et al., 2004;Rogers et al., 1994), experiments with potted tree seedlings, trees in whole-tree chambers or entire canopies in FACE experiments have yielded highly variable results (Gunderson et al., 2002;Nowak et al., 2004;Warren et al., 2010;Duursma et al., 2011).Two mechanisms are likely responsible for this variability in experimental outcomes: feedbacks from eC ainduced changes in plant size on water use and soil water balance; and acclimation of leaf physiology to long-term growth at eC a .We first study the response of plant carbon uptake and water use to a simulated dry-down, which is a useful baseline for comparison with actual experimental outcomes.We then study the effects of acclimation of two leaf physiology parameters on the C a × drought interaction.Finally, we apply the MAESPA model to a C a × drought experiment where a feedback of increased leaf area in the eC a treatment confounds the direct effect of C a , and show how analysis of experimental data within a model-based framework can extend basic empirical results.

Overview
In the following, we describe the MAESPA model, which is largely the product of the above-ground components of the MAESTRA model (Wang and Jarvis, 1990;Medlyn et al., 2007) and the water balance components of the SPA model (Williams et al., 2001a, b), with several modifications and additions (see Table 1).The MAESTRA model has its roots in

Leaf photosynthesis model
Leaf-level CO 2 assimilation (A), transpiration rates (E), leaf water potential ( ) leaf water potential ( L ) Fig. 1.Flowchart of MAESTRA, the above-ground model of the MAESPA model.Radiative transfer is calculated to a number of gridpoints (typically 72) in each target tree, which is used to drive the stomatal conductance and photosynthesis submodels.These leaf-level rates are then used to estimate whole-stand water use and carbon uptake (see Sect. "Total canopy transpiration").
an early study on radiative transfer by Norman and Welles (1983), and solidified by Wang and Jarvis, (1990).Since then, many improvements have been made, in particular the leaf gas exchange calculations (Medlyn, 1996(Medlyn, , 2004;;Medlyn et al., 2007), and the overall organization and dissemination of the code (see http://www.bio.mq.edu.au/maestra/,hereafter referred to as "the MAESTRA website").Our goal was to widen the applicability of the MAESTRA model by including detailed soil water balance routines and plant hydraulic constraints on gas exchange.
A basic overview of how MAESTRA estimates H 2 O and CO 2 exchange is given in Fig. 1.The processes included in the water balance sub-model are illustrated in Fig. 2. A full list of symbols is presented in Appendix A.
The MAESPA model, like MAESTRA and SPA, runs at typically a half-hourly time-step, although it can also be run at hourly or shorter arbitrary time-steps (up to every minute).

Canopy processes
Here, we provide a brief description of the canopy processes represented in MAESPA, which are largely unchanged from the original MAESTRA model (with the exception of stomatal conductance).Our goal is not to provide an in-depth description of the entire MAESTRA model (see Wang and Jarvis, 1990;Medlyn, 2004;Medlyn et al., 2007 and the detailed manual on the MAESTRA website).

Radiative transfer
The radiation routines are described in detail by Wang and Jarvis (1990).The canopy consists of individual tree crowns, which are described by a basic shape (one of several shapes, including ellipsoids, cylinders and cones), length, height to crown base, and width (in x and y directions).Radiation calculations are performed only for a set of target crowns specified by the user (see Sect. "Total canopy transpiration").A number of grid points are located in a target crown, and radiation (PPFD, NIR and long-wave) at those grid points is  Williams et al. (2001a, b) for the sources of those components, or the corresponding sections in the main text.

Model component Source
Radiative transfer MAESTRA Leaf energy balance " Leaf photosynthesis " Stomatal conductance (g s ), leaf and canopy transpiration " Additional models for g s Tuzet et al. (2003);  Campbell (1974) calculated based on shading within the crown, shading by neighbouring trees, the location of the sun, and whether radiation is direct or diffuse.Scattering of radiation is approximated following Norman (1979).Leaf area within crowns is assumed to be distributed randomly, or to follow a betadistribution in horizontal and/or vertical directions (Wang et al., 1990).At each grid point, leaf area is separated into sunlit and shaded leaf area (Norman, 1993), and the coupled stomatal conductance -photosynthesis model is run separately for each fraction (Fig. 1).

Photosynthesis
For each grid point in the crown, leaf net photosynthesis is modelled using the standard model of Farquhar et al. (1980).
where A c is the gross photosynthesis rate when Rubisco activity is limiting, A j when RuBP-regeneration is limiting, and R d the rate of dark respiration.A c and A j are saturating functions of the intercellular CO 2 concentration where * is the CO 2 compensation point without R d , and k 1 and k 2 are different parameter combinations for A c and A j .The details of these functions and the temperature dependence of the parameters are described elsewhere (e.g.Medlyn et al., 2002).

Stomatal conductance
For each grid point in the crown, leaf-level stomatal conductance to H 2 O (g s ) is modelled using a Ball-Berry type approach (Ball et al., 1987) (Eq. 2).
where A n is the leaf net assimilation rate (µmol m −2 s −1 ), g 0 the conductance when A n is zero, g 1 an empirical parameter, f (D) a function of the leaf-to-air vapour pressure deficit (D) and C s the CO 2 concentration at the leaf surface (µmol mol −1 ) (which is corrected for boundary layer effects, but this is usually a small correction).
We have implemented a number of options for the f (D) function, including the Ball-Berry model (f (D) = 1/ RH, where RH is relative humidity) (Ball et al., 1987), the Leuning model (f (D) = 1/(1 + D/D 0 ) (Leuning, 1995), and the model of Medlyn et al. (2011), who showed that f (D) = 1/ √ D follows from the assumption that g s varies to maintain δA/δE constant (cf.Cowan and Farquhar, 1977).The fourth option, which is used in the simulations shown in this paper, is a modified form of the Tuzet et al. (2003) model.This model takes into account the observation that the response of stomatal conductance to D and soil water deficit is controlled by the leaf water potential ( L ) (Comstock and Mencuccini, 1998).Because L depends on the soil water potential ( S ) as well as the transpiration rate (which depends on D), L summarizes the effect of both D and S on g s (Franks, 2004).Tuzet et al. (2003) use a flexible L modifier (in place of the f (D) function in Eq. 2), where s f and f are parameters; f is a reference water potential, and s f is the "steepness" of the response of f L to L .We assume a steady-state in water flow from the soil to the leaf, so that L at the gridpoint can be calculated from Ohm's analogy to water transport, where E L is the leaf transpiration rate (mmol m −2 s −1 ), itself dependent on D (and boundary layer effects, see Sect. "Calculation of leaf temperature and leaf water potential"), and k L the leaf-specific hydraulic conductance of the soil-to-leaf pathway (see Sect. "Hydraulics of the soil-to-leaf pathway").
Finally, the coupled model of g s and A can be solved by using the basic diffusion equation, where g C is stomatal conductance to CO 2 (= g s /1.6).The solution of the system of equations (Eqs. 1, 2 and 5) can be written as a quadratic equation, which yields the C i (see Wang and Leuning, 1998).We solve the Tuzet model (combination of Eqs.2-5) numerically, by finding the L that is the solution both to Eqs. ( 3) and ( 4).This way, a unique L is calculated for each grid point in the crown.Variation between trees arises due to shading effects, and k L is allowed to vary from tree to tree (as specified by the user).Also, k L is either assumed to be constant within the crown, or to follow a function of height above the ground as specified by the user.
It should be noted that Tuzet et al. (2003) used C i instead of C a in their model, but we use C a to be consistent with Medlyn et al. (2011), and to allow for much more straightforward numerical model solution.The argument by Tuzet et al. (2003) to use C i instead of C a is that stomata appear to respond to C i , not C a (Mott, 1988).However, we note that C i is still implicit in Eq. ( 2), because A n depends on C i through Eq. ( 1), and A n and g s are coupled by C i through Eq. ( 5).

Calculation of leaf temperature and leaf water potential
For the soil water balance, we use the Penman-Monteith equation applied at the canopy level (see Sect. "Total canopy transpiration"), to arrive at consistent calculations for stand energy balance, and boundary layer calculations for soil surface and canopy.At the leaf-level for each grid-point in the crown, transpiration (E L ) is calculated to yield L that is the solution to the Tuzet model of stomatal conductance (Eqs. 2 and 3), and to provide estimates of leaf temperature throughout the crown.For each grid point, E L is calculated with the Penman-Monteith equation (Eq.6), which takes into account boundary layer effects on the exchange of water vapour between leaves and the surrounding air.An iterative scheme is used that finds the leaf temperature that closes the energy balance of the leaf (following Wang and Leuning, 1998).Full details are provided in Medlyn et al. (2007).

Total canopy transpiration
In MAESPA, the soil water balance is calculated using the spatially-averaged transpiration rate by the stand.This spatial average is calculated for a sample of trees specified by the user (the "target trees"), to limit computing time.An example of a selection of target trees is given in Fig. 3.We have chosen not to simulate spatial variation in soil water content based on water uptake of the single trees in a stand, because such a model would be very difficult to parameterize and test.
This approach requires that we scale up estimates of transpiration to the stand level, and estimate total global (shortwave + longwave) radiation reaching the forest floor (for the soil heat balance calculations, see Sect.2.1.4).Absorbed radiation (photosynthetically active radiation (PAR), nearinfrared radiation (NIR) and thermal) is summed for each tree, and corrected for the difference in leaf area per tree for the sample trees, and that for the entire stand.Using the stand density, estimates of global radiation (W m −2 ) incident on the soil surface are obtained.Next, the average whole-tree conductance (mol tree −1 s −1 ) is calculated across the sample trees (weighed by their total leaf areas), and once again corrected for the difference in tree leaf area between sample trees, and the average tree in the entire stand.Stand transpiration is re-calculated using this average canopy conductance, and the intercepted radiation per tree, using the Penman-Monteith equation, where λ the latent heat of water vapour (J mol −1 ), E T is canopy transpiration (mol m −2 (ground) s −1 ), s the slope of the relation between saturation vapour pressure and temperature (Pa K −1 ), R n net radiation (W m −2 ), g B the boundary layer conductance, and g V the total conductance to water vapour.The boundary layer conductance (g B , mol m −2 s −1 ) for the canopy is calculated following Jones (1992), where k V von Karman's constant, u z the wind speed measured at a height of z H , and c converts to molar units.The parameters z D and z 0 are related to the extinction of wind speed above and within the canopy: z 0 is the "roughness length" related to the roughness of the canopy surface, and z D is a reference height.Both these parameters may be estimated from canopy height (Jones, 1992, p. 68), or from more detailed methods that account for stand structure (Schaudt and Dickinson, 2000).
The total conductance to water vapour is, where g C the canopy conductance (mol m −2 s −1 ), which is estimated by averaging total conductance (mol tree −1 s −1 ) across the sample trees weighed by their total leaf areas, and multiplying by stand density (m −2 ).

Overview
The soil water balance is calculated for a horizontally homogenous soil compartment, separated vertically in a number of layers (typically ca. 10, specified by the user).The advantage of a multi-layer soil model is that the root density distribution can be specified, as well as vertical variation in soil texture.Furthermore, experiments have shown that plant water potential tends to equilibrate with the wettest zones in the soil, not the average soil water potential (Schmidhalter, 1997).Only a multi-layer soil model can reproduce this observation.
In MAESPA, the soil water storage in each of the horizontally homogenous layers (s i , mm) is calculated from infiltration (I i ), drainage (D i ), root water uptake (E i ) and soil evaporation (E s,i ) (Eq. 9), at the same time-step as the aboveground processes (typically, half-hourly).Soil evaporation draws water only from the top layer.Drainage out of the lowest root layer is defined as deep drainage.

Infiltration
The rain that reaches the soil surface is assumed to infiltrate into more than just the top layer.This is accounting for rapid soil water flow through macropores (accounting for increases in soil water content at depth after heavy rainfall that cannot be accounted for by matric drainage rates).We use a simple exponential function (taken from the BROOK90 model, (Federer et al., 2003) (Eq. 10).
where φ is an infiltration parameter (0-1), z i the depth to the bottom of layer i (m), and Z the total soil depth (m).If φ = 0, all infiltration occurs in the top soil layer (no macropore flow).If φ = 1, all layers receive equal infiltration.

Root water uptake
Total water uptake is distributed among soil layers according to the fine root density and soil water potential in those layers.The fraction of roots in each layer can be specified manually, or assumed to follow Eq. ( 11), which is useful because Jackson et al. (1996) have summarized a large database on root distributions with this equation.
where F R is the cumulative fraction of fine roots to depth z (m), and β is a parameter that specifies the shape of the distribution (see Jackson et al., 1996).Soil depth is multiplied by 100, to be consistent with Jackson et al. (1996) who specified z in cm.We have modified the original equation so that the total fine root length (an input parameter) is reached at the maximum depth of the rooting profile.
The root system can be viewed as a combination of resistances that are coupled in parallel, consisting of a soil-toroot surface resistance (R sr ) and a resistance to water uptake by the root itself (the radial resistance, R rad ).Because water transport is largely passive, following gradients in water potential from soil to roots, the relative water uptake by different layers in the soil follows from the partitioning of the resistances to water uptake in these layers.Generally, water uptake in a soil layer (E i ) is, where S,i is the soil water potential in layer i, R,i is the root xylem water potential, R sr the soil-to root resistance, R lgi the longitudinal resistance to water flow, and R rad the radial resistance to water uptake (across the root epidermis to the xylem).This is a very difficult equation to solve simultaneously for multiple soil layers, because R,i , R rad and R lg are typically unknown, and could vary with depth of the layer (especially R lg,i , as the path length through the root becomes longer).To solve Eq. 12, we follow the assumptions by (Taylor and Keppler, 1975) who showed in an elegant experiment that (1) R lg is very small compared to R rad , so that R,i can be taken as a constant, and (2) R sr is small compared to R rad , except in very dry soil.Because R rad is inversely proportional to the total fine root length in a soil layer, it follows that, where R is a mean root water potential, and L v,i the fine root density (m m −3 ) in layer i.The constant of proportionality may be found from the fact that E i = E T , where E T is the total water uptake.Note that we do include R sr in the calculations of the overall resistance of the soil-to-leaf pathway (see Sect. "Hydraulics of the soil-to-leaf pathway"), but omit it here to simplify the fractional uptake of water from different soil layers.The mean root water potential is taken as the minimum value allowed for roots ( Rmin ), an input parameter; so that no uptake is possible in soil layers where the soil water potential is lower than Rmin ).When R > S , it is assumed that water does not flow from roots back into the soil, that is, we do not consider hydraulic redistribution in our model (see Domec et al., 2012).

Soil evaporation
The evaporation of water from the soil surface is estimated with a physical-based model developed by Choudhury and Monteith (1988), and modified by Williams et al. (2001a).Soil evaporation draws water only from the top layer, and assumes that water has to travel through a thin dry layer at the soil surface.
The rate of evaporation (mm) is determined from the difference between water vapour pressure in the soil pore space and the air above, and a conductance to vapour transfer (Eq.14), where G s,t is the total conductance from the soil air space to the air above the boundary layer (m s −1 ), e a the partial water vapour pressure of the air (kPa), and e s that of the soil pore space (kPa).The term k 1 converts from pressure units to volumetric units.G s,t is the total conductance, including the path through the soil (G ws ), and that through the boundary layer just above the soil surface (G am ).The latter is estimated following Eq.( 7) (Sect."Total canopy transpiration"), with the reference height z D set to zero.Next, G ws is estimated from the diffusivity of water vapour, which is soil temperature dependent, and the tortuosity of the soil air pathway, which determines the effective path length for vapour transfer.
where D eff is the effective diffusivity, ω s a tortuosity parameter, θ 1 the pore fraction of the top soil layer, L d the thickness of the dry layer at the soil surface, and D w diffusivity of water vapour (calculated following Jones (1992) accounting for soil temperature).The rationale for Eq. ( 15) (Choudhury and Monteith, 1988) is that G ws is proportional to the air space in the soil, but inversely proportional to the distance water vapour has to travel.This distance is a function of both the thickness of the dry layer (L d ), and the tortuosity of the dry layer (Eq.16).L d depends on the timing of the last rainfall, and the rate of soil evaporation, but is not affected by plant water uptake.Following Williams et al. (2001b), MAESPA keeps track of multiple dry layers to account for complex dynamics with short intermittent storms.A minimum value (L d,min ) is specified as a parameter to prevent very high rates of soil evaporation in wet soil.
The partial pressure of water vapour in the soil pore space is calculated from, where e sat is the saturated vapour pressure (calculated from temperature following Jones 1992), s,1 the soil water potential in the surface layer, V w the partial molal volume of water (m 3 mol −1 ), R the gas constant, and T s the soil surface temperature (K).

Canopy interception
The Rutter et al. (1975) where E w the wet evaporation rate (mm t −1 ), r 1 the free throughfall fraction (0-1), and r 2 and r 3 are canopy drainage parameters.This differential equation is integrated with a Runge-Kutta method to obtain canopy throughfall, canopy storage, and wet evaporation rates.

Drainage
The vertical drainage of soil water is estimated directly from the soil hydraulic conductivity, which itself is a function of the soil water content.After adding the infiltration of rainfall to the soil water content for each layer, the water balance for each soil layer becomes, where D i is the drainage from layer i.Because water is assumed to not travel upwards in the soil, this equation can be solved for one layer at a time, starting at the top and moving downward.Drainage is simply equal to the soil hydraulic conductivity (m s −1 ) in that layer, which is calculated from θ i .The system of equations is solved with a Runge-Kutta integrator (following the SPA model, Williams et al., 2001a, b).

Hydraulics of the soil-to-leaf pathway
We use the simple set of equations developed by (Campbell, 1974) for the dependence of soil water potential and the hydraulic conductivity on the soil water content.Soil water potential ( S , MPa) is given by Eq. ( 20), where θ the soil volumetric water content (m 3 m −3 ), θ sat the soil porosity, and e and b are soil texture dependent parameters (see Cosby et al., 1984).Each of these parameters can vary by soil layer (i subscripts are omitted for clarity).The soil hydraulic conductivity is estimated from Eq. ( 21), where K s the conductivity (mol m −1 s −1 MPa −1 ), e and b are the same parameters as in Eq. ( 20) (using the fact that conductivity and water potential are physically related), and K sat is the saturated hydraulic conductivity (mol m −1 s −1 MPa −1 ).
The soil-to-root resistance (MPa s m 2 (ground) mol −1 ) is estimated with the single root model of Gardner (1960) (see Williams et al., 2001a andDuursma et al., 2008 for more details).This model estimates the effective path length for water transport through the soil matrix to the root surface from the fine root density.The equation is, where r r the mean root radius (m), L v the total fine root density (m m −3 ) in the soil layer, H s the height of the soil layer, and r s the mean distance between roots (1/ √ π L V ).The total resistance for all soil layers combined (R sr,t ) is estimated by assuming that the resistances are coupled in parallel.The total leaf-specific hydraulic conductance (k L , mmol m −2 (leaf) s −1 MPa −1 ) from soil to leaf can now be found as, where L T total canopy leaf area index (m 2 (leaf) m −2 (ground)), and k P the plant component of the leaf-specific hydraulic conductance (mmol m −2 s −1 MPa −1 ), which is typically estimated from measurements of E L and L under well-watered conditions (see, e.g.Delzon et al., 2004).It is assumed that k P includes the root components of the plant pathway (defined in Eq. 12 by the resistances R rad and R lg ), and that these components do not change with the plant water potential.
Following Williams et al. (2001a), we calculate a weighted soil water potential for use in the calculation of L and stomatal conductance.The s in each layers is weighed by the maximum water transport possible in that layer, depending on R sr in that layer and a minimum root water potential ( R,min ),

Overview
The components of the soil surface heat balance are calculated to arrive at the soil surface temperature, and the vertical gradient of soil temperature.The soil surface temperature affects only the soil evaporation (see Sect. "Soil evaporation").If soil evaporation is not of interest, or can be assumed negligible, the soil heat balance need not be calculated (thereby simplifying the parameterization).

Soil surface temperature
The soil energy balance is the sum of four heat fluxes: net radiation (Q n ), soil surface evaporation (latent heat flux) (Q e ), soil heat transport (to deeper soil layers) (Q c ) and sensible heat flux (Q h ).Due to conservation of energy, these fluxes sum to zero at any time: All the components in Eq. ( 25) depend on the soil surface temperature (T s,1 ).The MAESPA model finds the soil surface temperature that provides closure in the soil energy balance.
Net radiation on the soil surface (Q n ) is global radiation (solar + downward thermal) minus long-wave radiation (Q L ) emitted by the soil surface.The latter depends on the soil surface temperature (Eq.26).
where ε is the emissivity (assumed to equal 0.95), σ the Stefan-Boltzmann constant (W m −2 K −4 ) and T s is in K. Global down-welling radiation is the sum of short-wave radiation (i.e.solar radiation minus that intercepted by the canopy), and long-wave radiation emitted by the canopy.
The soil latent heat flux (Q e ) is calculated from soil evaporation (see Sect. "Soil evaporation"), and the latent heat of evaporation (a function of temperature; Jones, 1992).
Soil heat transport (Q c ) follows from the difference in soil temperature between the first and second layer, where K th is the the soil thermal conductivity (W m −1 K −1 , see Sect."Soil surface temperature"), and z 1,2 the depth difference between the second and first layer (m).
Sensible heat flux (Q h ) is calculated from the difference between air temperature and soil surface temperature, and the total conductance to heat (Eq.28) where c p the heat capacity of air (J kg −1 K −1 , a constant), ρ the (temperature-dependent) density of air (kg m −3 ), G a,m the boundary layer conductance to heat transfer (assumed to be equal to that of water vapour transfer, see Sect."Soil evaporation"), and T s,1 -T air the temperature gradient between soil surface and the air.

Soil temperature profile
The transport of heat down the soil temperature gradient is calculated following the SPA model (Williams et al., 2001b).The Fourier heat transport equation is solved using the Crank-Nichols scheme (Press et al., 1990) resulting in a soil temperature profile and corresponding heat fluxes between soil layers.Inputs for this routine are the soil thermal conductivity (by layer), and soil heat capacity.
We used an empirical model developed by Lu et al. (2007) to estimate the soil thermal conductivity (in W m −1 K −1 ) from soil porosity, water content, temperature and organic matter content.Their model takes slightly different parameter values for fine and coarse textured soils; we used the soil texture parameter (b in Eq. 20) to determine which parameter set to use (based on their Table 1, we used b = 5.3 below which soils are considered "coarse textured").With this method, it is straightforward to set the top layer of the soil as a "litter layer" that effectively insulates the soil through a very low thermal conductivity.
The heat capacity of the soil is determined by separating the soil into solid (quartz) and water fractions, and finding the weighted average of their heat capacities (cf.de Vries, 1963; see also Ogée et al., 2001), and assuming that the soil air fraction has negligible heat capacity.

Parameterization
In MAESPA, the water balance is calculated for a horizontally homogenous soil, but the canopy consists of a collection of single trees.Canopy transpiration is estimated based on transpiration by the "sample trees", and it is therefore vital that these sample trees represent the canopy in terms of water use.It is recommended that MAESPA is run for a large number of sample trees, for example in an arrangement shown in Fig. 3.For all trees in the stand, estimates of leaf area, crown shape, crown width, crown length and height to crown base are needed.Stands may also consist of just one tree, in which case the product of plot size and rooting depth can be interpreted as the soil volume, or pot volume, available to the tree.
Environmental drivers need to be specified on a (half-)hourly time-step, and include air temperature, solar radiation, precipitation, relative humidity and optionally wind speed and CO 2 concentration.For the water balance, crucial parameters are total rooting depth and the soil water retention curve.The latter can be estimated from soil texture, and equations summarized by e.g.Saxton and Rawls (2006) and Cosby et al. (1984).A brief summary of the list of parameters needed to run the water balance component of MAESPA is given in Table B1.

Implementation and batch utility
MAESPA is written in Fortran, with simple text-based input and output files.An R (R Development Core Team, 2011) package is available, Maeswrap, which aides sensitivity analysis and other computer-intensive simulation studies.This also includes a utility to graph the stand in 3-D (Fig. 3 was produced with the Maeswrap package).The compiled model and the source code are available on the MAESTRA website (see Sect. 2.1).

Application of MAESPA: a case study on the interaction between C a and drought
In the following, we present three brief case studies that illustrate the application of the MAESPA model to studying complex interactions between environmental drivers and plant parameters.All case studies analyse the interaction between atmospheric CO 2 concentration (C a ) and drought.

Dry-down simulations
We simulated the response of several plant variables to a drydown in order to establish a clear picture of the baseline expectation of the interaction between C a and drought, that is, the interaction when C a -related feedbacks and acclimation are ignored.We used a hypothetical stand with total leaf area index of 3.8 m 2 m −2 , with a default leaf physiology parameter set (see Table 2 for parameters and their values), and the Tuzet model of stomatal conductance (Eqs. 2 and 3).Ambient C a was set to 380 ppm, and elevated C a was chosen as 620 ppm (Barton et al., 2010).Weather data was generated for a typical sunny summer day, using the Bristow and Campbell algorithm (Bristow and Campbell, 1984), and was the same for each day during the dry down, as was the solar angle (although it did change during the day).

Effect of acclimation of leaf parameters on C a × drought interaction in water use and CO 2 uptake
Acclimation of leaf physiology to long-term growth at eC a is one possible explanation for why many real-world experiments deviate from baseline expectations of C a × drought interactions.We test the impact of acclimation of two leaf physiology parameters, k L and f .A number of studies have found decreases in the leaf-specific hydraulic conductance (k L ) in plants grown in eC a (Atkinson and Taylor, 1996;Eamus et al., 1995;Eguchi et al., 2008;Heath et al., 1997), with reductions in a wide range of 10-100 %.Berryman et al. (1994) found a higher sensitivity of g s to decreasing water content of excised leaves in Maranthes corymbosa, which can be interpreted as a less negative f in the Tuzet model (Eq.3).This observation is also consistent with a higher sensitivity to abscisic acid (ABA) in eC a found in a number of species (Dubbe et al., 1978;McAdam et al., 2011), and the fact that ABA concentrations increase in low L (Pierce and Raschke, 1981).To test the sensitivity of the C a × drought interaction to these parameters, we ran two additional dry-down simulations, one with a 50 % reduction in k L at eC a , and one with a 0.4 MPa increase in f .Both these changes in parameter values are within observed ranges of change with long-term growth at eC a , but are arbitrary and only chosen to illustrate the effect on the C a × drought interaction.All other settings and parameters were as specified in the dry-down simulation (Sect.2.2.1).

Drought × C a interaction in cherry seedlings
To illustrate the importance of whole-plant feedbacks in treatment responses, we applied the model to a C a × drought experiment.Centritto et al. (1999a, b) describe an experiment where cherry seedlings were grown in ambient C a concentration (aC a ) (350 ppm) and elevated C a (eC a ; 700 ppm) treatments in well-watered conditions until half the plants were subjected to a dry-down.The analysis of their results was complicated by the fact that total leaf area of eC a seedlings was higher than that of aC a seedlings, compensating for lower water use per unit leaf area, so that total water use was similar between treatments.We re-analysed their dataset in a model-based framework where we can integrate effects of leaf area, C a , and leaf-level physiology parameters on wholeplant interactions between C a and drought.
We estimated MAESPA parameters for the cherry seedling based on the published information as much as possible, and found the remaining parameters by fitting to observed data (see below).Details and estimated parameters are presented in Table 3.We constructed a weather dataset based on the latitude of the study and the reported mean air temperature and relative humidity.Daily incident PAR was estimated from air temperature using the Bristow and Campbell algorithm (Bristow and Campbell, 1984).Using the fitted parameter set for aC a plants, we simulated water use by the eC a plants, with the only difference that leaf area was increased as observed in the experiment.

The interaction between C a and drought
A simulation was carried out to establish baseline behaviour of MAESPA during a dry-down.
At ambient C a , simulated total water use (E T ) declines earlier and more rapidly than total carbon uptake (A T ) (Fig. 4a and b), implying that C i /C a declines (Fig. 4c) and A T /E T increases as the dry-down progresses (Fig. 4d).As the soil water content declines (Fig. 4e), midday leaf water potential ( L ) decreases steadily, and continues to decrease because of cuticular water loss (Table 2).As g s decreases during the dry-down, the difference between leaf and air temperature increases (Fig. 4g), and the depth of water uptake gradually shifts to deeper layers (Fig. 4h).
These simulations also summarize our baseline expectations, in the absence of feedbacks or acclimation, for a C a × drought interaction.Under eC a , E L is initially lower (Fig. 4a), which leads to less negative L (Fig. 4f), and a higher soil water content (Fig. 4e).The daily integrated transpiration efficiency (A T /E T ) is higher under eC a , and increases more rapidly under eC a as compared to aC a (Fig. 4d).This latter prediction is sensitive to the assumed value for g 0 (the cuticular conductance), so that a lower g 0 leads to Fig. 4. Simulation of the effect of a dry-down on whole-plant fluxes and water balance under ambient and elevated C a .The simulation was performed for a single tree with 35 m 2 leaf area in a stand of identical trees.Fluxes are expressed as averages over the canopy of that single tree.See Table 1 for the parameter set used in the simulation.Shown are the decrease in total water use (E T ), total daily net photosynthesis (A T ), midday C i /C a , the daily integrated transpiration efficiency (A T /E T ), volumetric soil water content (θ), midday leaf water potential for the sunlit leaves ( L ), the difference between leaf and air temperature (T leaf -T air ), and the depth of root water uptake (d 50 , the depth above which 50 % of water is taken up).For panels (a)-(d), the ratio of eC a to aC a is shown with a grey line.Note the difference in scale for this ratio (right y-axis) for panel (b).
a more pronounced increase in A T /E T as the drought progresses (not shown).Both E T and A T show a three-phase response when expressed as the ratio eC a to aC a (Fig. 4a and  b).At first, when both plants have sufficient water, the ratio is constant.The ratio then increases due to higher soil water content in the eC a treatment.When this saved water store is exhausted (by day 40-45), the ratio declines, as both E T and A T are increasingly controlled by the cuticular conductance.MAESPA predicts a strong positive interaction between C a and drought for two reasons.The first is the "water savings" effect described above: photosynthesis of elevated C a trees is high for longer due to higher soil moisture.Secondly, even when both treatments are at the same soil water content, photosynthesis is stimulated more by eC a in dry soil compared to wet soil (Fig. 5a, see also Fig. 4).This response arises because C i declines as the soil dries out (Fig. 4c), and A n is more sensitive to C i at low values of C i , due to the saturating response of A n to C i .This "C i effect" is demonstrated in Fig. 5a   Leaf area (start-end) 0.06-0.08m 2 (aC a ), 0.1-0.115m 2 (eC a ) Calculated from Centritto et al. (1999a), their Table 5 and Fig.  shows total canopy transpiration rate (E T ) for the same simulations.The solid line shows the interaction when all plant parameters are unchanged due to growth at eC a ("no acclimation").Note that in very dry soil, transpiration rates are similar in aC a and eC a , but the stimulation of A T due to eC a is much higher than in wet soil (a positive interaction of drought and eC a ).These effects are much reduced for a higher sento L in the stomatal conductance model (" f +0.4 MPa"), which leads to a negative interaction between C a and drought (the C a effect is less in dry soil than wet soil, for both E T and A T ).A lower plant hydraulic conductance (k P ) reduces E T in wet soils, and greatly reduces the positive interaction of C a and drought that was found without acclimation.
We studied the effect of acclimation of two leaf physiology parameters, the leaf-specific hydraulic conductance (k P ) and the sensitivity to L ( f , Eq. 3), on the "C i effect".When k P is reduced by 50 %, the strong positive interaction with soil water content as observed in the baseline simulation is much reduced (Fig. 5a), so that the eC a stimulation of A T is only moderately dependent on soil water content.When f is increased by 0.4 MPa, the interaction even reversed, so that A T decreases with eC a in very dry soil (Fig. 5a) and E T is reduced by eC a in very dry soil (Fig. 5b).This occurs because a higher f leads to reduced stomatal conductance in eC a compared to aC a as the L declines in dry soil.
Plant size may change following long-term growth at eC a , which can feed back to drought responses and modify C a × drought interactions.We simulated a C a × drought interaction in an experiment where plant leaf area increased in response to eC a , fully compensating the lower water use per unit leaf area.MAESPA was able to simulate the decrease in soil water content and L (Fig. 6a and b), after calibrating parameters related to stomatal conductance (Table 3) using the aC a treatment only.Using this parameter set, however, the measured L was overestimated early in the dry-down for the eC a treatment, and underestimated late in the drydown (Fig. 6b).The relative response of E T to eC a (Fig. 6) illustrates that the eC a plants used more water early in the dry-down, which led to a more severe water stress, so that water use was substantially less under eC a towards the end of the dry-down.This response was also predicted by MAESPA (Fig. 6c), because L was lower in the (simulated) eC a treatment.
Using the calibrated model, it is possible to tease apart contributions of different mechanisms on the overall interaction between C a and drought, by running simulations with different settings.For this experiment, we ran simulations exploring the contributions of changes in leaf area vs. changes in leaf-level water use.If plant leaf area was assumed unchanged between aC a and eC a treatments, we observed a very strong positive interaction between C a and drought on E T , as expected from the baseline simulations (Fig. 6c, dashed line).If leaf area was assumed to increase in the eC a treatment, but leaf-level water use was assumed not to change (dot-dashed line in Fig. 6c), E T declined much more quickly, indicating that the leaf-level response to eC a did have a substantial ameliorating effect on the response of E T to drought.
Finally, we quantified the drought impact on total photosynthesis in the cherry experiment, by calculating total photosynthesis over the entire dry-down under different assumptions, and expressing it relative to simulated total photosynthesis under well-watered conditions (Fig. 7).At aC a , drought reduced total photosynthesis by 22 %.For the eC a treatment, three simulations are summarized.The first two are calculated drought responses if leaf area is the same between aC a and eC a treatments, and the third is the model prediction for the actual experimental conditions of the cherry experiment.The first simulation, Run 1, is the "baseline", which includes both the "water savings" and the "C i effect" (see above).The second simulation, Run 2, is at aC a , but with reduced leaf area to match the pre-drought eC a water use, and therefore is equivalent to the water savings effect of eC a .In the third simulation, Run 3, leaf area was increased similarly as in the cherry experiment, so that pre-drought water use was the same as in the aC a .In the eC a simulations, drought reduced total photosynthesis by 10 % (without feedbacks; Run 1), and 13 % (Run 2, water savings only), respectively.With the leaf area feedback (Run 3), there was a larger reduction in total photosynthesis (30 %), which counteracted the positive interaction between C a and drought.In summary, there was a positive eC a × drought interaction, because drought reduced photosynthesis less in eC a than in aC a .This positive interaction was largely the result of the "water savings" effect, and disappears completely when leaf area is increased in eC a , as was the case in the cherry experiment.

Discussion
We have presented a new soil-plant-atmosphere model, MAESPA, that can be applied to both individual plant and whole stand scales.The model includes detailed radiation transfer and leaf physiology routines from the MAESTRA model, and mechanistic water balance and hydraulics from q q q q q q q q q q q q q q q q q q 0 10 20 q q q q q q q q q q q q q q q q q q (b)  (1999a, b).Model parameters were based on reported values in the original study, or calibrated to yield a satisfactory fit to the data of the aC a treatment only (see text and Table 2).The optimized parameter set was then used to predict water balance in the eC a treatment, and taking into account the observed increase in leaf area in the eC a seedlings.(a) Decline in average soil water content over the rooting zone (θ) as the drought progressed.(b) Decline in the midday leaf water potential for sunlit foliage ( L ).Note a relatively poor fit for the eC a treatment: L is over-predicted early in the drought, and under-predicted towards the end of the drought.(c) The ratio of total water use in eC a to aC a during the dry-down.Higher leaf area in eC a initially leads to higher water use, but this leads to lower θ (a) and L (b), so that eC a seedlings were more water-stressed toward the end of the drydown than their aC a counterparts.The solid line shows the simulations where eC a and the increase in leaf area in the eC a treatment were taken into account, dashed lines show either the direct C a effect only (without leaf area feedback) or the leaf area feedback only (no eC a treatment).Effect of the drought treatment on simulated total photosynthesis (summed over the entire 47 day simulation) for the dry-down in the cherry experiment (Fig. 6).The aC a simulation shows a 22 % reduction in A T in the dry-down treatment.For eC a , three simulations are shown.Run 1: full simulation but without the leaf area feedback, Run 2: simulation at aC a , but with reduced leaf area to match total water use in the eC a simulation (Run 1), Run 3: including the observed leaf area increase in the eC a treatment.
the SPA model.We have shown that the model gives realistic predictions of the response of several plant variables to drought.As an example application, we used the model to study the interaction between atmospheric [CO 2 ] (C a ) and drought.

Understanding controls on the C a × drought interaction
In this paper, we illustrated the use of the MAESPA model by quantifying interactions between C a and drought under several potential experimental scenarios.First we simulated our "baseline" expectation of the C a × drought interaction, in the absence of feedbacks or acclimation of plant properties to long-term growth at eC a .Experimental outcomes can be compared to simulations like these, in order to evaluate whether the results are in quantitatively in line with current understanding of biophysical and physiological controls on whole-plant gas exchange and water balance.
The MAESPA model was able to simulate the effects of a soil dry-down on several variables in line with published observations.During the dry-down, C i /C a steadily decreased, so that A T /E T increased, which is consistent with published studies where non-stomatal limitations to carbon uptake are minimal (Brodribb, 1996).Midday leaf water potential ( L ) decreased steadily, as typically observed (Sperry, 2000), and leaf temperature increased as a result of lower stomatal conductance in dry soil (Jones, 1992;Triggs et al., 2004).Finally, the depth of root water uptake gradually shifted to deeper layers (cf.Rambal, 1984;Duursma et al., 2011).
The baseline simulations predict that total CO 2 uptake (A T ) is enhanced more by eC a in dry soil (Fig. 5a), which is in line with previous predictions (Grossman-Clarke et al., 2001), and follows directly from the nonlinearity of the dependence of A n on C i (the "C i effect").The magnitude of this drought-enhanced eC a response depends on the parameters used (in particular V cmax , g 1 ), as well as the soil water content (Fig. 5a).Although many studies on agricultural crops have demonstrated that biomass growth or net canopy carbon uptake is more enhanced by eC a during drought (Rogers et al., 1994), a great number of studies, particularly on trees, fail to demonstrate this effect (see Wullschleger et al., 2002;Nowak et al., 2004;Duursma et al., 2011;Warren et al., 2011).It should be noted that the interpretation of biomass growth in drought conditions is not straightforward because biomass growth may be temporally uncoupled from photosynthetic CO 2 uptake (Körner, 2003;Sala and Hoch, 2009).However, this is unlikely an explanation for the lack of observations for expected C a × drought interaction in trees, because many studies report CO 2 uptake as well as biomass growth.We showed here that two plant parameters, that are frequently observed to be affected by acclimation to eC a , can reduce or even reverse this expected interaction (Fig. 6).Some studies have found this "reverse" response: Schäfer et al. (2002) found that, in a FACE study on Pinus taeda L., E L was only reduced by eC a when soil was dry, which counters the baseline expectation of C a responses (Fig. 6).In a FACE study on Liquidambar styraciflua L., Gunderson et al. (2002) found a higher sensitivity of g s to S in eC a trees, which also has the potential to reverse the expected interaction between C a and drought.We need to quantify the effect of such leaf acclimation to eC a and investigate the degree to which it can explain experimental outcomes that diverge from baseline predictions.A model such as MAESPA is an essential tool to quantify the contribution of such mechanisms.
Baseline expectations are that lower leaf-level water use in eC a will lead to a higher soil water content, thus delaying the onset of drought (Morgan et al., 2004).A drought treatment should therefore have less impact on total carbon uptake in eC a than in an aC a treatment, because drought stress is postponed (the "water savings" effect).However, testing this hypothesis against data from C a × drought experiments can be difficult because increased plant size in eC a often compensates for lower leaf-level water use, so that a clear "water savings" effect is often not directly observed (Morison and Gifford, 1984;Roden and Ball, 1996;Centritto et al., 1999a).We applied MAESPA to an experiment in cherry, where increased leaf area fully compensated for lower leaf-level water use.In this experiment, soil water content declined at similar rates in aC a and eC a treatments (Fig. 6), despite lower water use at the leaf-level in in eC a , because leaf area was ca.50 % higher in the eC a plants.The MAESPA model successfully simulated the compensation of total plant water use by increased leaf area, using one parameter set for both C a treatments, and the measured leaf areas.Because of the leaf area feedback, Centritto et al. (1999a) concluded that there was no positive interaction between drought and C a , a conclusion that followed from a standard empirical analysis of the results.Experiments like this can be further analysed in a quantitative framework like MAESPA, because plant leaf area can be quantitatively accounted for in the simulations, which is much more difficult to accomplish in a purely empirical analysis.Using the parameterized model, it is then possible to separate the various contributing factors to the overall C a × drought interaction, and to estimate the strength of the interaction between C a and drought.
Using the parameterized MAESPA model, we showed that there were interactions between C a and drought in the cherry experiment.As the drought progressed, total plant water use declined more rapidly in the eC a treatment (Fig. 6c), an experimental observation that was roughly matched by the model simulation (Fig. 6c).However, a poor fit to L (Fig. 6b) was necessary to match the larger reduction in E T in eC a as the dry-down progressed (Fig. 6c).This mismatch is possibly because L was sampled on a few sunlit leaves that were not representative of the entire canopy.Without the leaf area feedback, there was a strong simulated positive interaction of C a and drought (Fig. 6c): water use could continue much longer in the eC a treatment due to initial water savings.This analysis was therefore able to separate effects of leaf area and leaf-level processes on the response of plant water use to drought and C a .
The interaction between C a and drought is expected to particularly affect plant CO 2 uptake (A T ), because of the "C i effect": photosynthesis is more responsive to C a at low stomatal conductance, as is the case during drought (Fig. 5a).But what is the expected strength of this interaction, and is it more important than the "water savings" effect?For the cherry experiment, we calculated total A T over the entire drying cycle, and expressed it as a ratio of droughted to irrigated control (Fig. 7).Drought reduced total A T by a smaller fraction in the eC a treatment compared to the aC a treatment (10 % vs. 23 %).This relatively small difference is perhaps one reason that the C a × drought interaction in experiments is often not significant, because there may be insufficient power to detect effects of this size.The simulation analysis demonstrated that the positive interaction was mostly a water savings effect in this case, the "C i effect" was very small (Fig. 7).It is possible that the C i effect is larger in other experiments, because it depends on the shape of the A-C i curve, the degree of drought stress, stomatal conductance, and the length of the drought period.
Large scale simulations of C a effects on vegetation water use and carbon uptake do not account for acclimation or feedbacks of plant processes to long-term growth at C a (e.g.Cramer et al., 2001;Luo et al., 2008), and as such yield predictions of a positive C a × drought interaction in line with the baseline predictions shown here (Fig. 4).However, actual experimental outcomes yield varied results, making it difficult to inform model formulation and parameterization with experimental data.Here, we showed that by taking into account the observed feedback of plant leaf area in one experiment, it is possible to study the C a effect on plant water use and carbon uptake, had the feedback not occurred.A model like MAESPA can also be used to evaluate alternative explanations for the deviation from experimental outcomes from the expected theory, such as acclimation of plant hydraulic parameters (e.g.k P , f ), and to evaluate whether responses at the leaf level match the responses at the whole-canopy scale.

Possible applications of MAESPA
While this study focussed on the interaction between C a and drought, there are a number of possible applications of a soilplant-atmosphere model that can be applied to whole-plant and forest stand scales.Analysis of complex experiments where data are collected at leaf-level, whole-tree or canopy level, and in the soil, can be strengthened if all data are integrated in the parameterization of a soil-plant-atmosphere model (Williams et al., 2001b;Medlyn et al., 2005;Duursma et al., 2007Duursma et al., , 2009)).Because MAESPA can be applied to potted plants, it may be used to generalize experimental results from the vast number of experiments on potted plants, that are typically confounded by changes in plant size with experimental treatments (e.g.Damesin et al., 1996).
There is a growing interest in the effects of stand structure on ecosystem functioning, because the spatial distribution of leaf area index (LAI) in sparse or dense crowns affects radiation interception, energy balance, and total water use (Chen et al., 2008;Yang et al., 2010), even though LAI is the primary driver for fluxes of water and carbon.The MAESPA model is well suited to study effects of canopy structure and grouping of foliage in tree crowns on whole-canopy performance, and to evaluate simplified approaches.
All currently available soil-plant-atmosphere models can only be applied to entire canopies, restricting their use to studying stand-level processes.The advantage of MAESPA is that single plants can be studied.For example, models of vegetation water use are typically tested against scaled-up sap-flux measurements (Hanson et al., 2004;Williams et al., 2001a;Zeppel et al., 2008).An individual-based model such as MAESPA can be used to address questions of resource distribution among plants of different size and species within a canopy (cf.Binkley et al., 2010), in particular regarding the use of soil water and response to soil drought.

Uncertainty in process representation
The quantitative understanding of a number of processes in the soil-plant-atmosphere continuum is limited, so that www.geosci-model-dev.net/5/919/2012/Geosci.Model Dev., 5, 919-940, 2012 improvements to models like MAESPA are certainly possible.Below, we discuss uncertainty in three components of MAESPA, but they apply to any soil-plant-atmosphere model: root water uptake, variation in plant hydraulic conductance, and non-stomatal limitations to CO 2 uptake during drought.These three components are just examples where progress in process understanding will improve soil-plantatmosphere models; there are certainly others.
Predicting the distribution of root water uptake with depth in the soil is an old problem (Gardner, 1964), and surprisingly little progress has been made since the simple model advanced by Taylor and Keppler (1975), which is used in nearly all root water uptake models (Feddes et al., 2001).Although this approach seems relatively successful in predicting relative uptake of water from different soil layers (Markewitz et al., 2010), it is not useful in predicting the reduction in total water uptake when only a part of the root system is accessing wet soil, as is the case for chronically droughted trees that have few roots at great depth (Calder et al., 1997).A better understanding of the root hydraulic conductance and how it varies with depth in the soil, and the partitioning of the resistance between radial and longitudinal components of the root pathway are needed to improve on this model component.
It is typically assumed, as in MAESPA, that the decline in CO 2 uptake during drought is the result only of reduced stomatal conductance, which simply limits the diffusion of CO 2 into leaves.However, there is ample evidence that photosynthetic capacity (A n at a given C i ) also declines during drought, albeit highly dependent on species and possibly only during more severe water stress (Lawlor and Cornic, 2002).Recently, Keenan et al. (2009) and Grant and Flanagan (2007) both showed that accounting for the reduction in photosynthetic capacity with drought stress improved model predictions of canopy fluxes.A more general understanding of non-stomatal limitations and how they develop during drought stress will improve models such as MAESPA.
In MAESPA, the hydraulic conductance of the plant pathway (k P ) does not decline during drought, and does not vary among shaded or sunlit portions of the canopy.Although k P typically does decrease during drought due to formation of air-filled vessels (Sperry ,2000), Duursma et al. (2008) showed that a model assuming a fixed k P was successful in predicting the response of plant water use to water limitation, in part because the soil resistance becomes limiting to water transport (Fisher et al., 2006).Nonetheless, gradual reduction in k P during drought is often an important determinant of plant water use (Sperry et al., 1998;Hacke et al., 2000).
Differences in k P between shaded and sunlit leaves in the canopy may exist because of shorter path length to shaded leaves, which increases k P in shaded leaves relative to sunlit leaves (as assumed in the SPA model, Williams et al. ,1996), or more conductive tissues connecting sunlit leaves to the roots, which increases k P in sunlit leaves (Lemoine et al., 2002).It is advantageous for plants to increase k P to more productive parts of the crown (Katul et al., 2003), but it is as yet unclear how k P actually varies within crowns.Recently, Peltoniemi et al. (2012) argued that it is optimal for plants to distribute total hydraulic conductance in a way that k P is proportional to average PAR for any leaf in the canopy.If this hypothesis is confirmed with measurements, it greatly simplifies a difficult problem, and it can be readily implemented in MAESPA.

Conclusions
We have implemented a new tool to study single-plant or canopy-scale interactions between environmental drivers, canopy structure, weather, and soil water balance.The usefulness of a single-tree model has already been demonstrated by the broad user-base of the MAESTRA model.Here, we have widened the applicability by introducing detailed water balance components, and hydraulic constraints on water use and CO 2 uptake.The new model incorporates a finer level of mechanistic detail than simplified water balance models (e.g.Granier et al., 1999), while still being relatively straightforward to parameterize (see Table B1 for a list of required parameters for the water balance component).
We suggest a way forward in integrating diverse experimental results, by evaluating experimental outcomes in a quantitative framework that summarizes our understanding of the soil-plant-atmosphere continuum.We showed that even relatively straightforward interactions like the C a × drought interaction can be highly variable, because they are dependent on feedbacks of plant size on the soil water balance and acclimation of plant properties due to long-term growth at C a .Quantitative evaluation of the role of such feedbacks is essential if we are to advance our understanding of plant responses to environmental change.Too often in the current literature on C a experiments, responses are simply presented as significant/not significant, rather than being compared quantitatively to expectations based on current theory.As argued by Phillips and Milo (2009), we need to move from asking "Was there a change?" to asking "How large was the change, and is that what we expected?"Table A1.List of symbols, their definition and units.The boundary layer conductance to heat transfer (mol m −2 s −1 )

Appendix B
Below we have summarized the list of parameters needed to simulate the various components of the soil water balance and the response of plant water use to soil water deficit.Parameters that have very little influence or are likely to be near constant are not listed, nor are variables that are derived from the inputs and other constants.An example is fine root density, which has little influence on water uptake apart from the relative distribution of roots in different layers.

Fig. 2 .
Fig. 2. Flowchart of the water balance components of MAESPA, which are taken from the SPA model.The soil compartment is horizontally homogenous, and vertically divided into an arbitrary number of layers.

Fig. 3 .
Fig. 3. Example of a stand of trees as represented in the MAESPA model.Water use and carbon uptake are calculated for a sample of target trees, here shown in red, and then added to give the totals for the stand.It is recommended to select a set of target trees that are representative of the stand, and to account for edge effects. 7

Fig. 5 .
Fig.5.Illustration of the influence of acclimation on the interaction between C a and drought.Panel (a), relative C a effect on total canopy CO 2 uptake (A T ) as a function of soil water content (θ) during a simulated dry-down.Panel (b) shows total canopy transpiration rate (E T ) for the same simulations.The solid line shows the interaction when all plant parameters are unchanged due to growth at eC a ("no acclimation").Note that in very dry soil, transpiration rates are similar in aC a and eC a , but the stimulation of A T due to eC a is much higher than in wet soil (a positive interaction of drought and eC a ).These effects are much reduced for a higher sento L in the stomatal conductance model (" f +0.4 MPa"), which leads to a negative interaction between C a and drought (the C a effect is less in dry soil than wet soil, for both E T and A T ).A lower plant hydraulic conductance (k P ) reduces E T in wet soils, and greatly reduces the positive interaction of C a and drought that was found without acclimation.

Fig. 6 .
Fig.6.Application of MAESPA to an experiment on droughted cherry seedlings byCentritto et al. (1999a, b).Model parameters were based on reported values in the original study, or calibrated to yield a satisfactory fit to the data of the aC a treatment only (see text and Table2).The optimized parameter set was then used to predict water balance in the eC a treatment, and taking into account the observed increase in leaf area in the eC a seedlings.(a) Decline in average soil water content over the rooting zone (θ) as the drought progressed.(b) Decline in the midday leaf water potential for sunlit foliage ( L ).Note a relatively poor fit for the eC a treatment: L is over-predicted early in the drought, and under-predicted towards the end of the drought.(c) The ratio of total water use in eC a to aC a during the dry-down.Higher leaf area in eC a initially leads to higher water use, but this leads to lower θ (a) and L (b), so that eC a seedlings were more water-stressed toward the end of the drydown than their aC a counterparts.The solid line shows the simulations where eC a and the increase in leaf area in the eC a treatment were taken into account, dashed lines show either the direct C a effect only (without leaf area feedback) or the leaf area feedback only (no eC a treatment).
Fig. 7.Effect of the drought treatment on simulated total photosynthesis (summed over the entire 47 day simulation) for the dry-down in the cherry experiment (Fig.6).The aC a simulation shows a 22 % reduction in A T in the dry-down treatment.For eC a , three simulations are shown.Run 1: full simulation but without the leaf area feedback, Run 2: simulation at aC a , but with reduced leaf area to match total water use in the eC a simulation (Run 1), Run 3: including the observed leaf area increase in the eC a treatment.
Model inputs and constantsJ maxMaximum rate of electron transport (µmol m −2 s −1 ) V cmax Maximum rate of Rubisco activity (µmol m −2 s −1 ) I RFraction inhibition of R d in the light (-) C a Atmospheric CO 2 concentration (µmol mol −1 ) γ Shape parameter of the light response of electron transport (-) g 0 g s when A n is zero (residual stomatal conductance) (mol m −2 s −1 ) g 1 Slope parameter of a Ball-Berry-type model of g s (units depend on units of f (D)) min Minimum leaf water potential (Ball-Berry model) (MPa) f L where f L = 0.5 (Tuzet model) (MPa) s f Steepness of the f L function (Tuzet model) (-) k P Plant component of the leaf-specific hydraulic conductance (mmol m −2 s −1 MPa −1 ) u z Above-canopy wind speed (m s −1 ) z H Measurement height of wind speed (m) Soil depth to the bottom of layer i (m) Z Total soil depth (m) β Root distribution parameter (-) Rmin Minimum root water potential (no uptake below this value) (MPa) ω s Tortuosity of the soil air space (-) L d,min Minimum thickness of the dry soil surface layer (m) L v Fine root density (varies by soil layer) (m m −3 ) r 1 Fraction throughfall of rain through canopy (-) r 2 ,r 3 Canopy drainage parameters (mm and -) θ sat Soil porosity (water content at saturation) (varies by layer) (m 3 m −3 ) K sat Saturated soil hydraulic conductivity (varies by layer) (mol m −1 s −1 MPa −1 ) e Parameter for the soil water retention curve (MPa) b Parameter for the soil water retention curve (-) r r Mean fine root radius (m) L T Canopy leaf area index (m 2 m −2 ) T air Air temperature ( • C) Model variables, constants, and outputsA cRubisco-activity limited gross leaf photosynthesis rate (µmol m −2 s −1 ) A jRuBP-regeneration limited gross leaf photosynthesis rate (µmol m −2 s −1 ) A nNet photosynthetic CO 2 uptake rate (µmol m −2 s −1 ) A T Total canopy CO 2 uptake rate (µmol m −2 (ground) s −1 ) R d Dark respiration (µmol m −2 s −1 ) C i Intercellular CO 2 concentration (µmol mol −1 ) C s CO 2 concentration at the leaf surface (µmol mol −1 ) Q Photosynthetic photon flux density at the leaf level (µmol m −2 s −1 ) JElectron transport rate (µmol m −2 s −1 ) * CO 2 -compensation point in the absence of dark respiration (µmol mol −1 ) g s , g C Stomatal conductance to H 2 O (g s ) and CO 2 (g c ) (mol m −2 s −1 ) -air vapour pressure deficit (kPa)k L Leaf-specific hydraulic conductance (mmol m −2 s −1 MPa −1 ) λLatent heat of water vapour (J mol −1 ) sThe slope of the relation between saturation vapour pressure and temperature (PaK −1 ) R n Net radiation (W m −2 ) g B Boundary layer conductance (mol m −2 s −1 ) g V Total conductanceto water vapour (mol m −2 s −1 ) s i Soil water storage in layer i (mm) I i Infiltration into layer i (mm) D i Drainage out of layer i (mm) E i Root water uptake (canopy transpiration) out of layer i (mm) E s Soil surface evaporation (mm) E w Canopy wet evaporation rate (mm timestep −1 ) transpiration rate (mmol m −2 s −1 ) E T Canopy transpiration (for use in water balance) (mol m −2 (ground) s −1 ) R sr Soil-to-root surface hydraulic resistance (MPa s m 2 mol −1 ) G s,t Total conductance to water vapour from soil air space to air above the soil boundary layer (m s −1 ) e s Partial water vapour pressure of the soil pore space (kPa) e a Partial water vapour pressure of the air above the soilk boundary layer (kPa) G ws Conductance of water vapour through the soil pore space (m s −1 ) D eff Effective diffusivity of the soil pore space (m 2 s −1 ) L d Thickness of the dry layer at the soil surface (m) D w Diffusivity of water vapour (function of soil temperature) (m 2 s −1 ) k V von Karman's constant (= 0.41) (-) T i Soil temperature of layer i ( • C) W can Water storage of the canopy (mm) θ i Soil volumetric water content of layer i (m 3 m −3 ) K s Soil hydraulic conductivity (mol m −1 s −1 MPa −1 ) r s Mean distance between roots (m) Q n Soil surface net radiation (W m −2 ) Q e Soil latent heat flux (W m −2 ) Q h Soil sensible heat flux (W m −2 ) Q c Soil heat transport (W m −2 ) K th Soil thermal conductivity (W m −1 K −1 ) ε Soil surface emissivity (-) σ Stefan-Boltzmann constant (W m −2 K −4 ) c pHeat capacity of air (constant) (J kg −1 K −1 ) ρ Density of air (kg m −3 ) G a,m the leaf-specific hydraulic conductance (mmol m −2 s −1 MPa −1 ) min Minimum leaf water potential (only when using Ball-Berry model) (MPa) f L where f L = 0.5 (only when using Tuzet model) (MPa) s f Steepness of the f L function (only when using Tuzet model) (-) Z Total soil depth (m) β Root distribution parameter (-) (or layer-wise specification of rooting density) r 1 Fraction throughfall of rain through canopy (-) r 2 , r 3 Canopy drainage parameters (mm and -) θ sat Soil porosity (water content at saturation) (varies by layer) (m 3 m −3 ) K sat Saturated soil hydraulic conductivity (varies by layer) (mol m −1 s −1 MPa −1 ) e Parameter for the soil water retention curve (MPa) b Parameter for the soil water retention curve (-) ω s Tortuosity of the soil air space (-) (only needed for soil evaporation) L d,min Minimum thickness of the dry soil surface layer (m) (only needed for soil evaporation)

Table 1 .
Summary of origin of the various components of the MAESPA model.For the many references for the components of the MAESTRA model, see Sects.2.1.1 and 2.1.2,and the MAESTRA website.For the components taken from SPA, see

Table 2 .
Cosby et al., 1984)s for the simulation of a dry-down on daily whole-plant gas exchange and plant water relations.The soil type was a loamy sand (parameters fromCosby et al., 1984).Halfhourly weather data were estimated from daily T air amplitude.

Table 3 .
Important parameter settings for the simulation of the C a × drought interaction in cherry.

Table B1 .
List of parameters needed to simulate the water balance.