Changes in the natural dynamics of Nothofagus dombeyi forests: population modeling with increasing drought frequencies

Drought-induced episodes of tree mortality can determine forest dynamics and structure, particularly in forests dominated by single species. Shortand mid-term climate projections indicate that strong changes in annual precipitation may strike more often in northern Patagonia. Data for recruitment, growth, and survival of Nothofagus dombeyi tree individuals were collected at several sites across the Nahuel Huapi National Park in Argentina. We combined mathematically all these different demographic stages into an Integral Projection Model to simulate 100-yr projections of simulated stand structure under different frequencies of extreme drought episodes. We projected total basal area and the number of individuals for three different initial stand types (i.e., young, medium, and old) and for varying drought frequencies (i.e., from 1 to 5 drought events every 100 years). Recruitment into the dbh ≥ 10 cm size class under normal conditions (i.e., without drought) was higher than under episodic drought conditions. In addition, survival under normal conditions was higher than under drought conditions, especially for small trees. Differences in growth were also important, with trees growing more vigorously under normal than under drought conditions. Our simulations predicted that N. dombeyi populations would experience a reduction in tree density in the mid-term if, as predicted by the IPCC projections, the frequency of future drought events increased. The simulations also showed that in those cases, young stands should suffer the most. Drought-mediated changes may induce a decline in the development of N. dombeyi forests in the midand long term by a drastic reduction in tree density.


INTRODUCTION
Climate change is modifying tree growth pattern and stand structure in many forests (Feeley et al. 2007, Reich and Oleksyn 2008, van Mantgem et al. 2009), altering ecosystem services such as carbon balance, water regulation, and nutrient cycling (Breshears andAllen 2002, Anderegg et al. 2012). In fact, climate projections indicate that in some regions, strong changes in annual precipitation will take place (see Figure SPM.7, IPCC 2014), and heat waves will occur more often and will last longer (see Section 2.2, IPCC 2014). This will likely have an impact on the growth dynamics of forests in general. Tree dieoff is already occurring in forests across the world associated with chronic or episodic conditions of strong drought (Allen et al. 2010. In many temperate forests affected by these drought-induced die-offs, dynamics is driven by a single or few dominant species, such as Pinus sylvestris in central and southern Europe (Galiano et al. 2010, Rigling et al. 2013, Pinus edulis (Breshears et al. 2005) and Populus tremuloides (Worrall et al. 2010, Anderegg et al. 2013 in western North America, Nothofagus dombeyi in Andean Patagonia (Suarez et al. 2004, Suarez andKitzberger 2008), and Eucalyptus sp pl. in Australia (Fensham andHolman 1999, Matusick et al. 2013). Thus, there is a need to assess the performance of those tree populations in scenarios in which drought episodes are likely to become more frequent.
Nothofagus dombeyi (Mirb.) Blume is a dominant, evergreen species growing in a wide range of edaphic conditions, from mesic to cool temperate rainforests, in the Patagonian Andes (Chile and Argentina) in South America. At the eastern side of the Andes (Argentina), and as a consequence of the rain shadow effect caused by mountains, N. dombeyi forests develop in sites with annual precipitation ranging from 4500 to 1400 mm/yr. As the precipitation declines toward east, only 80 km far, drought conditions in summer have a stronger effect on the vegetation and can extend for several weeks (Villalba and Veblen 1998, Villalba et al. 2003, Suarez and Kitzberger 2010. In Argentina, the species dynamics is determined by small-scale gap-phase regeneration following tree fall across the climatic gradient, and comprises both monospecific stands and mixed forests with the evergreen conifer Austrocedrus chilensis (Veblen 1989). The absence of strong competition of shade-tolerant species favors N. dombeyi regeneration even in small gaps, and also in the presence of A. chilensis, a more shadetolerant species (Suarez and Kitzberger 2008). The effect of climate variability on vegetation dynamics acts at different scales across the strong precipitation gradient prevailing in northern Patagonia. Thus, in mesic forests, strong but relatively short climatic fluctuations impact forest structure through direct effects on entire tree cohorts (Suarez and Kitzberger 2010). In contrast, in more humid forests, climatic-induced mechanical disturbance (tree fall) prevails, driving mortality and subsequent growth and establishment at individual level (Suarez and Kitzberger 2010). Like other forests types growing in northern Patagonia, tree demography of N. dombeyi forests is also indirectly affected by inter-annual climatic variability via fire regimes , Veblen et al. 1999.
Recently, climate variability has been directly affecting the demography and the dominance of N. dombeyi. Massive die-off and mortality events have occurred in response to extreme droughts, such as the 1998-99 episode in northern Patagonia (Bran et al. 2001). In heavily affected stands, N. dombeyi tree density diminished as much as 57% in adults and 30% in saplings (Suarez et al. 2004, Suarez andKitzberger 2008). In contrast, co-occurring A. chilensis adults and saplings remained mostly alive, favoring A. chilensis dominance in the sapling cohort in at least 10% of stands (Suarez and Kitzberger 2008). In addition, drought-induced gaps force changes in the physical microenvironment, likely reducing recruitment and survival rates of N. dombeyi saplings and promoting the recruitment of more droughtand shade-tolerant species (Suarez and Kitzberger 2008). However, even lower N. dombeyi sapling density may assure gap filling owing to higher growth rates in drought gaps with open shrub layer (Suarez and Sasal 2012).
Although barely explored, previous results suggest that healthy trees may survive to drought by shutting off a substantial amount of their upper canopy. Resilience of N. dombeyi forests to drought episodes would rely on this partial crown survival and on recruits established in resulting gaps (Lloret et al. 2004(Lloret et al. , 2012. Thus, the dynamic scheme in northern Patagonia is the coexistence of both partially killed and dead standing trees following drought events. Nonetheless, predicted changes in the intensity and periodicity of largescale climatic features (IPCC 2014) are expected to modify the frequency and intensity of droughtrelated mortality, the rate of partial crown dieback, and, ultimately, the negative effect of successive crown dieback events on tree resistance. Considering that these recent extreme droughts had the strongest negative impact on N. dombeyi populations at the eastern distributional limits, a change in the distribution of this species as a consequence of the projected climate change should be expected. Therefore, N. dombeyi constitutes an appropriate target to study changes in forest structure associated with drought episodes, likely due to climate change, which may be affecting temperate ecosystems. Moreover, the study of this species can also address the question of whether a higher occurrence of drought events will hinder the capability of forests to recover from those episodes, as determined by their resistance and resilience (Lloret et al. 2011). In spite of the existing information about the demography and growth patterns of this key species in austral forests, we lack an integrative approach with which to evaluate the future fate of these forests in a context of changing climate. This assessment could be obtained by integrating mathematically the different parameters that determine the stand dynamics of N. dombeyi populations into a simulation model.
Integral Projection Models (hereafter, IPM) consist of a set of mathematical and statistical tools to study the population dynamics of plants or animals under different scenarios (see Easterling et al. 2000). Unlike Matrix Population Models (hereafter, MPM), which require the population to be discretized into classes (see Caswell 2001), the IPM methodology allows a structured population to be defined naturally by a continuous variable (e.g., size). Furthermore, most population processes that form part of the IPM algorithm (usually, survival, growth, and recruitment) can be calculated even if only a limited dataset is available, in contrast to the MPM methodology. The properties of IPM are particularly adequate when information is scarce in terms of replication or relative to relevant demographic processes (e.g., recruitment). That is the case of forested areas without systematic forest inventories, such as Andean Patagonian forests. On the other hand, one of the main drawbacks of the IPM methodology is that it needs an integral to be solved instead of the simple matrix multiplication operation of the MPM. Moreover, in most cases that integral cannot be solved analytically, and it has to be determined by means of appropriate numerical quadrature methods.
Plant population studies have traditionally investigated the long-term, asymptotic dynamics in stands or plots. A battery of diagnostics is usually applied to the resulting population stable states and their growth rate. However, short-and mid-term dynamics (see, e.g., Stott et al. 2011) are important on its own. As stated by Caswell (2001:18), "Transient analysis, focusing on shortterm behavior, may be more relevant than asymptotic analysis in characterizing the response of a population to perturbations." Thus, in studies facing the dynamics associated with disturbances, it is quite appropriate to shift our focus from stable to transient states (Fukami and Nakajima 2011). In particular, stable-state solutions may be of little use if the starting assumptions of the simulations cannot be guaranteed to hold in the midto long term.
In this study, the main hypothesis we tried to assess was that resilience of N. dombeyi may be limited by recurrent episodes of drought, as expected in climate change scenarios for the eastern side of the Andean cordillera. Thus, to test this hypothesis, we undertook the following tasks: (1) We collected field data from transects and plots in N. dombeyi forests to determine tree recruitment, growth, and mortality under either normal or episodic drought cases, and (2) we implemented those processes in a IPM to project and evaluate the ensuing mid-term (~100 yr) N. dombeyi forest structure (i.e., size distribution, basal area distribution, and tree density per hectare) in those two cases. The design and implementation of the IPM is described in detail in the Materials and Methods section, whereas the calculated projections for the N. dombeyi stands are shown in the Results section. This outline was chosen for the sake of readability and easier comprehension of the present work, although the particular IPM implementation we present in this manuscript can be considered on its own as part of the results of the study.

Study area, species description, and climatic record
The datasets used in this study come from different surveys, all conducted in the Nahuel Huapi National Park (c. 40°36 0 to 41°32 0 S and 71°16 0 to 40 0 W; hereafter NHNP) in northern Patagonia (Argentina) between the years 2000 and 2008. In NHNP, the vegetation established on the foothills of the Andes is strongly determined by the precipitation gradient imposed by the rain shadow effect of the cordillera. Mean precipitation decreases abruptly from~3000 mm/yr at the main Andean cordillera to <500 mm/yr only 80 km eastward (De Fina 1972). Parallel to this precipitation gradient, vegetation changes abruptly from Valdivian temperate rainforests, which are moist forests dominated by tall evergreen N. dombeyi, to semiarid steppe dominated by cushion shrubs and bunchgrasses. As precipitation decreases to ca. 1500 mm/yr, N. dombeyi and A. chilensis codominate extensive stands. Nearer the ecotone with steppe, the conifer A. chilensis grows as open woodland interspersed with sclerophyllous shrublands. Field sampling was done in pure N. dombeyi and mixed N. dombeyi-A. chilensis forests, which formed even-aged populations following stand-devastating fires that occurred about one century ago. After a stand-devastating fire of natural or anthropic origin, both species establish simultaneously if seed sources are close to the burned area (both species are obligate seeders). However, the higher growth rate of the more shade-tolerant species, N. dombeyi, implies that this broadleaf species tends to dominate exclusively in the main canopy, whereas A. chilensis populations remain partially suppressed in the sub-canopy (Veblen and Lorenz 1987).
Climatic records for Patagonia are scarce and unevenly distributed in time and space. The Bariloche Airport station is the only meteorological station close to the studied N. dombeyi populations (<50 km), providing century-long records, including the 1998-1999 drought period. This station is representative of the regional climatic condition (Villalba et al. 2003, an intense drought occurred as a consequence of a strong La Niña event. Rainfall in Bariloche Airport station was 2.6 standard deviations below the historical average record. Unusual high summer temperatures (mean daily temperatures above 7°C and maximum daily temperatures above 28°C during a continuous 19-day period in January 1999) concurred with low precipitation (see Appendix S2: Fig. S1 for total annual precipitation data). As a result, more than 11,000 ha of forests (~10% of the N. dombeyi forests) growing in NHNP exhibited die-off (Bran et al. 2001). For a detailed specification or particular description of sampling, sites, and sampled plots, see Suarez (2001Suarez ( , 2010, Suarez et al. (2004), Kitzberger (2008, 2010), and Suarez and Sasal (2012).

Data acquisition
Transect data.-In 2004-2006, six transects were established in six areas (i.e., one transect per area) of N. dombeyi-dominated forests (Fig. 1). Those six localities covered different precipitation regimes, from three pure monospecific N. dombeyi forests growing under 2500-3000 mm/yr of rain, to three mixed N. dombeyi-A. chilensis forests with precipitation of 1200 mm/yr. Transect survey consisted of choosing sampling points (approximately every 50-100 m) along each transect. Transects were established in representative, homogeneous areas of well-developed and undisturbed (last fire happened >150 years ago) stands dominated by N. dombeyi. In each sampling point, increment cores (as low as possible, <30 cm above the ❖ www.esajournals.org ground) from the four live and dead N. dombeyi adult trees closest to the sampling point location were obtained using increment borers. Those samples were then processed using dendrochronological methods (Stokes and Smiley 1968), which allowed for a precise determination of the establishment year, mortality dates, as well as the retrospective analysis of the growth pattern. The date of an annual ring was assigned to the calendar year when growth began (Schulman 1956). Tree-ring widths were measured under a stereomicroscope, using a HENSON measuring bench (0.01-mm resolution). The quality of the measurements was examined using the software COFE-CHA (Holmes 1983). The establishment date corresponded to the calendar date of the innermost ring of the core. Live trees were used to build standard chronologies with the ARSTAN program for each study site (Cook and Krusic 2006). Next, mortality dates were obtained using the chronologies developed for each study site using COFE-CHA as master information. Further details about establishment and mortality date determination and correction factors employed can be found in Suarez and Kitzberger (2010). Diameter at breast height (hereafter dbh) was also recorded during sampling and the corresponding bark thickness was registered following the extraction of the increment core. Data on dbh were obtained for 591 live and 487 dead trees. For growth inspection, 537 live trees could finally be cross-dated. In addition, 485 dead trees were cross-dated to determine their death date. These transect data were then used to estimate growth parameters and mortality/survival rates (see Mortality/survival and Growth for a description of the different methodologies).
Plot data.-In 2000, 26 pairs of mixed N. dombeyi-A. chilensis 20 9 20 m plots were chosen based on a stratified random sampling procedure. Each pair of plots consisted of one "dead" (with >50% of dead trees) and one "live" (with <25% of dead trees) sites located nearby. We thus obtained 26 sites encompassing the entire area affected by mortality. The affected zone was identified using aerial photographs (Bran et al. 2001), and then, a 50 9 50 m square grid was built to locate paired plots. Altitude, aspect, and proximity to rivers were similar among dead and live sites, although dead sites tended to be located in steeper, rockier slopes (Suarez et al. 2004). In each plot, we measured the dbh of live and dead trees with dbh ≥ 10 cm.

Data processing
All statistical analyses were carried out within the R (2014) environment.
Bark width vs. size determination.-Points of observed bark width as a function of dbh (see Appendix S2: Fig. S2) approximately followed a parabola in a log-dbh scale. At the same time, the dispersion in the thickness data seemed to become greater as the mean thickness grew larger. Consequently, we employed a generalized linear algorithm (glm function in R) with a gamma distributional model and an identity link function such as: to fit the observed bark dataset, where x was dbh and a Φ , b Φ , and c Φ were the parameters to be determined. The generalized linear regression fit yielded parameters a Φ and c Φ that were statistically significant (P < 10 À16 ). Parameter b Φ was not significant (P = 0.442), but was anyhow maintained in the regression after a stepwise Akaike model selection scheme (stepAIC function in R) chose the full model as the best one for prediction purposes (Burnham and Anderson 1998). Notice that Eq. 1 was calculated and applied solely in situations where dbh ≥ 10 cm, thus ruling out extrapolations below 10 cm. Initial tree size distribution determination.-We matched the observed distribution of adult tree dbh values from the observed 52 stands with a gamma distribution with shape (j) and scale (h) parameters (see Appendix S2: Fig. S3): where I is a proportionality constant which indicates the total number of individuals per ha in the observed initial distribution. The shape and scale parameters were straightforwardly computed by the method of moments. The main drawback was that no statistical significance level could be provided. Nevertheless, the resulting curve shown in Appendix S2: Fig. S3 displays a remarkably close match to the shape of a histogram of the observed data.
The gamma distribution was chosen due to its flexibility to represent different tree stand types simply by changing the shape parameter. Starting from the observed distribution, decreasing or increasing its shape parameters j would shift the resulting gamma distribution toward smaller or larger diameters, respectively, creating younger or more mature simulated tree distributions. Other distributional options, not explored in the present work, may be, for example, Weibull or beta (Jaworski andPodalski 2012, Ogana et al. 2015).

Integral projection modeling
Since the study area had never been surveyed systematically and periodically, we used ring growth and bark data from field sampling to characterize trees in the plots as they were before the sampling. Ingrowth, survival, and growth functions, as used in the IPM, were then derived from these data, as explained below. Simulations were subsequently carried out under different frequencies of extreme drought episodes that were equivalent in severity to the one observed in 1998-1999. Consequently, no distinction was made between datasets from plots and/or transects in drier or wetter areas.
Methodology.-The collected demographic data enabled us to determine the demographic stages (i.e., recruitment, survival, and growth) used to build up an IPM. In the IPM methodology, changes to the tree population at time t, structured by a continuous state variable x, are subsequently determined by means of an integral equation (e.g., Easterling et al. 2000): (3) where y is, in turn, the value of the continuous variable at a later time t + D. The term pðx; yÞ in the integrand of Eq. 3 is usually written as survival sðxÞ times growth (conditioned on survival) gðx; yÞ: In addition, f ðx; yÞ is identified with a recruitment term (i.e., the generation of new offsprings with continuous variable y generated by parent trees with continuous variable x; in other studies, it is also called fecundity). The integral is then evaluated between appropriate lower and upper limits. In this study, as is often done in tree studies, continuous variables x and y represent the diameter at breast height (dbh) of N. dombeyi individuals at times t and t + D, respectively. That is, we used size, not age, to describe the structure of tree stands at every time step.
Due to the characteristics of our sampling procedure, recruitment could not be determined as a function of size of parent tree. However, it was still possible to outline an alternative strategy to estimate recruitment from the datasets (see Recruitment/ingrowth for a detailed description).
To that aim, we simplified Eq. 3 above by dropping the dependence of f ðx; yÞ on x. We thus rewrote it as follows: (5) where FðyÞ was the new recruitment function, which hereafter we will denote also as ingrowth. Consequently, the total influx or ingrowth of new individuals at each time step was assumed to be only a function of size at time t + D, y, and it did not depend on any characteristics of the stand at time t. This and other simplifying assumptions (like leaving out competition effects) were forced upon us due to field sampling constraints and may arguably work when projecting the model into the near future. However, it will very likely be overly inadequate for very long time projections. Consequently, we limited our calculations to mid-term projections of 100 yr. This period of time harbors the immediate forest dynamics associated with the currently ongoing climate change, as investigated by Suarez and Kitzberger (2010). Results for size distribution and derived quantities were always given in units per ha.
Throughout our study, time interval D was 10 yr. This number was chosen as a trade-off, so that populations would show significant changes (i.e., changes in survival, growth, and recruitment of the populations were large enough to be measured with confidence) between t and t + D. At the same time, however, it was ensured that those changes were not too large that the dynamics of the population could not be appropriately explained by the functional expressions of the selected demographic stages. Similar D values have typically been chosen in national inventories elsewhere (e.g., D = 10 in Spain, Vallejo and Villanueva 2002, D = 10 in Austria, Schadauer et al. 2007, D = 5-12 in France, Vidal et al. 2007).
The integral term on the right-hand side of Eq. 5 above determined the dynamics of adult N. dombeyi individuals. We defined an adult individual to be any tree with dbh ≥ 10 cm, which therefore became the lower limit of the integral (i.e., L = 10) and of the recruitment function, FðyÞ. The lower limit of dbh ≥ 10 cm was imposed upon us by field work procedures. In addition, maximum dbh adopted was 250 cm (i.e., U = 250), which was approximately its maximum observed value. The integral in Eq. 5 was impossible to compute analytically and we were thus compelled to use numerical quadrature procedures to solve it, as described in Appendix S1.
We sought to determine the sðxÞ, gðx; yÞ, and FðyÞ functions in Eq. 5 that characterized, respectively, the survival, growth, and ingrowth of N. dombeyi during a 10-year period. Those functions were in turn computed under either normal or severe drought conditions separately. We did so by choosing a time interval with average rain and temperature values in all observed areas (i.e., 1988-1997), as well as another time interval with severe episodic drought conditions (i.e., from 1998 onward). We assumed that a severe drought would only affect N. dombeyi survival, growth, and recruitment during the very first years (≤10) after that event, in agreement with studies from temperate forests (e.g., Lloret et al. 2011). For example, a 26% growth reduction was detected during the immediate years following the 1998-99 drought event in N. dombeyi adult trees growing at the eastern distributional limit (Suarez 2009). Furthermore, survival of N. dombeyi trees with higher percentage of crown damage was considerably lower during the same period.
We compensated, when needed, for the different lengths of the ring-width series from all stations (ring series ended between 2004 and 2007, and in that case, a complete set of 10 consecutive years, starting in 1998, was often unavailable) by multiplying for a proportional factor to make all functions equivalent to those of a 10-year period.
Recruitment/ingrowth. -From the sampled stands, we selected live trees, as well as trees that were unambiguously known to have died in 1998, with dbh ≥ 10 cm. That dbh was then numerically shrunk with the help of the equations described in section Growth below to reflect their diameter 10 years before sampling. We next assumed that those individuals whose shrunk dbh was smaller than 10 cm could then be counted as 10-yr ingrowth. Moreover, dbh was shrunk differently to reflect normal or severe drought conditions. Assuming that it was likely that the number of ingrowth trees approximately followed a simple exponential expression, such that there were fewer trees at larger sizes, we represented the distribution of ingrowth trees as: where parameters a F and b F were determined as the number of ingrowth trees and their mean dbh, respectively. Mortality/survival.-We could not calculate tree survival as a function of x directly from our datasets. As an alternative methodological route to determine survival, however, we proceeded by first formulating Bayes' theorem. The survival probability function we were seeking, sðxÞ, could be formulated as follows: sðxÞ ¼ 1 À pðD=x; AÞ ¼ 1 À pðx; A=DÞ Á pðDÞ pðx; AÞ Labels A and D denote trees alive 10 years before a given date and dead at that date, respectively. That is, sðxÞ is one minus the probability for a tree, alive and with size x, to die during the next 10 years, pðD=x; AÞ. Term pðx; AÞ is the size distribution of live trees at the beginning of the measurement interval. We adopted the simplifying assumption that pðx; AÞ was identical to the current size distribution of live trees. With this simplification, we could compute pðx; AÞ from the current sampled distribution of live trees. In turn, pðDÞ was simply the current proportion of dead trees from sampled plots. Finally, term pðx; A=DÞ represented the probability for a dead tree at a given date to be alive 10 years before that date and with size x. This term could be calculated as follows: We inspected the tree-ring series of our sample of ring-dated dead trees to select only the fraction that had been alive sometime during the last 10 years. The dbh of these trees was then shrunk with the help of the equations put forward in the following section to represent their diameter 10 years before. Then, we calculated a histogram of the dbh values from ❖ www.esajournals.org those shrunk trees. Finally, those counts per dbh interval in the histogram were scaled to represent proportions with respect to the total sample of current dead trees used in these calculations.
To obtain a functional form for adult tree survival sðxÞ, we logit-transformed those proportions and fitted a straight line by least squares. The fitted parameters were then used to express survival as follows: where a s and b s were determined from the fit.
The 10-yr survival function sðxÞ was calculated under either normal or episodic drought conditions by using inverse growth equations (see Section Growth below) for the 1988-1997 and 1998-onward time intervals, respectively. Growth.
-We defined tree growth as the increment in dbh during a time interval of 10 years. The absence of a complete tree inventory obtained at two different time intervals did not allow a determination of dbh increment directly from consecutive diameter measurements. However, treering data were available for a number of live trees, whose dbh had also been determined at the time of measurement. Consequently, we could discount accumulated ring width for individual ringdated trees, starting from its current dbh data. However, discounting ring widths from dbh had to be done including bark width into the calculation. Given the tree-ring radial decrement dR (determined from tree-ring measurements for each individual tree), plus function UðdbhÞ in Eq. 1, we could simply write for each ring-dated tree: where D t stands for tree diameter without bark at time t. Then, at a past time tÀD, we had: as the new diameter without bark at time tÀD. Finally: or: represented the new dbh at a past time tÀD. The Φ function has a nonlinear dependence on dbh.
Therefore, to solve for dbh tÀD for each individual tree, we had to implement an iterative Newton-Raphson scheme. Diameter increment dbh t Àdbh tÀD was plotted vs. initial tree dbh tÀD . The nlgamlss function in the gamlss.nl R package (Rigby and Stasinopoulos 2005) was then employed to fit a negative exponential function to the response variable with dbh tÀD as predictor, assuming a log-normal distributional model (which conveniently forced diameter increment to be nonnegative, as observed in the field).
These calculations were carried out for either normal or severe drought conditions by using tree-ring data from either the 1988-1997 or the 1998-onward time intervals, respectively. Notice that the derivation above would be exactly equivalent to using dbh t+D Àdbh t (i.e., yÀx) vs. dbh t (i.e., x) instead. Growth function gðx; yÞconsequently became: where l was a negative exponential that depended on size x as follows: and r, a g , and b g were parameters to be determined. So-called inverse growth equations, employed during this study to compute past tree dbh values, were similarly derived by carrying out an analogous fit to the diameter decrement dbh tÀD Àdbh t .
Simulations. -Computing 10-yr survival, growth, and recruitment under either normal (i.e., from 1988 to 1997) or severe drought (i.e., from 1998onward) conditions allowed us to simulate the dynamics of N. dombeyi populations under either circumstance with varying frequency. To that aim, we chose the probability p D with which drought conditions could randomly hit the simulated tree stands in a period of 100 years (in steps of 10 yr) at each simulation. The decision whether to select survival, growth, and ingrowth functions affected by drought was then taken at each time step by drawing a random number z 2 ½0; 1 such that z < p D . For example, choosing p D = 0.1 meant that, on average over all simulations, there would be one 10-year severe drought period every 100 years. After choosing a given p D and a starting ❖ www.esajournals.org stand type (young, medium, or more mature), a simulation run consisted in applying the IPM algorithm iteratively, drawing a number z at each step and subsequently selecting the corresponding (i.e., normal or drought) survival, growth, and ingrowth functions.
Appendix S2: Fig. S4 shows the flowchart of the complete algorithm, for one simulation run. For each p D value, we then ran a total of 100 different simulations, which yielded a representative sample of all the possible trajectories for the N. dombeyi stand dynamics. Tree size distribution, basal area, and total number of trees were calculated at every time step. The main parameters of the IPM implementation are specified in Table 1.
We selected three different j, I value pairs to derive three different initial dbh distributions for the simulations with the help of Eq. 2. The three pairs of parameters were meant to be representative of young, observed (which we assumed to be medium age), and more mature stands. The scale parameter h was not modified and was equal to the observed value for the average stand, ϑ = 9.9. Resulting N. dombeyi size distributions are shown in Fig. 2. The set of values j = 4.6, I = 267.3 reproduced the average observed tree size distribution per ha shown in Appendix S2: Fig. S3. We also simulated a tree stand with smaller dbh values (i.e., a young tree stand) and many more individuals (j = 0.9, I = 534.6). Finally, parameter values were chosen to simulate a tree stand with large dbh values (i.e., a more mature tree stand) and fewer individuals (j = 9.1, I = 80.2). Table 2 lists the model parameter values, computed as described above and used in the IPM algorithm. A few of the parameters turned out not to be statistically significant, but were nevertheless kept in the model, as explained below.

IPM functions
Ingrowth. -Fig. 3 shows the calculated exponential distributions of ingrowth trees that entered the class of adult trees during a period of 10 years. A Mann-Whitney U test showed that the means (i.e., exponent b F ) of the two observed dbh distributions under normal and drought conditions did not have statistically significant differences. Therefore, the two curves only differed in the total number of ingrowth trees a F , which we computed from observational data.
Survival.-The results (Fig. 4) depict the observed (dots) and fitted logistic (lines) 10-yr survival probability curves under normal and severe drought conditions. Both survival curves appear very similar for larger trees. However, for smaller trees, survival becomes lower under drought than under normal conditions.
Growth.-dbh increment (dbh t + 10 Àdbh t ) as a function of dbh t is shown in Fig. 5. These results show that growth is lower under drought conditions. Further, growth also decreases monotonically under both normal and drought conditions as dbh t becomes larger. These facts demonstrate the appropriateness of the negative exponential No. of trees θ = 9.9, κ = 0.9, = 534.6 θ = 9.9, κ = 4.6, = 267.3 θ = 9.9, κ = 9.1, = 80.2 Ι Ι Ι Fig. 2. Initial dbh distributions that were used in this study. Curves were defined by continuous gamma curves with distinctive combinations of parameters as defined in the text, describing three different types of initial tree stands: (a) young stand with smaller trees (dotted line), j = 0.9, I = 534.6, h = 9.9; (b) observed stand (solid line), j = 4.6, I = 267.3, h = 9.9; and (c) older stand with larger trees (dashed line), j = 9.1, I = 80.2, h = 9.9. choice as a model for size increment, as is shown by the resulting fits, drawn as solid lines in Fig. 5.

IPM simulations
The results for the 100-yr simulation runs for the three initial simulated tree stands are shown in Fig. 6. Droughts that are more frequent lead to a general drop in tree numbers, especially at small sizes. Moreover, simulation results are more variable (i.e., the gray area in Fig. 6a-c is wider) when the IPM algorithm starts from a younger stand with smaller individuals.
The temporal dynamics of the basal area results in the 100-yr simulations (Fig. 7) show that area increment is faster for those stands that initially have smaller individuals. Its approximately linear increment is due to the lack of a limiting factor that could account for intra-species competition.
However, even though competition was unaccounted for in our simulations, we can notice that basal area results for p D = 0.5 (lighter gray band in Fig. 7a-c) seem to flatten out, suggesting that frequent severe droughts may take a heavy toll on stand dynamics, compromising their future. In addition, young stands show the largest relative differences between the results for normal and for drought conditions.
As pointed out for the basal area case, the differences displayed in total number of trees for the normal and drought results are largest for young stands (Fig. 8). The effect of higher mortality for lower sizes is conspicuous for young No drought: Eqs. 13 and 14 a g = 1.450 ÃÃÃ b g = 0.005 Ã r g = 0.670 ÃÃÃ Drought: a g = 0.998 ÃÃÃ b g = 0.003 r g = 0.724 ÃÃÃ Notes: The significance levels for r g parameters actually correspond to the significance of its log-transformed value. Equation numbers used in the text are included for reference. Significance levels are given as Ã (P < 0.05), ÃÃ (P < 0.01), and ÃÃÃ (P < 0.001).  Fig. 4. Tree survival as a function of dbh, under normal and drought conditions. Solid and dashed lines are logistic curves fitted to observed normal (empty bullets) and drought (asterisks) data points, respectively.
( Fig. 8a) and observed (Fig. 8b) stands. An increase is detected only in stands that have an initial size distribution composed mainly of large trees (Fig. 8c). In addition, in this latter case, the results for both the basal area (Fig. 8c) and the number of trees (Fig. 8c) level off for highfrequency droughts.

DISCUSSION
Our study reveals that the combined effect of drought episodes on mortality, growth, and recruitment would give rise to noticeable midterm changes in N. dombeyi forests. A great deal of previous research has addressed the impact of drought on tree mortality (e.g., Fensham and Holman 1999, Bigler et al. 2006), growth (e.g., Camarero et al. 2011, Lloret et al. 2011), and recruitment (e.g., Olano and Palmer 2003, Redmond and Barger 2013 separately. However, there is a need for Fig. 5. dbh increment (growth) as a function of tree dbh. Growth under (a) normal conditions and (b) drought conditions was computed. Empty bullets correspond to observed data, whereas solid line represents the log-normal fit (Eq. 14). Notice that the response variable in the fit was the increment in dbh between t and t + 10, rather than dbh itself. Fig. 6. Initial (dashed line) and 100-yr simulations of dbh distributions. Shaded bands illustrate the interquartile range of simulated results for the different 10-yr occurrence probabilities of a severe drought event. Those probabilities are shown in the legends within the plots. The solid line delineates the 100-yr simulation result for a tree stand under normal conditions (i.e., without drought events). Plots, from top to bottom, correspond to different parameter values, chosen to describe distinctive initial dbh distributions: (a) young stand with smaller trees, j = 0.9, I = 534.6, h = 9.9; (b) observed stand, j = 4.6, I = 267.3, h = 9.9; and (c) older stand with larger trees, j = 9.1, I = 80.2, h = 9.9. studies to examine the dynamics of tree stands associated with die-off events in an integrative way. This approach is essential, due to the existence of possible mechanisms of compensation between demographic processes (Lloret et al. 2012). In N. dombeyi forests, drought episodes may separately enhance mortality, particularly in small trees, and reduce growth, whereas its effect on recruitment is arguably more complex. When implemented together in a population dynamics model as the one described above, the fate of N. dombeyi populations appears hampered in the long term by a reduction in tree density, especially if the frequency of future drought events increases as predicted by some models. The relevance of our results is enhanced when one considers that IPCC Regional Climate Projections (2014) indicate that unusual negative precipitation changes may be expected in and around our study area. Therefore, the resilience to a single drought event demonstrated by the modeled N. dombeyi stands may not be enough to cope with a strong increase in drought recurrence.  As described in Williams et al. (2012), IPM implementations may suffer from so-called eviction effects, whereby individuals may simply "leave" the model at a given time step by shrinking or growing beyond the lower or upper dbh limits of the IPM integral, respectively. Not accounting for this eviction effect may bias the population growth rate and other diagnostics. However, there is no universal "cure" for eviction (Rees et al. 2014). By choosing a relatively short projection time of 100 yr, we eliminated eviction of large individuals for all practical purposes. In addition, by choosing the lognormal distribution as our distributional model for dbh growth, we readily eliminated eviction for small trees.

Drought effect on stand type
In the reported N. dombeyi populations, drought-induced mortality was higher in the smaller size classes. This has long-term consequences for the affected stands. Under recurrent droughts, population projections that started from stands dominated by young trees showed a greater reduction in tree density, as opposed to the ones dominated by older trees (Fig. 7). Juvenile and small trees in other species have also been reported to exhibit higher impact by drought events elsewhere (Mart ınez-Vilalta and Piñol 2002, Lloret et al. 2004, Galiano et al. 2010. This is particularly so in N. dombeyi populations (Suarez et al. 2004), which may be attributable to their lower ability to store reserves, or to build root systems to uptake water. In contrast, other studies have found higher mortality in larger trees (Slik 2004, Muller-Landau et al. 2006, likely due to higher metabolic and physiological costs (Meinzer 2003). In general, tall trees may be particularly vulnerable to drought if they exhibit isohydric stomatal regulation . In the case of the isohydric N. dombeyi, the absence of studies monitoring the ecophysiological behavior of adult trees during drought periods hampers our efforts to elucidate the mechanisms underlying plant-level responses. However, studies carried out on seedlings subjected to gradual drought reveal that drought influences the processes associated with carbon gain in N. dombeyi, and the synthesis of reserve compounds (mainly starch) prevails over the consumption of reserves during dry conditions, while photosynthesis is reduced (Piper 2011).
The simulated effect of severe drought events can be noticed even in the case of stands with large trees. When the probability of a 1998-type drought event hitting the N. dombeyi forests is as high as five events per century, stand basal area vs. time tends to become flat (Fig. 7a, b) or even shows hints of a slight drop. Thus, as expected, our simulations revealed the cumulative depletion of both tree density and size after successive drought events. Delayed effects of past episodes of low-growth periods on growth and survival may be complex (Lloret et al. 2011) due to stand modifications imposed by these past mortality events Veblen 1998, Camarero et al. 2011), in addition to genetic, phenotypic, and adaptive characteristics of individual trees (Lloret et al. 2012). Not all this variability could be considered in our study. However, dendrochronological studies on N. dombeyi populations have revealed that previous drought events induced a divergence in growth patterns that eventually caused trees with lower growth rates to become more vulnerable to additional drought events (Suarez et al. 2004, Suarez andKitzberger 2010).
The analysis above showed that the total number of ingrowth individuals was larger under normal conditions than under drought episodes. However, no further conclusions could be drawn from the recruitment datasets alone other than computing the mean size of the ingrowth trees. Generally speaking, the intrinsic spatial and temporal variability of the recruitment stage (e.g., seed germination and emergence, seedling antagonisms) should prevent us from extrapolating those results to study short-and mid-term patterns of recruitment. Nonetheless, our approach enabled us to integrate those recruitment parameters into a model to calculate mid-term projections of size distributions for adult individuals. It is known that gap opening resulting from adult tree deaths is a common mechanism for the recruitment of new cohorts in forests. In our study system, growth release of understory vegetation, including N. dombeyi saplings, has also been documented under dead trees (Suarez and Sasal 2012). Our simulations point to the relevance of large dead trees for the establishment of new cohorts in the mid-and long term, as well as to the vulnerability of new recruits to drought events, particularly if drought recurrence is high (Figs. 7 and 8).

Limitations of the present study
Although our simulations suggest that N. dombeyi forests may decline under recurrent drought episodes, the difficulties we faced when collecting the datasets limited our approach. Consequently, N. dombeyi tree ingrowth, survival, and growth were assumed throughout the study to depend only on dbh. The absence of competition as a limiting factor in the ingrowth, survival, and growth functions did not enable us to draw stronger conclusions about the future dynamics of the simulated tree stands. If competition effects were included, one should likely see enhanced mortality and slower growth rate at higher tree densities, such that curves of basal area per ha (Fig. 7a-c) would gradually reach an asymptotic upper value. Effects on recruitment, on the other hand, are more difficult to identify. Total amount of seeds produced should be approximately proportional to the number of adult individuals available in the stand. However, seedlings would arguably face increasing competition for light and nutrients that might, on the other hand, reduce their germination and/or establishment rates. In addition to competition, other factors (e.g., slope, orientation, altitude, soil characteristics, understory composition and structure, site quality) may likely exert an impact on some or all the processes that determine the evolution of a tree stand, but could not be measured and included in our model. These problems could only be tackled by acquiring better datasets, a task that is left for future studies.
In general, N. dombeyi forests show a rather complex structure with dominant, subdominant, and small trees. In these forests, gap dynamics of those non-dominant individuals may become important when large trees fall and competition for light and nutrients is modified locally. Saplings may also benefit when light conditions improve after an adult tree dies and falls. However, such a complex gap dynamics could not be included in the model. Furthermore, due to data limitations, recruitment has been coarsely determined, which may limit the accuracy of the ingrowth rate. Finally, throughout the study, we assumed that all stands were monospecific. In our experience, this assumption is not unrealistic, although in the field mixed stands of N. dombeyi and A. chilensis often occur. Future work with better datasets may allow us to model those mixed stands with the IPM methodology.

CONCLUSIONS
In the present work, we have integrated basic demographic information on recruitment, growth, and survival into a population dynamics model. In addition to implementing such methodology, our objective has been to study the evolution of N. dombeyi forests in a context of increasing drought episodes. We have applied such approach to dominant N. dombeyi tree stands, found typically in austral temperate ecosystems and for which information about projections of forest dynamics is lacking. Although it would probably be too speculative to suggest a risk of extinction in these populations, the trends uncovered in the simulations point to a situation in which frequent droughts may spell trouble for these populations at stand level. Our models show that, in those cases that concur with future climate scenarios (characterized by an increase in the frequency of severe drought episodes), forests may go through a dramatic decline, from which full recovery may be difficult.