Appendix 1. Key physiological parameters in the hybrid model
Symbol / Parameter / Value / Units / SourceBDBRST / Budburst / 578 / degree-days / Thomson and Moncrieff (1982)
CSTRANS / Max clear-sky transmissivity / 0.74 / - / Bristow and Campbell (1985)
CWSC.LAI / Canopy water storage capacity (LAI) / 0.03 / cm m-2 / Keim (2004)
CWSC.BAI / Canopy water storage capacity (BAI) / 0.025 / cm m-2 / Keim (2004)
GMAX / Max stomatal conductance / 150 / μmol m-2 s-1 / Bond and Kavanagh (1999)
LTCC / Conversion of leaf stomatal conductance
to canopy conductance / 0.13 / - / Tan and Black (1976)
LWP.THOLD / threshold leaf water potential / -2.1 / - / Bond and Kavanagh (1999)
KL / leaf-specific hydraulic conductance / 0.11 / μmol m-2 s-1 / Bond and Kavanagh (1999)
Kb / direct radiation extinction coefficient / function of stand density, LAI, and solar angle / - / Smith (1993)
Kd / diffuse radiation / function of LAI / - / Campbell and Norman (1998)
Vcmax.25 / maximum Rubisco activity at 25°C / 50.6 / μmol m-2 s-1 / Manteret al. (2005)
Jmax.25 / electron transport capacity at 25°C / 122.8 / μmol m-2 s-1 / Manteret al. (2005)
Rdark / dark respiration rate / 1.78 / μmol m-2 s-1 / Manteret al. (2005)
Max.Ci / maximum ratio of intercellular CO2 to ambient CO2 / 0.66 / - / Ripulloneet al. (2003)
Ci.a / parameter describing the relationship between intercellular CO2
and stomatal conductance / 20 / μmol m-2s-1 / Katulet al.(2000)
Ci.b / parameter describing the relationship between intercellular CO2
and stomatal conductance / 0.1 / μmol m-2 s-1 / Katulet al.(2000)
G.critical / critical stomatal conductance / 50 / μmol m-2 s-1 / Fuchs and Livingston (1996)
ANET.2 / net photosynthesis reduction for foliage age class 2 / 0.72 / - / Woodman (1971)
ANET.3 / net photosynthesis reduction for foliage age class 3 / 0.50 / - / Woodman (1971)
ANET.4 / net photosynthesis reduction for foliage age class 4 / 0.30 / - / Woodman (1971)
ANET.5 / net photosynthesis reduction for foliage age class 5 / 0.12 / - / Woodman (1971)
BGC.fol / foliage maintenance respiration rate in FOREST-BGC / 0.002 / kg kg-1 / Korolet al. (1996)
BGC.sap / sapwood maintenance respiration rate in FOREST-BGC / 0.001 / kg kg-1 / Korolet al. (1996)
BGC.brnch / branch maintenance respiration rate in FOREST-BGC / 0.001 / kg kg-1 / Korolet al. (1996)
BGC.root / root maintenance respiration rate in FOREST-BGC / 0.0002 / kg kg-1 / Korolet al. (1996)
CENW.resp / daily respiration rate in CenW / 0.35 / kg kg-1 N / Kirschbaum (1999)
CENW.tr / time constant for the respiratory acclimation response in CenW / 5 / - / Kirschbaum (1999)
Rfrac / fraction of GPP that is lost to respiration / 0.47 / - / Waringet al(1998)
WD / Wood density / 500 / kg m-3
FAC / foliage absorption coefficient / 0.8 / - / Campbell and Norman (1998)
CRC / canopy reflection coefficient / 0.06 / - / Campbell and Norman (1998)
SRC / soil reflection coefficient / 0.04 / - / Campbell and Norman (1998)
PRPC / precipitation interception coefficient / 0.5 / - / Paulet al.(2003)
Appendix 2: Detailed description of the hybrid model
Branch, Crown, And Canopy Simulator (BCACS)
BCACS first estimates the expected number of years to breast height based on the user supplied site index value and equation [6] from Nigh and Mitchell (2003). The height of each tree at the breast-height age 50 is predicted by the dominant height growth equation of Bruce (1981). Height of each annual whorl is then estimated based on the reformulation of the Bruce (1981) dominant height growth equation as presented in Nigh and Mitchell (2003). The total number of live + dead branches, number of live branches, number of whorl branches, and number ofsouth-facing interwhorl branches for each annual stem segment were predicted from segmentage, length, height, and depth into crown and tree crown length (equations [1a]-[1d] in Weiskittelet al.(2007b). The relative heights of branches within the annual segment are assignedwith the empirical distribution presented by Maguire et al. (1994), with whorl branches located at the tip of the annual segment. The first whorl branch is given a random azimuth and the remaining whorl branches are distributed uniformly around the stem circumference(e.g. Doruska and Burkhart, 1994). Interwhorl branches are given a random azimuth consistent with their designation as north- or south-facing. Maximum branch diameter for the annual segment is estimated from the height and relative height of the associated whorl, as well as tree DBH, HT, crown length, and height to crown midpoint (equation [2] in Weiskittel et al. 2007b). Relative branch diameters within an annual segment area function of segment age, relative depth of the segment in the stand canopy, and relative depth of the branch in the annual shoot(equation [2] in Maguireet al. (1994)). Absolute branch diameter is the product of its relative diameter and the maximum diameter for the annual segment. Individual branch angle and length (total and foliated) are predicted from branch diameter, branch height, relative branch height and depth into crown; tree DBH, HT, crown length and height to crown midpoint; and site index (equations [3], [4], and [5] in Weiskittelet al.(2007b)). Branch foliage (eqn. [1]) and woody biomass (eqn. [2]), foliage age class structure (eqn. [3]), and specific leaf area (eqn. [4]) are estimated using equations built from the data presented in Weiskittelet al.(2006; 2008).
[1]FOLBRNCH = 11.1824*BD1.9444*DINC0.3069*BHT-1.1964*(BHT/HT)2.4636*FOLRET0.1539
[2]WODBRNCH = 0.2952*BD2.4449*BHT-0.2907
[3]
[4]
where FOLBRNCH is foliage branch biomass (g), WODBRNCH is branchwood mass (g), BD is branch diameter (mm), BHT is branch height above ground (m), HT is total tree height (m), FOLRET is plot mean foliage retention, CFi is the cumulative proportion of foliage in age class i, DIN_CAN is depth into the canopy (m; maximum tree height in the stand – BHT), MAX.LAI is an indicator variable for maximum leaf area (1 if max leaf area, 0 otherwise), CL is tree crown length (m), SLA is specific leaf area (cm2 g-1), and FOLi is an indicator variable for foliage age class i. Branch surface area is predicted at the tree level from DBH, crown length, crown ratio, and sapwood area (equation [4]in Weiskittel and Maguire (2006a)).
Branch-level estimates of foliage biomass and leaf area by age class are summed to the individual tree- and plot-levels. Estimates of aboveground biomass at the plot-level were used to estimate belowground biomass with the following equation based on 38 published case studies from Douglas-fir (Ranger and Gelhaye, 2001):
[5]BBstand = 0.208*ABstand -1.405
where BBstand is stand belowground biomass (kg ha-1) and ABstand is stand aboveground biomass (kg ha-1). Plot percent canopy cover was estimated by determining maximum crown width of each tree implied by branch length and angle, assuming a circular crown projection, summing crown projections, and correcting for crown overlap (Crookston and Stage, 1999):
[6]CAN_COV =
where MAX.TCRD is maximum tree crown radius (m) and EXPF is tree expansion factor.
Net primary production (NPP)
The required weather stream for driving DF.HGS includes daily records of: maximum air temperature, minimum air temperature, mean air temperature, total precipitation, vapor pressure deficit (VPD), shortwave solar radiation (SW.RADday), and daylength (DAYLEN). Shoot and needle extensionand LAI increase were assumed to occur linearly for 30 days after bud burst (BDBRST), defined as the date that degree-days (>2.78°C) exceed 578 (Thomson and Moncrieff, 1982).
Solar declination, solar zenith angle (SZA), solar azimuth (SOLAZM), solar noon, and half daylength are predicted as standard functions of latitude, longitude, and Julian date (pages 115 – 118 in Gates (1980)). Potential shortwave radiation (PSR) is a standard function of latitude, solar declination, and daylength (Gates 1980). Incoming shortwave radiation was separated into direct and diffuse with the equation of Bristow and Campbell (1985):
[7]
where %RADDIR is the percent direct radiation.Instantaneous direct and diffuse shortwave radiation are calculated using a cosine-approach developed by Wanget al. (2002):
[8]
[9]
[10]
[11]
where SW.RADDIR is instantaneous shortwave direct radiation (W m-2), SW.RADDIF is instantaneous shortwave diffuse radiation (W m-2), SW.RAD.NOONDIR is shortwave direct radiation at noon (W m-2), SW.RAD.NOONDIF is shortwave diffuse radiation at noon (W m-2), SW.RAD.DAYDIR is the daily total shortwave direct radiation (= %RADDIR∙ SW.RADday; W m-2), SW.RAD.DAYDIF is the total shortwave diffuse radiation (= [1-%RADDIR] ∙ SW.RADday; W m-2), SZA is solar zenith angle estimated from latitude, solar declination and hour angle (Gates, 1980), and SZAnoon is solar zenith angle at noon. Shortwave radiation (W m-2) was converted to photosynthetically active radiation (μmol m-2 s-1) by multiplying by 4.55 (PARdir=4.55 SW.RADDIR and PARdif=4.55 SW.RADDIF).
The effective radiation canopy extinction coefficient was estimated from an equation in Smith (1993) and converted to the direct radiation extinction coefficient using an assumption from Goudriaan (1988):
[12]
where kdir is the direction radiation extinction coefficient, RDDF is relative stand density (Curtis, 1982), and all other variables have been defined above. The assumption from Goudriaan (1988) is that the total extinction coefficient is equal to the direct radiation extinction coefficient multiplied by the square root of the of one minus the canopy reflection coefficient. The diffuse radiation extinction coefficient was calculated from a linear relationship with LAI based on a spherical distribution of foliage (Campbell and Norman, 1998):
[13]
where kdif is the diffuse radiation extinction coefficient. Canopy absorbance was calculated separately for direct and diffuse radiation and corrected for incomplete canopy closure:
[14]
where Q is total amount of absorbed PAR (μmol m-2 s-1), Qdir is the total amount of absorbed direct PAR (μmol m-2 s-1), Qdif is the total amount of absorbed diffuse PAR (μmol m-2 s-1), PARdir is total direct incoming PAR (μmol m-2 s-1), PARdif is total incoming diffuse PAR (μmol m-2 s-1), Ω is the estimated canopy clumpiness factor, andthe proportion of radiation reflected by canopy and soil were assumed to be 0.06 and 0.04, respectively. Clumpiness of foliage influences absorption and transmission of solar radiation, so is estimated from the ratio of stand mean crown depth to crown diameter,modified by solar zenith angle (Kucharik et al., 1999):
[14]
where Ω0 is the clumping factor when SZA is 0 (assumed to be 0.7), CL is crown depth (m), and CD is crown diameter (m). The canopy was divided into sunlit and shaded LAI based on the approach of Wanget al. (2002):
[15]
where LAIsun is the amount of leaf area index in full sunlight, LAIshade is the amount of leaf area index in the shade, LAIslope is the horizontal LAI adjusted for percent slope (%SLOPE), SOLAZM is the solar azimuth, ASPECT is site aspect (0° north, 180° south), ASPECTrad is site aspect in radians (north positive, south negative), and kdir is defined above. LAIsun and LAIshade are multiplied by the sunlit and shaded net photosynthesis rate to scale leaf photosynthesis to the canopy.
Water constraints on photosynthesis are determined by both available soil water and vapor pressure deficits. Daily throughfall precipitation and canopy storage of precipitation are determined by stand canopy cover and the surface area of foliage and branches:
[17]
where TF is daily canopy throughfall (mm), CS is maximum canopy storage of precipitation (mm), PRCP is total daily precipitation (mm), LAI is leaf area index (projected), BAI is branch area index (all-sided), and CAN_COV is defined above. Net daily precipitation is then taken as the maximum of two values to account for the difference between precipitation falling on a saturated or unsaturated canopy:
[18] NET_PRCP = max [TF, TF + (PRCP-TF-CS)]
Soil water holding capacity (kg H2O m2) was calculated from soil texture with equations provided inSaxton and Rawls (2006). Canopy transpiration and soil evaporation were estimated on a daily basis from the approach outlined in Watt et al.(2003). Transpiration is calculated using a simple diffusion equation, while evaporation is driven by available energy beneath the canopy and a coupling coefficient. Potential transpiration and evaporation were reduced as soil water decreased or days when rainfall occurred. Plant available water was updated daily by adding precipitation to the plant available water left from the previous day and subtracting evapotranspiration,while not allowing it to exceed the soil water holding capacity. Soil water potential was estimated from soil texture and plant available water with equations from Cosby et al.(1984).
Daily mean leaf water potential was calculated similar to equation [3] in Manteret al (2003), but setting leaf-specific hydraulic conductance to 0.11mmolm−2 s−1MPa−1:
[17]
where ψleaf is leaf water potential (MPa), ψsoil is soil water potential, GMAX is daily maximum canopy conductance (mmol m-2 s-1), and VPD is the vapor pressure deficit (Pa). Stomatal conductance (gs; μmol m-2 s-1) is calculated as:
[18]
where SUNRISEis the local time of sunrise and all other variables have been defined above.
Daylengthwas divided into five normalized Gaussian distances (i.e. 0.0469101, 0.2307534, 0.5, 0.7692465, 0.9530899) and converted to local time, then air temperature was estimated with standard equations of Goudriaan and van Laar (1994):
[19]
where Temp is instantaneous temperature (°C), Tmax is daily maximum temperature (°C), Tmin is daily minimum temperature (°C), TIME is local apparent time, and DAYLEN is day length (s). Vapor pressure deficit (VPD) was then estimated at these same five times from temperature and daily minimum temperature (Tmin).
Net photosynthesis was estimated with the Farquharet al. (1980) biochemical equation. Originally, the necessary parameters for this model were estimated from equations presented in Manteret al. (2005) as a function of foliar nitrogen. The important equations of this calculation were:
[20]
where A1n is the net photosynthesis rate for the 1-yr old foliage age class (μmol m-2 s-1), Ci is the intercellular CO2 pressure (Pa), Ca is the ambient CO2 pressure (Pa), Rdark is the temperature adjusted dark respiration rate (μmol m-2 s-1), V is the temperature adjusted rate of maximum rubisco carboxylation (Vcmax; μmol m-2 s-1), J is the temperature and radiation adjusted value of maximum rate of electron transport (Jmax; μmol m-2 s-1), K is an enzyme kinetics constant (a function of Michaelis constant for CO2 and the inhibition constant for O2), and Г is the partial pressure of O2 (Pa). IntercellularCO2 pressure was conditional on stomatal conductance and ambient CO2 pressure (equation [15] from Katul et al. (2000)).Temperature-modified V and J were estimated withequation [13] and Figure 1 from Manter et al.(2003). Other parameters such as CO2 compensation point in the absence of mitochondrial respiration and the Michaelis-Menten constant of Rubisco are estimated from equations [37] and [39]-[41] from Schwalm and Ek (2004). The dark respiration was estimated at optimal foliar nitrogen using the equation from Figure 5 in Manteret al. (2005) and adjusted for temperature with the equations of Medlyn et al.(1999). The intercellular CO2 pressure was estimated with the approach of Katul et al.(2000)with the critical stomatal conductance set at 50 µmol m-2 s-1(Fuchs and Livingston, 1996) and the maximum ratio of ambient to intercellular CO2 set at 0.66 (Ripulloneet al., 2003):
[21]
The 1-, 2-, 3-, and ≥4-year-old foliage age classes were assumed to be 75, 50, 30, and 10%, respectively, of the net photosynthetic rates of current year foliage calculated from the Farquhar et al. (1980) equation and were inferred from information available in the literature (Woodman, 1971; Manteret al., 2003; Ethieret al., 2006).These net photosynthesis rates are weighted by the amount of leaf area in each age class and multiplied by the corresponding Gaussian weight for each relative time (i.e. 0.1184635, 0.2393144, 0.2844444, 0.2393144, 0.1184635). This mean value was multiplied by the daylength (s), converted to biomass, and summed for an entire year to estimate gross primary production (GPP; kg ha-1). Respiration was initially estimated in several different ways, but none improved the correlation between NPP and observed growth over assuming that net primary production is one-half of GPP.
Annual foliage litterfall (LFALLFOL) is predicted from stand density, age, foliage retention, crown density, slope, and aspect with equation [4] from Weiskittel and Maguire (2006b). Nitrogen mineralization was assumed to occur in the upper 80% of the soil to a maximum of 1 m and was estimated with the approach of Paulet al.(2002):
[22]
where Nmin is daily nitrogen mineralization rate (mg N kg-1), Tempsoil is the soil temperature (°C) estimated with the approach of Paulet al.(2004), PAW is the daily plant available water (kg H2O m-2), WHC is the soil water holding capacity (kg H2O m-2), kn is the optimum nitrogen mineralization rate based on data presented in Chappellet al.(1999), and LFALLFOL is annual foliage litterfall (kg ha-1). To scale this value to the plot-level, soil bulk density is estimated with an equation presented by Kauret al.(2002):
[23]
where BDsoil is the soil bulk density (kg m-3), SOIL_C is the soil carbon content estimated as a function of canopy nitrogen concentration (13.66*CAN_N – 9.99;Perakiset al.(2006)), SOIL_CLAY is the soil clay percentage, and SOIL_SILT is the soil silt percentage.
Daily nitrogen uptake by the trees was estimated from the following function modified from Thornley (1991):
[24]
where Nuptake is nitrogen uptake (kg N day-1), ROOTfine is fine root biomass (kg ha-1), ROOTact is root activity that is dependent on soil temperature, Naval is soil available nitrogen (kg ha-1), ROOT.Cfine is the fine root carbon concentration (assumed to be 0.5), and ROOT.Nfine is the fine root nitrogen concentration (assumed to be equal to CAN_N). Daily Nmin and Nuptake were summed to estimate annual values of each.
ALOGRO
The amount of net primary production (NPP; kg ha-1) allocated to the stem was estimated with the approach of Landsberg and Waring (1997) corrected for bias inherent to the mean tree approach (Duursma and Robinson, 2003):
[25]
where pFS is the foliage to stem partitioning ratio, pNPPstem is the proportion of NPP allocated to the stem, pNPProots is the proportion of NPP allocated to the roots, pNPPwood is the proportion of NPP allocated to the branches and bark, DBHmean is the stand mean diameter at breast height (cm), and DBHstdev is the stand standard deviation of diameter at breast height (cm). The amount of NPP allocated to bark and branches was predicted using the approach of Kirschbaum (1999). The amount of NPP allocated to the roots was estimated using the approach of Landsberg and Waring (1997), except that the soil site fertility variable was estimated as a linear function of canopy foliar nitrogen (e.g. Swenson et al., 2005). NPP allocated to the foliage is taken as the remainder of the NPP not allocated to the other biomass compartments.
Disaggregating stand-level projections of growth to individual trees has been widely debated in the empirical growth and yield model literature (e.g. Qin and Cao, 2006). Process-based models have generally concentrated at the whole-stand level, so little work has been done on mechanistically disaggregating stand NPP to individual trees. Models such as BALANCE (Grote and Pretzsch, 2002) and MAESTRO (Wang, 1990) make highly detailed physiological calculations at the individual-tree level, but a simpler approach was desired for this model. Confronted with the same dilemma, Korol et al. (1995) developed a sophisticated algorithm to disaggregate stand NPP and allocate it to individual trees. Their approach was to divide the canopy into layers of equal leaf area index and within each layer, individual tree radiation absorbance was estimated as a function of tree leaf area (Korol et al., 1995). Allocation to each tree was based on the amount of radiation absorbed by that particular tree as well as an adjustment for hydraulic architecture (Korol et al., 1995).
An alternative approach was also explored as Brunner and Nigh (2003) ran extensive simulations with a highly detailed light interception model to simulate individual tree growth. A regression model was fitted to the data presented in Brunner and Nigh (2003). The final model predicted the mean proportion of incoming light absorbed by an individual crown as a function of its size, relative social position, and stand density. This proportion was then multiplied by the tree leaf area, similar to the weighted leaf areas estimated by Brunner and Nigh (2003). Finally, whole-stand NPP was allocated to each tree based on its proportion of the total stand weighted leaf area: