Application of a Hybrid Forest Growth Model to Evaluate Climate Change Impacts on Productivity, Nutrient Cycling and Mortality in a Montane Forest Ecosystem

Climate change introduces considerable uncertainty in forest management planning and outcomes, potentially undermining efforts at achieving sustainable practices. Here, we describe the development and application of the FORECAST Climate model. Constructed using a hybrid simulation approach, the model includes an explicit representation of the effect of temperature and moisture availability on tree growth and survival, litter decomposition, and nutrient cycling. The model also includes a representation of the impact of increasing atmospheric CO2 on water use efficiency, but no direct CO2 fertilization effect. FORECAST Climate was evaluated for its ability to reproduce the effects of historical climate on Douglas-fir and lodgepole pine growth in a montane forest in southern British Columbia, Canada, as measured using tree ring analysis. The model was subsequently used to project the long-term impacts of alternative future climate change scenarios on forest productivity in young and established stands. There was a close association between predicted sapwood production and measured tree ring chronologies, providing confidence that model is able to predict the relative impact of annual climate variability on tree productivity. Simulations of future climate change suggest a modest increase in productivity in young stands of both species related to an increase in growing season length. In contrast, results showed a negative impact on stemwood biomass production (particularly in the case of lodgepole pine) for established stands due to increased moisture stress mortality.


ForWaDy Description
The ForWaDy (Forest Water Dynamics) model was constructed around the foundation of FORHYM, a general water balance model designed by [2] for the simulation of water fluxes through temperate forest ecosystems. It includes a representation of the vertical flow of water through canopy and soil layer compartments ( Fig. S1). Movement of water through each soil layer is regulated by its physical properties that dictate moisture holding capacity, permanent wilting point moisture content, and infiltration rate. Subsequently, a passive competition for available soil moisture is simulated through the use of an algorithm that combines species-specific root occupancy information with energy-limited transpiration and evaporation demands. Canopy water stress is determined as a function of energy-limited canopy transpiration demand and soil-limited actual canopy transpiration through the calculation of a cumulative water stress index. This index represents a dynamic measure of tree water stress over a given time period and may be used as a parameter to limit tree growth rates based on light and nutrient availability in models such as FORECAST and FORCEE. The simulation of snowfall and snowpack dynamics in ForWaDy is based on the RHYSSys Snow Model [4].
All equations describing water flow within the model are solved using a simulation dt of 0.25 days. The use of a dt less than one day enables the model to divide the flux of water from a particular soil reservoir among competing outflows.

S1.1. Data Requirements
One of the goals in model development was to produce a model that is portable and suitable for use in forest management; thus, climate and site-specific soil and vegetation data requirements were kept to a minimum. The use of parameters that must be calibrated for each site was also avoided when possible. Data requirements are shown in Table S1. species) 1 Solar radiation may be estimated from max and min air temperature, elevation, latitude, slope and aspect using published radiation models. 2 Only annual data required.

S1.2 Radiation interception
Radiation interception in ForWaDy may be estimated from a host model as is the case with FORECAST Climate or it may be done using the default equations within ForWaDy as follows. The energy available for driving evapotranspiration is estimated separately for the canopy, understory, and forest floor based upon the proportion of adjusted total daily solar radiation intercepted by each layer. The interception and extinction of radiation as it passes through the canopy and subsequently through the understory is calculated as a function of LAI using a simple formulation of Beer's Law with an extinction coefficient of 0.4 ( Figure S2). Forest floor radiation is calculated as the total remaining radiation following simulated canopy and understory interception.
For simplification purposes, the transmittance of net radiation (Sn) is assumed to be equal to that of shortwave radiation [10] following the reflection of a portion of total incident shortwave radiation according to a surface albedo (a) (Eq. S1). The default albedo (a) for the canopy and understory is set at 0.12 which is representative of typical values reported for forests [18,9] and should be used as a reasonable approximation in the absence of a measured value. The forest floor albedo is determined as a function of solar zenith angle and moisture content according to [22].

S1.3 Energy-limited evapotranspiration
where α is an experimentally determined coefficient, s, γ, and L are respectively, the slope of the saturation vapor pressure curve, the psychrometric constant, and the latent heat of vaporization of water each evaluated at the daily average air temperature and Sn, G and M represent the daily (24-hour) values of net radiation, soil heat flux and energy storage in the canopy. The daily value of M is ignored in this model as it generally represents less than 5% of Sn in coniferous forests [13]. Thus, using the Sn values calculated for each layer according to the radiation interception submodel, PET is calculated separately for the canopy, understory and forest floor using equations (S3-S5), respectively.  Typically, the use of (Eq. S2) has been limited by the fact that, for dry canopy conditions, α must be calibrated against actual evapotranspiration rates for a particular site using Bowen ratio / energy balance techniques [19] which can be prohibitively expensive. For a range of dry forest conditions, α has been shown to vary from 0.6 to 1.1 [7]. However, for a wide range of relatively smooth, freely evaporating (wet) surfaces α is generally equal to 1.26 [16,20,5]. Values of α less than 1.26, typical in most forests with dry canopies, are thought to be the result of surface control of evaporation through stomatal resistance [11]. In an effort to eliminate the need for α to be experimentally determined for each new site, we suggest the use of a canopy resistance term (RCan) based on mean literature values of α for dry canopies of various forest types (Table S2). Using this method α is set at 1.26 for wet canopy conditions and adjusted under dry canopy conditions (where evaporated water must pass through stomata) through the use of RCan. The use of an approximated RCan will undoubtedly lead to some error in the calculation of canopy PET, but it should provide a reasonable estimate suitable for the intended use of the model.

S1.4 Canopy water, throughfall and canopy evaporation
The amount of water held in the canopy at time t (CWat(t)) is calculated as the difference between incoming rainfall (P) and the sum of throughfall (TF) and canopy evaporation (ECan) evaluated at each time step (dt) (Eq. S7). The fraction of rainfall intercepted by the canopy is determined as a function of canopy vegetation area index (VAI), where a maximum canopy storage term is calculated based on the assumption that at saturation the upper surface of leaves and branches can hold a film of water 0.2 mm in depth (Rutter, 1975), thus CMax = 0.2 * VAI.

S1.5 Soil water storage and drainage
Hydrologic dynamics in the forest floor and rooting zone are simulated using a multi-layered approach in which inflows and outflows are estimated sequentially for each soil layer ( Figure S1). Water storage in and vertical movement through the mineral soil layers are simulated using a "tipping bucket" type algorithm based on total porosity adjusted for coarse fragment content, field capacity and permanent wilting point boundaries (inputs to the model). Water stored in the humus and soil layers between field capacity and permanent wilting point boundaries is considered to be available for plant uptake. Water storage in the litter layer is calculated as a function of fine litter mass (i.e. not including coarse woody debris) per area and is assumed to be unavailable for plant uptake. Lateral flow or interflow of water from soil layers into stream channels is only allowed to occur when the water content of a given soil layer is greater than its estimated field capacity and is based on a coefficient related to slope.

S1.6 Soil-limited evapotranspiration
As the forest floor and soil layers begin to dry, evapotranspiration rates are limited by the availability of soil water. Transpiration loss is calculated separately from each layer for both the canopy and understory. Rooting depths of canopy and understory species are used to determine the layers from which water may be extracted. In order to represent the limited capacity of a drying soil to supply water to roots, particularly during periods of high PET, the daily, energy-limited transpiration rate is scaled back through the use of a relative transpiration rate (RTR) term. RTR is estimated as an empirical function of the percent of available water in each layer and PET [6] (Fig. A3). A root fraction term is also introduced for both the canopy (RFCan) and understory (RFUS) as a scalar to account for the fraction of each layer occupied by roots.
This term is calculated from root depth occupancy data as well as soil layer depths provided by the user. The daily extraction of water from the various layers via transpiration occurs sequentially from the top downward, beginning with the humus layer and continuing down through the soil B layer. The total amount of energy available for canopy transpiration (CanT(total) ) is estimated after accounting for the canopy PET energy consumed by canopy evaporation and canopy resistance (Eqs. S10-S11). CanT(total) is then used to drive transpiration loss from each of the soil layers using an energy budget approach, moving from the humus layer downward as the amount of extractable soil water is depleted in each soil layer (Eqs. S12-S14). The daily extractable fraction of available water in a given layer is largely a function of the estimated RTR value and is generally less than the total plant available water in that layer. The transpiration algorithm continues until either all layers have been depleted of extractable water or the energy available for daily transpiration is depleted. Understory transpiration is calculated simultaneously using the same method.
Drying proceeds from the top downwards (Eqs. S15 & S16) and evaporation from the humus layer is only allowed to occur when the litter layer has dried out.

S1.7 Canopy water stress
Typically, models of forest growth which include water stress as a feedback to forest growth have quantified tree water stress using a summed measure of soil water deficit. One of the problems in a using a summed water deficit for a measure of water stress is that it frequently fails to capture the dynamic interactions between different components of the hydrological cycle and net effect of such interactions on tree water stress. For example, as part of the Biology of Forest Growth Study near Canberra, Australia, [12] demonstrated that cumulative soil water deficit may not be well correlated with tree water stress as defined by pre-dawn xylem potential, a commonly used physiological indicator of tree water stress. The weakness of the technique is even greater when monthly summaries of climate data are used to calculate the water deficit.
In order to evaluate tree water stress more dynamically we propose the use of a cumulative index termed the transpiration deficit index (TDI) (Eq. S17-S18). TDI provides a means of summarizing the net effect of several factors including canopy evaporation, understory transpiration and surface evaporation on actual tree transpiration.
Furthermore the index is calculated on a daily timestep to better capture the cumulative effect of short term water stress events.
CanT (actual) = CanT (humus) + CanT (soilA) + CanT (soilB) (S18) In (Eq. S17) the difference between total canopy transpiration demand (CanT(total)) and actual canopy transpiration (CanT(actual)) is divided by CanT(total) to normalize the effect of leaf area index on the relative canopy transpiration deficit. This adjustment effectively converts the total canopy water deficit to a per unit leaf area measure of water deficit which should be more representative of the water stress of individual trees.