Forest disturbances and climate constrain carbon allocation dynamics in trees.

Forest disturbances such as drought, fire, and logging affect the forest carbon dynamics and the terrestrial carbon sink. Forest mortality after disturbances creates uncertainties that need to be accounted for to understand forest dynamics and their associated C-sink. We combined data from permanent resampling plots and biomass oriented dendroecological plots to estimate time series of annual woody biomass growth (ABI) in several forests. ABI time series were used to benchmark a vegetation model to analyze dynamics in forest productivity and carbon allocation forced by environmental variability. The model implements source and sink limitations explicitly by dynamically constraining carbon allocation of assimilated photosynthates as a function of temperature and moisture. Bias in tree-ring reconstructed ABI increased back in time from data collection and with increasing disturbance intensity. ABI bias ranged from zero, in open stands without recorded mortality, to over 100% in stands with major disturbances such as thinning or snowstorms. Stand leaf area was still lower than in control plots decades after heavy thinning. Disturbances, species life-history strategy and climatic variability affected carbon-partitioning patterns in trees. Resprouting broadleaves reached maximum biomass growth at earlier ages than nonresprouting conifers. Environmental variability and leaf area explained much variability in woody biomass allocation. Effects of stand competition on C-allocation were mediated by changes in stand leaf area except after major disturbances. Divergence between tree-ring estimated and simulated ABI were caused by unaccounted changes in allocation or misrepresentation of some functional process independently of the model calibration approach. Higher disturbance intensity produced greater modifications of the C-allocation pattern, increasing error in reconstructed biomass dynamics. Legacy effects from disturbances decreased model performance and reduce the potential use of ABI as a proxy to net primary productivity. Trait-based dynamics of C-allocation in response to environmental variability need to be refined in vegetation models.

INTRODUCTION

Forests are a major component of the global carbon balance and environmental disturbances are among the most influential factors determining forest carbon dynamics (Anderson‐Teixeira et al., 2021; Friedlingstein et al., 2020; McDowell et al., 2020; Seidl et al., 2017). Currently, forests suffer increased disturbance regimes and can express high vulnerability to drought and hot spells (Brodribb et al., 2020; Choat et al., 2018; McDowell et al., 2018, 2020). The negative effects of global change can be particularly evident in ecosystems exhibiting negative land‐use legacies (Gea‐Izquierdo et al., 2021; Kannenberg et al., 2019; McDowell et al., 2018, 2020). Forest degradation can be caused by recurrent and intense anthropogenic disturbances including logging, fire or overgrazing. The negative effects of anthropogenic disturbances on forests are intimately linked to those of natural disturbances such as drought (Anderegg et al., 2020; Lin et al., 2015; McDowell et al., 2020; Seidl et al., 2017). Changes in disturbance regimes, climate, and global factors such as atmospheric [CO2] affect the carbon balance controlling gross (GPP) and net primary productivity (NPP) (Chen & Luo, 2015; Gower et al., 2001; Vogel et al., 2008; Walker et al., 2019). In the long‐term, net increments in biomass production under enriched [CO2] are generally restricted to young tree stages or constrained by some sink limitation like nutrient availability (De Kauwe et al., 2021; Jiang et al., 2020; Körner, 2015; McDowell et al., 2020; Norby & Zak, 2011). Yet, there are uncertainties at different organizational levels regarding the effect of CO2 particularly under increased water stress (De Kauwe et al., 2021; Friend et al., 2014; Pugh et al., 2016; Sperry et al., 2019; Walker et al., 2021).

Modifications in microenvironmental conditions, plant‐plant competition and stand age after disturbance affect NPP (Anderegg et al., 2020; Besnard et al., 2018; Kannenberg et al., 2019; Zhu, 2020). Disturbances modify stand conditions and the access of surviving trees to light, water and nutrients. Plant species have developed different functional traits to maximize resilience to disturbances, including sprouting behavior (Bond & Midgley, 2001; Pausas et al., 2016; Zeppel et al., 2015). Anthropic disturbances such as thinning aim to enhance individual growth and drought resilience of post‐disturbance remaining trees (Bottero et al., 2021; Fernández‐de‐Uña et al., 2016; Sohn et al., 2016; Tsamir et al., 2019). However, tree responses to stress are scale dependent and individual responses do not necessary apply when upscaled to the stand level (De Kauwe et al., 2021; Ingrisch & Bahn, 2018; Kannenberg et al., 2019). Species resilience depends not only on individual performance but also on demography and regeneration success (Albrich et al., 2020; Gessler et al., 2020; Xu et al., 2017). Consequently, enhanced individual growth resilience to disturbances may not readily translate into enhanced species‐level or ecosystem‐level resilience and enhanced long‐term forest sink capacity (Albrich et al., 2020; Gessler et al., 2020; Ingrisch & Bahn, 2018).

Terrestrial carbon sink dynamics are in part driven by dynamical redistribution of carbon within plants, the carbon use efficiency, the carbon turnover and the carbon residence time (Friend et al., 2014; Malhi et al., 2015; Muller‐Landau et al., 2021; Trugman, Anderegg, Sperry et al., 2019). Different mechanical and hydraulic mechanisms explain how plants modify carbon allocation dynamically in response to environmental forcing (Franklin et al., 2020; Reich et al., 2014; Trugman, Anderegg, Wolfe, et al., 2019; Zuidema et al., 2018). Plants preferentially allocate more biomass to the organ that acquires the most limiting resource (Kannenberg et al., 2019; Reich et al., 2014; Weng et al., 2015). Optimality principles expect plants under stress to allocate more carbon belowground than aboveground (Franklin et al., 2020; Jiang et al., 2020; Schiestl‐Aalto et al., 2015). Therefore, invariable allocation schemes in vegetation models are inadequate because plant allometry depends on resource availability following a functional equilibrium (Lehnebach et al., 2018; Merganičová et al., 2019; Potkay et al., 2021). Implementation of accurate dynamic carbon allocation schemes is key in vegetation models (Franklin et al., 2020; Potkay et al., 2021; Trugman, Anderegg, Wolfe, et al., 2019), particularly because models yield very different output depending on the carbon allocation schemes used (Reich et al., 2014; Sitch et al., 2008; Weng et al., 2015; Zuidema et al., 2018).

There are different allocation principles and theories (Franklin et al., 2020; McMurtrie & Dewar, 2013; Merganičová et al., 2019; Trugman, Anderegg, Sperry et al., 2019). The carbon distributed to growth has been often treated in vegetation models as a constant ratio of photosynthesis (Merganičová et al., 2019). Such simplified allocation rule follows to some extent the pipe‐model and related proportional allocation theories (Chiba, 1998; Le Roux et al., 2001; Mäkelä, 2002; Shinozaki et al., 1964). To address variability in C‐allocation, models add empirical coefficients or mechanistic functions linked to tree age, stand competition or to environmental forcing (Gea‐Izquierdo et al., 2015; Guillemot et al., 2017; Hayat et al., 2017; Peng et al., 2002). It is important to minimize the number of local parameters to increase global applicability while addressing plasticity in functional traits (Anderegg & Venturas, 2020; Sperry et al., 2019). Complex models without single algebraic solutions (i.e., equifinality) need to constrain their parametric spaces within biologically sound ranges (Prentice et al., 2015; Trugman, Anderegg, Sperry et al., 2019; Zuidema et al., 2018). The representation of physiological processes in models has to be continuously refined through model‐data fusion, and model reliability can be improved by fusing with data sets of forest carbon dynamics with wide spatio‐temporal coverage (Babst et al., 2021; Heilman et al., 2022; Peng et al., 2011; Zuidema et al., 2018).

Stand‐level forest carbon estimates can be accomplished using different techniques (Baldocchi, 2014; Walker et al., 2021; Weng et al., 2015). One of them estimates annual series of woody biomass production using allocation equations applied on data from: (1) forest inventories where past mortality is taken into account (permanent sampling plots, PSPs); and, (2) tree‐rings from dendroecological plots (DPs) (Babst et al., 2018; Chen & Luo, 2015; Helcoski et al., 2019). Combining these two data sources can help to overcome specific biases of each sampling type (Bottero et al., 2021; Teets et al., 2018; Walker et al., 2021). DPs follow a C‐oriented sampling scheme to assess stand productivity (Babst et al., 2014, 2018; Brienen et al., 2018; Hember et al., 2019). Woody NPP estimates from DPs can be used to benchmark vegetation models with data in NPP units (g C m−2 year−1) to avoid bias potentially arising when using normalized tree‐ring growth estimates to calibrate NPP model estimates (Gea‐Izquierdo et al., 2017; Muller‐Landau et al., 2021; Zuidema et al., 2018). Additionally, it is necessary to assess bias added by past mortality in the analysis of NPP calculated using tree‐rings (Brienen et al., 2018; Gower et al., 2001; Hember et al., 2019; Walker et al., 2021). Decreasing reliability of tree‐ring data with increasing time before the present (“fading record” problem, Swetnam et al., 1999) can be explained by the impact of past mortality both through missing NPP and modification of carbon distribution in the remaining trees. Yet, tree‐rings are the most effective data‐driven tool to reconstruct forest carbon history particularly when DPs and PSPs are combined (Bowman et al., 2013; Brienen et al., 2018, 2020; Heilman et al., 2022; Hember et al., 2019).

Understanding plant carbon allocation patterns in response to the changing environmental conditions and disturbance regimes is essential to assess the carbon budget under global change (Franklin et al., 2020; Reich et al., 2014). Plants generally follow the Surplus Carbon Hypothesis, namely plant growth is usually constrained by soil nutrients, water or temperature rather than by availability of photosynthates (Prescott et al., 2020). Therefore, C‐sink limitations are the most common limiting plant growth (Körner, 2015). Consequently, models have evolved from the C‐limitation to the C‐surplus hypotheses to better characterize time variability in the C‐sinks (Fatichi et al., 2019; Hayat et al., 2017; Zuidema et al., 2018). To achieve this goal, allocation rules need to be dynamically set as functions of environmental variability as we do in the model‐data fusion approach in this study. Temporal dynamics of aboveground woody biomass increments (ABI) combining DPs and PSPs, which experienced different levels of disturbance were used to benchmark a vegetation model, which in turn was used to analyze variability in patterns of carbon allocation in forests. We analyzed woody biomass trajectories, their response after disturbance and estimated the bias added by past‐mortality. Specifically, we addressed the following: (1) How do climatic constraints and environmental forcing dynamically drive carbon allocation? (2) are carbon allocation rules a function of plant–plant competition and modified by disturbance intensity?

MATERIALS AND METHODS

Study sites: Dendroecological plots within permanent sampling plots

DPs were sampled within PSPs located at different locations in Spain to study dynamics in aboveground woody biomass (Appendix Figure 1). The PSPs were rectangular of different size and originally set to study the effect of thinning (i.e., a stand disturbance) on tree growth and forest productivity (Table 1). We selected PSPs from six different sites corresponding to two different forest types in Mediterranean mountains, namely: (1) coniferous Pinus sylvestris L. stands, a nonresprouting species; (2) hardwood deciduous stands of Quercus pyrenaica Willd. and Quercus faginea Lam, both being resprouting species. All stands were monospecific and coetaneous (Table 1). P. sylvestris is a boreal conifer with a wide Eurasian distribution that finds in the Mediterranean mountains its southern limit. The two Quercus species are ring‐porous sub‐Mediterranean oaks that can form monospecific or mixed stands with other species including Pinus sylvestris, which has lower drought‐tolerance than the oaks (Gea‐Izquierdo et al., 2021; Martín‐Gómez et al., 2017). All plots were located on acidic soils, except for Q. faginea, which were calcareous. Life‐history strategy was different in the oak and pine stands: pines were of seedling origin whereas oaks were established after the last coppicing for firewood in the mid‐1900s (Table 1). All trees were unistemmed except in Barriopedro that were pluristemmed.

TABLE 1

Site characteristics. n = number of trees within each plot. Year0 = year of plot establishment. The subindex in the thinned plots refers to percentage of basal area removed. Mean ABI is for the period 2000–2016. Standard deviations are given between parentheses. Alt, altitude; Diff, LAIthinning‐LAIcontrol; LAI, mean stand tree leaf area index in 2016; Rad, plot radius

Species	Site	Year0	Plot	Alt (m)	n	Rad (m)	dbh (cm)	Height (m)	Age (years)	Plot BA (m2 ha−1)	Density (trees ha−1)	Tree canopy LAI		ABI (g C m−2 year−1)
												(m2 m−2)	Diff (%)
Quercus faginea	BarrioPedro	1980	Thinning50	873	43	14.3	14.2 (3.5)	6.6 (1.5)	64.6 (4.8)	11.3	666.9	0.86	−13.1	29.6 (5.4)
			Control		48	7.3	8.7 (3.1)	4.6 (1.5)	63.3 (7.1)	19.4	2917.4	0.99		44.9 (6.5)
Quercus pyrenaica	Navasfrías	2004	Thinning35	898	35	13.9	21.2 (4.1)	15.2 (2.3)	69.0 (5.1)	21.0	574.0	2.57	−7.6	88.2 (24.2)
			Control		35	10.5	20.0 (4.4)	15.7 (2.2)	67.4 (6.2)	33.1	1005.9	2.78		136.0 (37.9)
	Rascafría	1994	Thinning50	1298	36	14.3	18.8 (3.9)	13.4 (1.9)	60.1 (5.0)	16.1	557.9	2.62	−15.5	79.4 (21.0)
			Control		39	6.9	14.8 (5.2)	12.4 (3.7)	58.6 (6.3)	50.1	2596.9	3.10		172.2 (38.7)
Pinus sylvestris	Navafría	1971	Thinning25	1636	40	13.1	25.6 (3.5)	18.0 (1.5)	68.2 (5.0)	38.9	742.4	2.84	−11.8	175.2 (17.2)
			Control		48	9.3	21.9 (6.0)	19.1 (2.3)	70.7 (8.6)	71.9	1779.8	3.22		246.8 (29.5)
	Covaleda	1968	Thinning25	1530	44	14.0	28.0 (4.2)	17.1 (0.9)	92.9 (4.7)	44.5	706.9	2.29	−0.4	117.7 (16.1)
			Control		40	9.8	21.6 (5.6)	16.0 (1.8)	86.5 (7.2)	51.3	1320.5	2.30		121.6 (13.1)
	Duruelo	1968	Thinning50	1206	41	14.1	27.6 (4.7)	16.8 (1.2)	79.5 (3.3)	37.1	631.2	2.24	−20.6	117.4 (17.5)
			Control		40	9.0	20.9 (5.1)	16.4 (1.9)	77.0 (7.4)	58.4	1612.1	2.82		153.7 (23.7)

We sampled 12 plots from the six sites: six control plots with natural stand dynamics, and six plots within the maximum thinning intensity applied at each site (Table 1). With these data we could assess the effect of two major disturbances affecting stand dynamics: (1) thinning at the six sites; (2) a snowstorm in winter in 1996 in one site (Navafría). Thinning was applied at the moment of plot establishment (between 1968 and 2004. Table 1; Figure 1) and woody stems over 7.5 cm in diameter at breast height (DBH) were identified and tagged. Tree diameter and height were repeatedly measured at varying intervals every 5–10 years. PSP data ranged between 3 and 9 censuses, with a total monitoring interval between 11 and 45 years until the last census. Data ranged between 14 and 49 years from the PSPs establishment to the DPs sampling date in 2016. We fitted height‐diameter non‐linear models with DBH and height data from PSPs to be used in combination with species‐specific allometric equations to calculate biomass in PSPs and DPs. Using linear interpolations between periodic measurements within PSPs, we estimated annual time series of biomass loss accounted by past mortality and estimates of woody biomass increments. These estimates of biomass dynamics from PSPs included biomass gain by ingrowth of surviving trees and new recruits, and biomass loss due to tree mortality. Stand‐level biomass was calculated by summing the biomass of all trees within each sample plot for each census (Chen & Luo, 2015).

FIGURE 1

Aboveground woody (i.e., stem and branches up to 2 mm diameter) biomass increment (ABI, g C m−2 year−1) chronologies of the six control (blue) and thinning (red) plots. Gray polygons show the period after thinning. ABI calculated using just standing trees sampled in 2016 are depicted with thick lines, whereas estimates combining ABI and permanent plot data (ABIcorr) to take into account the influence of the “fading record” are depicted with thin lines. Blue (for control) and orange (for thinning) background polygons represent the estimated bias included when neglecting the influence of “ghost‐trees.” We indicate for Navafría the 1996 snowstorm (disturbance) producing a mortality event

Additionally, in 2016 we collected data from 12 ad‐hoc DPs located at the center of each PSP (Davis et al., 2009; Graumlich et al., 1989). DPs were circular and of variable radius where every tree with DBH greater than 7.5 cm was mapped, cored twice at 1.3 m and measured for DBH and height (Table 1). DPs included between 35 and 48 trees, which is considered representative for biomass dynamics in the studied forests (Babst et al., 2014, 2018; Graumlich et al., 1989). Finally, we systematically estimated leaf area index (LAI) at five equidistant points along one diagonal of each PSP using a plant canopy analyzer LICOR 2200C (i.e., we report plant area index; Table 1).

Assessment of unbiased woody biomass dynamics combining PSPs and DPs

Cores from DPs were processed to enhance ring‐width (RW) identification, visually and statistically cross‐dated and annual radial growth measured (Fritts, 1976). Tree individual radial growth data (Appendix Figures 2 and 3) was used to calculate annual estimates of C biomass increments to calculate ABI (see below), and C biomass of woody roots components up to 2 cm in diameter used together with ABI to calculate cumulative woody biomass. We reconstructed tree diameters from ring‐widths, proportionally to measured DBHs and bark ratios (Babst et al., 2014; Foster et al., 2014; Graumlich et al., 1989). Then, we calculated individual woody biomass annual growth considering a 50% of carbon content, using calculated diameters and estimated heights as inputs in allometric equations of DBH or DBH and height (Ruiz‐Peinado et al., 2011, 2012). Tree annual time series of aboveground biomass increments were integrated at the plot level and are reported along the manuscript as ABI (g C m−2 year−1) (Graumlich et al., 1989). ABIs are estimates of aboveground carbon accumulated in woody biomass, that is, proxies of net aboveground woody productivity (NPPAW, Malhi et al., 2015; Knapp et al., 2017). To characterize trends in woody biomass allocation we estimated peak biomass increment and stand age at peak biomass using 20‐year splines on the ABI chronologies (Foster et al., 2014).

Biomass and stand estimates from unique sampling schemes include growing bias (underestimation) with increasing time into the past resulting from unknown past mortality (“ghost trees” or “fading record”; Foster et al., 2014; Swetnam et al., 1999). Observations from PSPs were used to estimate and correct this bias (Chen & Luo, 2015; Helcoski et al., 2019; Teets et al., 2018). Based on PSP data we built time series of estimated woody biomass loss back in time by estimating the percentage of biomass lost every year from both natural mortality and thinning. Then annual woody biomass and ABI estimated from DPs were corrected (ABIcorr). Finally, combining ABI and ABIcorr (Figure 1), we assessed bias in ABI calculated from DPs and estimated the C biomass lost back in time to discuss the influence of climate and disturbance on C‐allocation variability (Davis et al., 2009; Foster et al., 2014; Helcoski et al., 2019).

Exploring sensitivity to climate of carbon allocation to woody biomass in trees

For the different analyses we combined daily climatic data from http://www.meteo.unican.es/datasets/spain02 (Herrera et al., 2016) for 1950–2007 and ERA‐interim reanalysis data (http://www.ecmwf.int/en/research/climate‐reanalysis/era‐interim) for 2008–2016. Additionally, we calculated the standardized precipitation‐evapotranspiration index (SPEI) (Beguería et al., 2014) with a 6‐ and 12‐month lag. All analyses apply to the shared 50‐year period 1967–2016. The six sites ranged in mean annual precipitation between 1131 mm and 532 mm, and in mean annual temperature between 8.3°C and 14.1°C (Appendix Table 1).

We standardized ABI and ABIcorr time series to be used in the climate response analysis and model calibration (Teets et al., 2018). To remove the long‐term trend related to stand age and time we used flexible 20‐year splines and ratios centered to mean ABI of the last 25 years (Hember et al., 2019. Appendix Figure 4). The 12 ABI chronologies did not include in the last 25 years any trend potentially related to juvenile growth and ontogeny or the maximum early peak in broadleaves (Foster et al., 2014. Appendix Figure 4). ABI and ABIcorr standardized estimates are both in g C m−2 year−1, therefore in NPPAW units. We explored the sensitivity to climate of carbon allocation to woody biomass in control plots using bootstrap Pearson correlations.

Dynamic carbon allocation as a function of climatic variability and disturbance

In addition to exploring the empirical linear relationship between ABI and climate, we assessed the role of climatic variability and disturbance on C‐allocation using a vegetation species‐specific process‐based model where the water and carbon balances are simulated at the stand level (MAIDEN; Gea‐Izquierdo et al., 2015; Misson, 2004). Following the functional balance principle, C‐allocation rules as functions of environmental variability in the model express sink‐limitations (Körner, 2015) directly related to temperature and moisture constraints, and indirectly reflecting internal functions such as cell turgor (Fatichi et al., 2019; Franklin et al., 2020; Merganičová et al., 2019). Transpiration, GPP, and NPP are calculated based on a classic coupled photosynthesis‐stomatal conductance model (Farquhar et al., 1980). Then NPP is allocated to four different plant compartments (canopy, stem, roots, and non‐structural carbohydrates, NSCs) using different dynamic non‐linear functions of temperature and moisture (soil and atmosphere) along five different phenological stages (see Gea‐Izquierdo et al., 2015 and Appendix Table 2 for more details). This is an extension of regularly used empirical estimates of secondary growth as a function of climate but instead of statically and empirically, C‐allocation to the stem is directly derived dynamically with biologically sound non‐linear functions of environmental variability that can therefore fit both C‐source limitations (C supply from photosynthesis) and C‐sink limitations (Fatichi et al., 2019; Körner, 2015). The model lacks specific allocation to root exudates or secondary C metabolites (Jiang et al., 2020), which are allocated within the bulk NSC pool.

We used the model to assess whether carbon allocation is a dynamic function of climatic variability and disturbance. In our approach, by fusioning the same process and data (i.e., modeled and estimated NPPAW, ABI and ABIcorr, in g C m−2 year−1), we improve model performance and avoid bias associated to using unitless normalized growth indices to calibrate NPP output from vegetation models (Babst et al., 2018; Gea‐Izquierdo et al., 2017; Prentice et al., 2015). Along the study, poorer model performance is interpreted as a departure from a direct relationship between carbon allocation and climatic variability. To assess the functional relationship between climatic variability and C‐allocation, and discuss data and model bias, we compared maximum likelihood (ML) and Bayesian methods (Gaucherel et al., 2008; Hartig et al., 2011) to calibrate nine allocation parameters in MAIDEN accounting for different phenological stages (Appendix Table 2). First, we calibrated the nine local parameters constrained within ecologically realistic ranges using a global optimization algorithm and ML principles as in Gea‐Izquierdo et al. (2017) both to ABI and ABIcorr from control plots. Then, we explored a posteriori distribution of the same nine parameters using a Bayesian inference and a Markov chain Monte Carlo (MCMC) DREAMsz algorithm, which has proofed to be best for complex multimodal models compared with other MCMC samplers (Hartig et al., 2019; Lu et al., 2017; Rezsöhazy et al., 2020). We ran a minimum of three MCMC chains of 100,000 iterations each with uninformative uniform priors within a sound range for each parameter (Gennaretti et al., 2017). We report the maximum a posteriori (MAP) estimates of the MCMC and check convergence using the univariate and multivariate Gelman–Rubin statistics () considering a threshold of 1.2 for convergence (Hartig et al., 2019; Lu et al., 2017). Then, using the same MAP and ML parameters calibrated for the control plots, we run the model in the thinning plots just modifying LAI to assess model performance and how C‐allocation varies after disturbance. The goodness‐of‐fit of models was assessed through the coefficient of determination (R 2), linear correlation, and bias (Appendix Table 2).

RESULTS

Trends in woody biomass depend on tree life‐history and disturbance

Woody biomass increment of the studied broadleaves exhibited different trends than those of the conifers (Figure 1). Age at peak biomass was more than twice older in pines than in oaks (Appendix Table 3), meaning that oaks expressed the maximum in ABI earlier (at mean age 21.3 years) than in pines (at mean age 54.8 years) (Figure 1; Appendix Table 3). After reaching maximum ABI in oaks, mean ABI declined for about 20 years before it stabilized around a sitewise mean value. In pines, after the peak in biomass increment, mean ABI stabilized or slightly decreased. This suggests the existence of allocation trends in ABI related to life‐history strategy or regeneration mode (resprouters vs. nonresprouters) that we had to take into account to analyze ABI time series. Yet, the observed ABI profiles did not suggest a single general allocation rule. For instance, in BarrioPedro there was a hump after thinning resulting on greater ABI during circa 10 years than in stands with higher basal area and leaf area (i.e., thinned vs. control plot). In contrast, a hump in ABI after thinning was not observed in the other sites, which followed a positive relationship between ABI and plot basal area (Figure 3a). Increases in ABI after reductions in basal area expressed that C‐allocation among plant compartments were not always proportional. For this reason, we standardized ABI and ABIcorr for the sensitivity to climate and modeling analyses.

FIGURE 3

Plot Basal Area and: (a) aboveground woody (stem + branches < 2 mm diameter) biomass increment (ABI); (b) modeled NPP; (c) leaf area index (LAI). Trends are depicted for all observations in black, in red for oaks and blue for pines. Coefficient of determination (R 2) for the different fits are shown. In (a) and (b) annual data are shown for the whole period since each permanent plot was established (see Table 1). In (c) observations correspond to the last year (i.e., LAI in 2016, BA from the last inventory available)

Woody biomass steadily increased along stand age but did not reach a plateau in any of the forests studied (Figure 2). Woody biomass and hence modeled NPP, steadily increased with increasing stand basal area (Figure 3). LAI increased non‐linearly with stand basal area, but this relationship did not reach a steady state within the observed data range (Figure 3c). LAI did not recover decades after thinning under the heaviest thinning treatments (Table 1). Lower LAI after thinning consequently resulted on lower simulated forest GPP and NPP in thinned than in control stands (Table 2). The carbon sink capacity of tree stands was negatively affected by disturbance intensity.

FIGURE 2

Woody biomass (WB, stem, branches up to 2 mm diameter and coarse roots) estimates for the six sites for control and thinning plots. Gray polygons show the period after thinning. Lower (blue or red thick lines) polygon borders depict WB calculated using trees present at the sampling date. Upper polygon borders (blue or red thin lines) show estimates including past mortality as recorded in the permanent plots to take into account the influence of the “fading record” in biomass estimates (i.e., WBcorr). Note that in all cases we assume nil mortality before thinning because we lacked those data. The (blue and red) polygons show increasing bias with time before sampling date if using just standing trees at the sampling date to estimate woody biomass time series. We indicate for Navafría the 1996 snowstorm (disturbance) producing a mortality event

TABLE 2

Simulated carbon balances in control and thinned plots since year 2004 (i.e., the most recent thinning, in Navasfrías; see Table 1). Standard deviations are shown between parentheses. Carbon use efficiency, CUE = NPP/GPP. ANPP = aboveground forest NPP (leafs and wood ingrowth). NPP and GPP refer only to the tree layer, without accounting for understory NPP

Site	Plot	GPP (g C m−2 year−1)	NPP (g C m−2 year−1)	ANPP (g C m−2 year−1)	CUE
BarrioPedro	Control	823.3 (71.3)	437.8 (37.9)	114.9 (7.4)	0.532
Navasfrías		1622.8 (169.3)	862.4 (89.9)	358.8 (23.2)	0.531
Rascafría		1603.7 (132.2)	856.4 (70.8)	425.8 (16.4)	0.534
Navafría		1778.9 (112.5)	966.0 (59.3)	369.0 (23.4)	0.543
Covaleda		1556.6 (132.8)	839.4 (70.1)	251.6 (14.6)	0.539
Duruelo		1883.6 (177.4)	987.8 (93.8)	306.2 (17.9)	0.539
BarrioPedro	Thinning50	770.4 (65.5)	409.7 (34.8)	103.9 (7.6)	0.532
Navasfrías	Thinning35	1392.1 (135.9)	739.7 (72.1)	220.2 (12.7)	0.531
Rascafría	Thinning50	962.6 (54.8)	514.0 (28.7)	157.9 (9.1)	0.534
Navafría	Thinning25	1245.9 (75.7)	674.3 (39.7)	238.0 (6.0)	0.541
Covaleda	Thinning25	1445.8 (122.0)	779.3 (64.5)	220.4 (7.7)	0.540
Duruelo	Thinning50	1358.1 (95.4)	730.1 (50.3)	189.5 (5.2)	0.537
Functional type or treatment
Broadleaves		1195.8 (390.2)	636.7 (207.8)	230.3 (133.7)	0.532 (0.001)
Conifers		1544.8 (246.5)	829.5 (126.7)	262.5 (64.7)	0.540 (0.002)
Control		1544.8 (374.1)	825.0 (199.4)	304.4 (110.1)	0.536 (0.005)
Thinning		1195.8 (270.3)	641.2 (146.7)	188.3 (50.1)	0.536 (0.004)
Control − Thinning	(trait units)	349.0 (248.5)	183.8 (130.8)	116.1 (91.7)	0.001 (0.001)
	(%)	+22.6	+22.3	+38.1	+0.1

The disturbance intensity was expressed by the estimates of neglected past mortality (“fading record”), which were higher after heavy thinning and after the 1996 snowstorm with more than 30% of forest biomass loss. Bias in past biomass estimates increased with increasing disturbance intensity. In average, neglected mortality was below 10% up to 18 years before sampling date, and increased to 20% in woody biomass reconstructions 25 years before sampling date (Figure 4). Figures 1 and 2 show bias in ABI and woody biomass chronologies. In the absence of a major disturbance, the past natural mortality observed in PSPs was reduced or almost nil for the studied period and forests, and focused on suppressed individuals with very slow growth and little influence on stand NPPAW. Past‐mortality was not a function of extant basal area (Figure 4).

FIGURE 4

Estimate of the influence of the “fading record” (neglected mortality) to depict the increasing bias with time included in dendrochronological stand biomass chronologies. (a) Estimated bias (in % of aboveground woody biomass, AWB) as a function of distance (time) before sampling date (year 2016, t = 0) at the six sites and 12 plots studied calculated using repeated inventories since plot set up. Same as in (a) but as a function of plot basal area for the 12 sites for the common period with permanent plot data (i.e., since 2004) (b); and only for the six coniferous plots (i.e., since 1971) (c). Vertical dashed lines in (a) depict years (ys) when plot chronology mean bias is below 1%, 5%, 10%, 12%, and 20%

Carbon allocation dynamics in response to climatic variability and disturbance

Climatic variability partly controlled carbon allocation dynamics to woody biomass, expressing both the limiting effect of winter low temperatures and water stress in the active spring‐summer season (Figure 5). Woody biomass increments were positively influenced by warmer winters and negatively by winter precipitation in the evergreen pines, whereas winter climate did not affect ABI in the deciduous oaks (Figure 5). Humid springs and humid March–August (as expressed by SPEI6Aug) increased ABI in four out of six of the studied forests. The relationships between ABI and drought at different time lags and seasons suggested phenological differences among sites. The most sensitive to drought was BarrioPedro, which is the only calcareous site, with the driest climate at the lowest elevation of all studied forests (Appendix Table 1). Additionally, this was the site expressing a differential woody allocation (ABI) hump after thinning (Figure 1). ABI was more reduced under drier conditions, at longer (12 months) or earlier time periods (previous year fall and winter) in oaks than in pines. ABI in pines from Navafria was particularly reduced by SPEI6April (Figure 5), which suggests a negative effect of cold and humid winters (i.e., heavier snowstorms) at higher elevations (Table 1). In this sense, ABI from the coniferous site at the lowest elevation (Duruelo) was the only expressing a constraining effect of spring‐summer heat on biomass growth and exhibited the highest sensitivity to drought among conifers (Figure 5).

FIGURE 5

Bootstrap correlations between ABI chronologies for control plots detrended to remove the juvenile trend and selected climatic covariates for periods of maximum growth response for the studied species (shared period 1966–2016). JFM, January–February–March; AMJ, April–May–June; MJJ, May–June–July; Jul, July; Nov, November; Apr, April; Aug, August. Significant thresholds are represented outside the background gray color. T max, maximum temperature; T min, minimum temperature; P, precipitation; SPEI, standardized precipitation evapotranspiration index for either a lag of 6 or 12 months

The sensitivity to climate expressed in ABI was concomitantly reflected by model performance (Table 3, Figure 6). Multiparametric convergence ( < 1.2; Table 3) was achieved in all models except for ABIcorr in Covaleda, although univariate diagnostics were all below 1.2. Just for this site and the parameters that did not converge otherwise we report results for three 100,000‐iteration MCMC runs using informative normal priors with means coinciding with the estimated ML parameters. The poorest model performance in control plots was observed in the only forest suffering a major natural disturbance (Navafría, a snowstorm). The best model performance was expressed at the most water stressed forest (BarrioPedro), where in addition the influence of “ghost‐trees” was the most reduced of all sites (Table 3; Appendix Table 4; Figures 1, 2, and 6). The goodness of fit between modelled and ABI data was similar using either ML or MCMC (Figure 6). Parameters fitted using ML and MCMC, and also MAP distributions estimated using MCMC for ABI and ABIcorr were often similar (Figure 7). The model was calibrated to data from control plots, meaning that model performance was better in control than in thinning plots where the model was applied with the same parameterization from control plots (Tables 1 and 3). Yet, ABI models performed well also for thinning plots except under high intensity disturbance. In all cases, model performance was enhanced when the influence of “ghost‐trees” (i.e., models for ABIcorr) was taken into account (Figure 6; Table 3).

TABLE 3

Goodness of fit statistics of fitted MCMC models to ABI and ABIcorr data from control plots (period 1967–2016) and models with same parameters as in control plots but applied in thinning plots just changing LAI accordingly for each PSP. ABI, aboveground woody biomass; ABIcorr, ABI corrected taking into account past mortality (“fading record”); n, number of forest sites; R 2, coefficient of determination; r, coefficient of correlation. See Section 2 for details on model fit and smoothing for coppice stands (oaks) and high forest conifers (pines). Significant r‐values at α = 0.05 are represented with an asterisk*.  = multivariate Gelman–Rubin statistic

MCMC			Control plots					Thinning plots
Model	Site	Species	R 2	r	Bias		R^	R 2	r
					(g C m−2 year−1)	(%)
ABI	BarrioPedro	QUFG	0.528	0.739*	0.981	2.375	1.05	0.256	0.526*
	Navasfrías	QUPY	0.155	0.396*	1.263	0.994	1.05	0.081	0.307*
	Rascafría	QUPY	0.280	0.530*	−0.661	−0.373	1.03	0.099	0.325*
	Navafría	PISY	0.069	0.451*	0.648	0.272	1.19	−0.208	0.104
	Covaleda	PISY	0.313	0.563*	−0.475	−0.349	1.19	0.239	0.508*
	Duruelo	PISY	0.359	0.601*	−0.707	−0.429	1.09	−0.023	0.256
ABIcorr	BarrioPedro	QUFG	0.525	0.725*	0.017	0.039	1.08	0.445	0.691*
	Navasfrías	QUPY	0.145	0.389*	2.088	1.624	1.03	0.175	0.420*
	Rascafría	QUPY	0.254	0.505*	−0.359	−0.187	1.03	0.157	0.396*
	Navafría	PISY	0.196	0.490*	1.707	0.578	1.02	−0.120	0.143
	Covaleda	PISY	0.359	0.605*	−0.449	−0.299	1.31	0.237	0.531*
	Duruelo	PISY	0.356	0.600*	1.158	0.692	1.10	0.204	0.453*

FIGURE 6

Boxplot with mean differences in correlations and coefficients of determination (R 2 as from Table 3) of models fitted to ABI and ABIcorr in control and thinning plots. Different letters express significant differences for the means

FIGURE 7

Estimated marginal probability density functions (PPDFs) of parameters (see Appendix Table 2) from MCMC with best values both from the maximum a posteriori likelihood (MAP) estimates and the estimated maximum likelihood parameters for ABI and ABIcorr

DISCUSSION

Forest disturbance intensity constrains the forest carbon dynamics

Plant growth is generally constrained by sink rather than by source limitations, meaning that excluding shaded conditions, forest growth is curtailed by water, nutrients or temperature, and not by photosynthesis (Fatichi et al., 2019; Körner, 2015; Prescott et al., 2020). Natural and anthropic disturbances reduce forest leaf area and modify micrometeorological conditions within stands, reducing the short‐term forest carbon sequestration capacity (Curtis & Gough, 2018; Lin et al., 2015; Riutta et al., 2018; Seidl et al., 2017). Recovery of predisturbance biomass levels is inversely related to disturbance intensity (Del Campo et al., 2019; Misson et al., 2005; Vesala et al., 2005; Volkova et al., 2018). After low‐intensity disturbances, forests may be able to quickly restore predisturbance NPP and recover their carbon sequestration capacity. Conversely, forest biomass may not recover after disturbances of high intensity (Amiro et al., 2010; Davis et al., 2009; Lin et al., 2015; Riutta et al., 2018). We observed little tree regeneration and both standing tree biomass and leaf area were still lower decades after heavy thinning. In the absence of disturbances, LAI can be considered stable or to experience small variations despite background mortality (Pappas et al., 2020). There was an overall positive relationship between NPP, stand basal area, and stand leaf area (Sperry et al., 2019; Zweifel et al., 2020). As a result, forest NPP and GPP were more than 20% lower in thinned stands where predisturbance NPPAW did not recover.

The observed lower tree biomass years after disturbance may not necessarily imply a reduced C‐sink capacity at the ecosystem scale (Misson et al., 2007; Pereira et al., 2007; Volkova et al., 2018). Ecosystem leaf area and net ecosystem productivity (NEP) can be enhanced by [CO2] during early stages in stand development or until there is limitation on complementary resources such as nutrients (De Kauwe et al., 2021; Knauer et al., 2017; McDowell et al., 2020; Norby & Zak, 2011). Additionally, lower forest canopy closure leads to an increase in understory cover, adding to aboveground and belowground carbon pools (Misson et al., 2007; Pereira et al., 2007; Terrer et al., 2021). Yet, a reduction in the woody to green biomass ratio concomitantly reduces C retention time and C residence times, while increases C turnover rate (Friend et al., 2014; Muller‐Landau et al., 2021; Walker et al., 2021; Weng et al., 2015). Therefore, regardless of whether NEP was modified or not after disturbance, the forest C‐cycle was modified and expressed a long‐lasting land‐use legacy in the study forests (Curtis & Gough, 2018; McDowell et al., 2020; Seidl et al., 2017).

Trees modify carbon allocation in response to forest disturbances and climate

Stand woody biomass stocks did not reach a steady state in any of our study forests, and self‐thinning did not preclude a steady increase in standing tree biomass for decades even in the oldest stand (Anderson‐Teixeira et al., 2021; Chen & Luo, 2015; Xu et al., 2012). Variability in woody biomass was partly explained by changes in C‐allocation related to differences in species life‐history, climate and disturbance intensity (Graumlich et al., 1989; Trugman et al., 2018; Trugman, Anderegg, Sperry et al., 2019; Xu et al., 2012); similar to the mixed non‐linear signal expressed within dendroecological growth data (Alexander et al., 2018; Foster et al., 2016; Fritts, 1976; Gea‐Izquierdo et al., 2015). Variability in ABI was strongly related to drought particularly in the driest site (Knapp et al., 2017). Stand woody biomass recovery after disturbance depends on ecosystem conditions, with more water‐limited ecosystems showing a slower recovery or even no recovery after disturbance (Anderegg et al., 2020; McDowell et al., 2020; Seidl et al., 2017; Zhu, 2020).

ABI increased for several decades until it peaked and then decreased. These dynamics are expected because of the age‐related trends in NEP and NPP (Anderson‐Teixeira et al., 2021; Besnard et al., 2018; Magnani et al., 2000; Stephenson et al., 2014). The observed biomass growth trajectories in pines were similar to those published in the literature for different species (Babst et al., 2014; Foster et al., 2014; Klesse et al., 2016) but differed from those of oaks. Biomass growth in resprouting oaks peaked on average over 30 years earlier than in pines. Resprouting is a key life‐history functional trait for many species in disturbance prone environments (Bond & Midgley, 2001; Pausas et al., 2016; Zeppel et al., 2015). Difference in ABI trajectories between oaks and pines could be explained by allocation strategy between resprouts and seedlings. The root system already existent in resprouts allows individuals to invest more in aerial tissues within the first years than in seedlings (Bond & Midgley, 2001; Pausas et al., 2016; Zeppel et al., 2015). Nevertheless, faster early growth could diminish the lifespan of trees if there is a trade‐off between growth rates and longevity, and different growing rates in seedlings and sprouts could influence the posterior life expectancy of stems (Brienen et al., 2020; Bugmann & Bigler, 2011; Munné‐Bosch, 2018; Piovesan & Biondi, 2021; Searle & Chen, 2018). Trends in C‐allocation related to ontogeny and age have implications for forest dynamics and the terrestrial carbon cycle (Foster et al., 2014, 2016; Franklin et al., 2020; Zhu, 2020).

In the absence of major disturbances tree biomass and leaf area adjust to the maximum potential of available resources within sites (Breda et al., 2006; Magnani et al., 2000; Zweifel et al., 2020). Overall, good performance of models calibrated for control plots and then applied in thinning plots shows that moderate thinning induced little changes to C‐allocation rules, which were proportional to the leaf area and thus greatly in accordance with the pipe‐model theory (Chiba, 1998; Lehnebach et al., 2018; Mäkelä, 2002; Shinozaki et al., 1964). Yet, with high disturbance intensity or under high stress, a direct proportionality between leaf area and C‐allocation can be lost (Lehnebach et al., 2018). Consequently, trees modified carbon allocation rules under the most intense disturbances, therefore with major changes in stand competition affecting radial growth rates, aboveground:belowground and leaf:wood ratios (Breda et al., 2006; Sohn et al., 2016; Trugman, Anderegg, Sperry et al., 2019; Zweifel et al., 2020). This was shown by the hump after thinning in BarrioPedro and by tree C pools not recovering years after a snowstorm in Navafria. At that site, both model performance and the correlation between ABI and climate were low. This suggests the existence of changes in allocation after disturbance that we did not successfully address with the dynamic allocation rules implemented. Additionally, negative effects of low temperature in winter on ABI in this forest could indirectly express slower acclimation and lower photosynthesintatic rates in winter, and the direct reduction in carbon pools produced by mortality after the snowstorm (Gea‐Izquierdo et al., 2010; Gessler et al., 2020). A similar response to winter climate in evergreen conifers seems to be the rule in Mediterranean mountains (Gea‐Izquierdo et al., 2021; Martin‐Benito et al., 2010; Martín‐Gómez et al., 2017).

Can time series of woody biomass estimates be used to model forest C‐dynamics?

Age and demographic biases embedded by ghost‐trees in NPPAW estimated using tree‐ring data must be carefully taken into account (Alexander et al., 2018; Bowman et al., 2013; Brienen et al., 2018). DP data were used to avoid bias in individual tree sampling and discuss woody C dynamics at the stand level (Babst et al., 2014, 2018). Some studies combining PSPs and DPs suggested that ABI is accurately represented by tree rings for periods of 40 years or longer (Dye et al., 2016; Klesse et al., 2016). In our case, mean bias exceeded 40% at 40 years before sampling date. Accurate representation of biomass growth back in time using data from DP very much depends on the existence of disturbances. Average bias in ABI was below 12% up to 20 years before sampling, a period that seems reasonable to assume as a time window where DP data can be used in the studied forests.

BAI and ABI share almost identical high‐frequency information, thus they can be both used as proxies to study interannual variability in woody productivity (Alexander et al., 2018; Babst et al., 2014, 2018; Bouriaud et al., 2005). Yet, ABI and BAI express different functional scales. This means that, unlike estimates of NPPAW using radial growth data (i.e., BAI or RW), estimates of NPPAW and C‐allocation from vegetation models calibrated using directly ABI are unbiased (Gea‐Izquierdo et al., 2017; Prentice et al., 2015). Most marginal distributions of parameters were leptokurtic and unimodal; therefore, model parameter estimations and consequently model simulations were robust. Despite some differences in parameters estimated with ML and MCMC, model performances were very similar. This shows that, in order to maximize model performance, both model structure and data used to benchmark models are indeed more important than calibration methods (Prentice et al., 2015; Zuidema et al., 2018). Vegetation models need improvement at different levels. For instance, we did not address sink limitations related to nutrients; hence, the model used may not properly address downregulation of different traits with N‐P limitation (Pappas et al., 2020). Additionally, a single global rule in our vegetation model was not enough to model the effect of life‐history traits in C‐allocation (Pausas et al., 2016; Zeppel et al., 2015). There are uncertainties in our analyses related to the model but also to the data, for example, the relationship between aboveground NPP and total NPP (Gower et al., 2001; Pappas et al., 2020). The model may not totally address the surplus carbon hypothesis since the C balance depends on the carbon use efficiency, which was modeled only with a temperature dependent relationship independently of other environmental factors or functional cues (DeLucia et al., 2007; Gea‐Izquierdo et al., 2017). Considering all these uncertainties, ABI was successfully used to improve model performance and representativity (Babst et al., 2018, 2021; Prentice et al., 2015; Zuidema et al., 2018). Further improvements are necessary to address the effect of the different factors modifying carbon allocation, including the effect of disturbances, life‐history traits and climate (McDowell et al., 2020). This should be done further addressing models to the ecohydrologic equilibrium, linking hydraulic mechanisms to C‐allocation rules (Anderegg & Venturas, 2020; Potkay et al., 2021; Sperry et al., 2019).

CONCLUSIONS

Disturbance intensity, life‐history traits, and climate modified the C‐allocation pattern of trees; hence, their effect needs to be implemented in C‐allocation rules within forest vegetation models. Much variability in C‐allocation among stands was explained by differences in leaf area. We successfully modeled woody biomass dynamics using a vegetation model with dynamic carbon allocation as a function of environmental variability particularly in forests expressing a high sensitivity to drought under light disturbance regimes. The effect of disturbance on biomass dynamics and C‐allocation increased with disturbance intensity, thus modeling success decreased in the most highly disturbed sites. The calibration technique did not affect the ecological interpretation of the observed changes in allocation. Yet, uncertainty in parameter estimation related to complex multiparametric spaces in vegetation models needs to be carefully addressed. Stand‐level and individual‐level responses to environmental variability do not necessarily match. Using ABI directly to benchmark vegetation models improved model performance and assessment of forest carbon dynamics.

CONFLICT OF INTEREST

The authors declare no conflict of interest.

Supplementary Material

Click here for additional data file.