Assessing the response of forest productivity to climate extremes in Switzerland using model–data fusion

Abstract The response of forest productivity to climate extremes strongly depends on ambient environmental and site conditions. To better understand these relationships at a regional scale, we used nearly 800 observation years from 271 permanent long‐term forest monitoring plots across Switzerland, obtained between 1980 and 2017. We assimilated these data into the 3‐PG forest ecosystem model using Bayesian inference, reducing the bias of model predictions from 14% to 5% for forest stem carbon stocks and from 45% to 9% for stem carbon stock changes. We then estimated the productivity of forests dominated by Picea abies and Fagus sylvatica for the period of 1960–2018, and tested for productivity shifts in response to climate along elevational gradient and in extreme years. Simulated net primary productivity (NPP) decreased with elevation (2.86 ± 0.006 Mg C ha−1 year−1 km−1 for P. abies and 0.93 ± 0.010 Mg C ha−1 year−1 km−1 for F. sylvatica). During warm–dry extremes, simulated NPP for both species increased at higher and decreased at lower elevations, with reductions in NPP of more than 25% for up to 21% of the potential species distribution range in Switzerland. Reduced plant water availability had a stronger effect on NPP than temperature during warm‐dry extremes. Importantly, cold–dry extremes had negative impacts on regional forest NPP comparable to warm–dry extremes. Overall, our calibrated model suggests that the response of forest productivity to climate extremes is more complex than simple shift toward higher elevation. Such robust estimates of NPP are key for increasing our understanding of forests ecosystems carbon dynamics under climate extremes.


| INTRODUC TI ON
Forests provide a wide range of ecosystem functions and services from global to local scale (Brockerhoff et al., 2017). It is therefore essential to understand how forest ecosystem productivity responds to climate extremes across environmental gradients Cramer et al., 2001;Reichstein et al., 2013) and how those responses feed back to the climate system (Humphrey et al., 2018). Climate change can affect forests on various levels, for example, by modifying the balance and interactions between direct abiotic constraints on tree growth (Cuny et al., 2019), shifting the timing of the growing season (Bigler & Bugmann, 2018), or altering disturbance regimes Sommerfeld et al., 2018). Large-scale variations in forest ecosystem productivity have been primarily attributed to interactions between environmental constraints, namely temperature, water availability and demand and radiation (Beer et al., 2010;Jung et al., 2017;Seddon, Macias-Fauria, Long, Benz, & Willis, 2016), rather than to a single one. In particular, global warming amplifies water limitation as a key constraint for global forest ecosystem productivity, and the spatial extent of drought-limited areas is increasing (Allen et al., 2010;Babst et al., 2019;D'Orangeville et al., 2018;Nemani et al., 2003).
A diverse set of methods is currently used to quantify and project the impact of changing environmental constraints on forest ecosystem productivity, including extensive collections of in situ observations Charney et al., 2016;Clark et al., 2001;Klesse et al., 2018;Shestakova et al., 2019), remote sensing data (Beer et al., 2010;Jolly, Dobbertin, Zimmermann, & Reichstein, 2005;Nemani et al., 2003;Piao et al., 2014), or dynamic vegetation models (DVMs, e.g., Huang, Gerber, Huang, & Lichstein, 2016;Rollinson et al., 2017;Zhang et al., 2018). These and other studies identified important differences in the response of forests to environmental constraints, depending on ambient climate conditions. Forests growing in cold environments at high elevations and latitudes may benefit from higher temperatures because their productivity is predominantly limited by temperature and particularly by a short growing season. In contrast, forests at lower elevations may increasingly suffer from lack of soil water because warming causes an increase in atmospheric water demand, even if precipitation does not decrease (Körner & Paulsen, 2004). Accordingly, a recent research focus has been on warm and dry extremes, whereas cold and wet extreme events (e.g., Figure S1) have received less attention, despite their importance (but see Chen et al., 2019;Vitasse et al., 2019). Moreover, the quantification of forest ecosystem productivity responses to climate extremes at high spatial and temporal resolution is still rare.
A challenge in this context is the synthesis of various and often heterogeneous datasets into a product that summarizes our best knowledge about the dynamics and climate sensitivities of forest ecosystems and that can be used for projections. DVMs can achieve this purpose, especially when coupled with various types of monitoring data (Hartig et al., 2012). We refer to this process as data assimilation (also "modeldata fusion" or "inverse modeling"; see Keenan, Davidson, Moffat, Munger, & Richardson, 2012). Data assimilation can help to better estimate the true ecosystem state, its dynamics, and the associated uncertainties (Keenan, Carbone, Reichstein, & Richardson, 2011;Lahoz, Khattatov, & Menard, 2010;Niu et al., 2014). Data assimilation can also reduce uncertainties in many areas of the modeling process, for example, via initial state updating (data assimilation in a narrow sense), parameter estimation (model-data fusion, model calibration), input updating, and error correction (Houser, De Lannoy, & Walker, 2010).
Despite the fact that this method allows us to combine multiple data sources and types, most studies have focused on the local scale. Hence, an important step forward is now to use large and diverse datasets in combination with DVMs at the regional scale (Cailleret, Bircher, Hartig, Hülsmann, & Bugmann, 2019;Fer et al., 2018;Minunno, Peltoniemi, et al., 2019;Thomas et al., 2017;Van Oijen et al., 2013).
We assimilated extensive and long-term forest ecosystem monitoring data into the 3-PG forest ecosystem model (Landsberg & Waring, 1997). With the parameterized model, we assessed how forest productivity responds to climate extremes across environmental gradients in Switzerland. Switzerland is a highly suitable case study for this purpose, because its elevational gradients span a range of bioclimatic conditions that are comparable to at least 1,800 km of latitudinal gradient in Europe (Halbritter, Alexander, Edwards, & Billeter, 2013), but within a small geographic area. This alleviates the need to control for different synoptic drivers, continentality, population genetic differences, etc. To constrain the parameter distributions of 3-PG and estimate their uncertainty ranges, we compiled monitoring data for two dominant European species (Picea abies (L.) H. Karst. and Fagus sylvatica L.) from 271 sites, totaling almost 800 observation years.
We then used the constrained model parameter distributions to test for shifts in forest productivity responses to climate extremes across environmental gradients. Specifically, we addressed the following questions: (a) What is the contrast in climate response at low versus Grant/Award Number: POIR.04.04.00-00-5F85/18-00; Swiss National Science Foundation, Grant/Award Number: 20FI20_173691 high elevation and in average versus extreme years? (b) How strong are NPP anomalies during warm versus cold extremes and what is the spatial extent of the affected area? Answering these questions helps us to better understand and anticipate possible trajectories of forest ecosystem productivity in a warmer and more variable future climate.

| Monitoring data
We used data from 271 permanent forest monitoring plots covering the actual habitat of P. abies (N = 237) and F. sylvatica (N = 34; Figure 1 (Forrester, Nitzsche, & Schmid, 2019), the Long-term Forest Ecosystem Research Network (LWF; Etzold, Waldner, Thimonier, Schmitt, & Dobbertin, 2014;Schaub, Dobbertin, Kräuchi, & Dobbertin, 2011;Thimonier et al., 2010), and one forest site from the Swiss FluxNet (Etzold et al., 2011;Zielis et al., 2014). We used eight variables that describe stand stocks and characteristics: stem biomass (SB), foliage biomass, root biomass, number of trees, average diameter at breast height (1.3 m; DBH), basal area (BA), leaf area index (LAI), and gross primary production (GPP). To calculate the stand-level stocks, we applied the biomass equations developed for European forests following Forrester, Tachauer, et al. (2017) for each measured tree, and summed it up to the stand level in Mg dry matter/ha. The first observations on each monitoring plot were used to initialize the 3-PG model runs (see below).

| National Forest Inventory
The Swiss NFI records the current state of forests on a regular grid of 1.4 km covering about 6,500 permanent monitoring plots that have been measured since 1983 (Brändli, 2010;Fischer & Traub, 2019). Each plot is remeasured every 10 years, with a one-time change in timing due to a switch from a periodic to a continuous survey in the fourth NFI phase (i.e., since 2009). The NFI plot design comprises nested circular plots, such that every tree with a DBH ≥ 12 cm is recorded within an inner 200 m 2 circle (radius = 7.98 m), and every tree with a DBH ≥ 36 cm is recorded within a 500 m 2 circle (radius = 12.62 m). For every individual tree, the position, DBH, status, and species are recorded. In addition, tree height (H) and crown length (H c ) are measured on a subset of trees.
Age is estimated based on a regression model that was fit to the data obtained either from counting tree rings or counting layers of whorling branches (for the P. abies trees) directly on the plot (Brassel & Lischke, 2001). Management and mortality on each individual monitoring plot were derived from inventory data. The specific year of management interventions and the timing of tree mortality were randomly assigned between two consecutive inventories, as the exact dates are unknown. The monitoring plots for this study were selected using the following criteria: (a) monospecific even-aged stands of either P. abies or F. sylvatica, (b) at least two consecutive remeasurements were available, (c) no ingrowth during the selected period, (d) no obvious measurement errors or missing measurement, and (e) stand age estimation was available for the first observation used. Based on these criteria, we retained 176 (P. abies N = 147, F. sylvatica N = 29) NFI plots, in total accounting for 451 observation years. The time span between the first and the last measurement ranged between 4 and 35 years.

| Experimental Forest Management
The EFM project has been collecting growth and yield data for more than a century (Forrester et al., 2019). The EFM network currently includes 459 permanent monitoring plots, which are measured every 5-12 years, depending on their growth rates, stand age, and research objectives. The EFM monitoring plots are of varying size with precisely defined boundaries, within which all individual trees with a DBH ≥ 8 cm are measured. For each tree, the position, DBH, status, and species are recorded. In addition, H and H c are measured for a subset of trees. Age is estimated based on the planting date in even-aged stands. Management (thinning type and timing) is recorded for each individual monitoring plot and is done in the same year as the measurements. The monitoring plots for this study were selected based on the same criteria used for NFI monitoring plots.

F I G U R E 1
Location of the 271 monitoring plots (dots) distributed across the potential habitats (colored background) of Picea abies and Fagus sylvatica dominated forests in Switzerland (a). Potential habitats are based on the MoGLI projections (Wüest et al., 2020). (b, c) Distribution of the selected plots (colored dots and contours) compared to all potential habitats in Switzerland (gray dots and contours) along the annual mean temperature and annual precipitation sum gradients Based on those criteria, 94 plots remained in our analysis (P. abies N = 89, F. sylvatica N = 5), in total accounting for 331 observation years. The time span between the first and the last measurement ranged between 5 and 30 years.

| Long-term Forest Ecosystem Research and Swiss FluxNet
In the LWF network, information on tree growth and crown condition as well as on the nutrient cycle and the ecosystem water balance is collected to assess the impact of environmental changes on forest functioning Schaub et al., 2011;Thimonier et al., 2010). The LWF network includes 19 permanent monitoring plots, on which monitoring data have been recorded every 1-5 years, since 1994. The LWF uses a standardized protocol for data collection based on the International Co-operative Programme on Assessment and Monitoring of Air Pollution Effects on Forests (UNECE ICP FOrests Programme Co-ordinating Centre, 2016). The monitoring plots for this study were selected based on the same criteria used for the NFI and EFM plots and one plot remained (P. abies N = 1). The retained LWF P. abies plot (Davos CH-Dav) is also part of the Swiss FluxNet ecosystem-scale CO 2 and H 2 O vapor eddy-covariance flux measurement network (Etzold et al., 2011), providing measurements of net ecosystem production, GPP, and ecosystem respiration.

| Dynamic vegetation model
3-PG is a process-based forest ecosystem model that consists of five submodels in a causal chain, starting with light absorption and assimilation, and ending with the conversion of biomass into output variables (Landsberg & Waring, 1997;Sands & Landsberg, 2002).
A simple structure, readily obtainable input data, and a low number of parameters have facilitated the widespread use of 3-PG in various forest types around the world (Gupta & Sharma, 2019). Initially developed for simulating evergreen, even-aged, monospecific forests, the model has recently been further developed for deciduous, uneven-aged, and mixed-species forests (Forrester & Tang, 2016). It is a cohort-based, non-spatially explicit model with a monthly time step.
Each cohort can be a different species and/or age class. Stand-level calculations avoid the propagation of potential errors when scaling up from higher resolution calculations (e.g., leaves or trees) while providing outputs at the level required for this study (Landsberg & Waring, 1997;Pretzsch, Forrester, & Rötzer, 2015).
The first of the five submodels predicts light absorption and GPP using a species-specific canopy quantum efficiency (α C ). The α C is reduced in response to limitations imposed by temperature, frost, vapor pressure deficit (VPD), soil moisture, soil nutrient status, atmospheric CO 2 , and stand age (Almeida, Landsberg, & Sands, 2004;Landsberg & Waring, 1997;Sands & Landsberg, 2002). NPP is calculated as a fixed fraction of GPP (Waring, Landsberg, & Williams, 1998) and is distributed to roots, stems, and foliage by the second submodel. Partitioning to aboveground versus belowground biomass depends on soil nutrient status, VPD, and soil moisture, while partitioning between stems and foliage depends on tree size, with larger trees partitioning a lower proportion of NPP to foliage compared to smaller trees (Landsberg & Waring, 1997;Sands & Landsberg, 2002). The third submodel simulates densitydependent mortality, which is calculated using the −3/2 selfthinning law by Yoda (1963). In the fourth submodel, which calculates the water balance, the Penman-Monteith equation is used to calculate transpiration and soil evaporation, which are added to canopy interception to predict evapotranspiration. Canopy conductance g c is determined using a species-specific maximum g c , LAI and limitations caused by VPD, soil moisture, atmospheric CO 2 , and stand age.
Changes in soil water storage are calculated as the difference between evapotranspiration and rainfall; any excess of the maximum soil water holding capacity is drained off (Sands & Landsberg, 2002).
The fifth submodel converts biomass into output variables such as mean tree diameter, height, BA, wood volume, etc., using allometric relationships. Different management strategies are specified using the residual stocking (trees/ha) after thinning at a nominal age.
Thinning from below or above is achieved by specifying the fraction of the foliage, root, and SB of an average tree that was thinned (Landsberg, Mäkelä, Sievänen, & Kukkola, 2005). All submodels were evaluated by comparing predictions of the given process against empirical data of that process for many different forest types (Gupta & Sharma, 2019;Landsberg & Sands, 2011), including central European forests (Forrester, Ammer, et al., 2017;Nolè et al., 2009).
For our simulations, we used a re-implementation of the 3-PG model programmed in Fortran 90 (Minunno, Hartig, & Trotsiuk, 2019). It was driven with time series monthly mean of daily minimum and maximum temperatures (T min , T max , °C), rainfall (P rcp , mm/month), monthly mean of daily solar radiation (S rad , MJ m −2 day −1 ), and the number of frost days (F days , days/month with T min < 0°C). We used spatially interpolated monthly meteorological data as input for 3-PG ( Figure S2). The interpolation (100 m spatial resolution) of the meteorological data was done by the Landscape Dynamics group (WSL, Switzerland) using data from MeteoSwiss stations (Swiss Federal Office of Meteorology and Climatology) by employing the DAYMET method (Thornton, Running, & White, 1997). Site-specific information on soil type and plant available soil water was retrieved from a digitized soil suitability map of Switzerland (scale 1:200,000 ;Frei et al., 1980;Swiss Federal Statistical Office, 2000).

| Parameter estimation
We used Bayesian inference to derive the parameter estimates and uncertainties of the 3-PG model. The approach accounts for observational uncertainties and to make use of multiple types of data at different temporal scales. We assumed uniform (i.e., non-informative) prior distributions for each of the 54 model parameters. The ranges of the priors (Tables S3 and S4) were set to the minimum (maximum) value found in the literature minus (or plus) half of the range for this parameter (following Augustynczik et al., 2017). The likelihood function was constructed to be robust against outliers by modeling the residual error as a Student's t distribution with sampled degrees of freedom (see Code S1; Lange, Little, & Taylor, 1989). We used the Differential Evolution Markov Chain Monte-Carlo algorithm (DEzs MCMC, ter Braak & Vrugt, 2008), implemented in the BayesianTools R package (Hartig, Minunno, & Paul, 2019) to estimate the joint posterior distribution for the model parameters. For each species, we ran three independent DEzs MCMC runs, each with three internal chains, and tested convergence by visual inspection of the trace plots and additionally using the Gelman-Rubin diagnostic (Gelman & Rubin, 1992), with convergence being accepted when the multivariate potential scale reduction factor was ≤1.1. Three independent DEzs MCMC chains with 2.4 × 10 7 (P. abies) and 1.7 × 10 7 (F. sylvatica) iterations were required to achieve convergence. All analyses and calculations were performed in the R language for statistical computing (R Core Team, 2018).

| Model evaluation and validation
To evaluate the skill of the model and generate model projections, we calculated posterior predictive distributions by running the model with 1,000 random samples from the parameters' posterior distribution.
Model performance was evaluated using the percentage bias (pBias), root mean squared error (RMSE), and normalized root mean squared error (NRMSE). We first calculated statistics on the plot level, and then averaged over plots for each of the 1,000 samples. For the validation, we only used the most recent set of observations at all permanent monitoring plots to maximize the time between initialization and validation, which ranged from 4 to 35 years. To perform cross-validation, we randomly split the full set of monitoring data into two equally sized groups, resulting in a calibration and a validation set.

| Model simulations
We simulated forest productivity (i.e., NPP) for the species' potential distribution range in Switzerland (Wüest, Bergamini, Bollmann, & Baltensweiler, 2020) on a 1 × 1 km grid for a total of 10,100 grid points for P. abies and 7,030 grid points for F. sylvatica. For this purpose, we first simulated the growth of P. abies and F. sylvatica monocultures with the average climate observed during the 1961-1990 period, until the age of 40 years (spin-up). The stands were simulated starting as 2-year-old plantations with an initial density of 10,000 trees/ha.

| Parameter estimation
The width of the posterior (measured by the 95% quantile range) was on average 59% (P. abies) and 32% (F. sylvatica) smaller than the prior range across all parameters (Tables S3 and S4). The largest reduction in uncertainty for both species was for parameters associated with allometric relationships, biomass partitioning, and stem mortality.
Monitoring data were least informative for the parameters associated with branch and bark fractions and soil fertility. The parameter u, controlling the number of degrees of freedom in the Student's t likelihood (code S1), was much smaller for the F. sylvatica compared to the P. abies monitoring plots, indicating heavier tails in the error for P. abies (Tables S3 and S4).

F I G U R E 2 Statistics on predictive error (percent bias [pBias]
, normalized root mean squared error [NRMSE], and root mean squared error [RMSE]) of the 3-PG model. The posterior predictive uncertainty was calculated by drawing 1,000 parameter combinations from the posterior distribution and calculating model predictions for these combinations. The dots represent the median value of the posterior predictive distribution, while the horizontal lines represent the 95% confidence interval. DBH, diameter at breast height Predictions based on the posterior distributions significantly improved compared to predictions based on the prior distributions ( Figure S3) for both P. abies and F. sylvatica (Figure 2). The NRMSE were below 8%, and the magnitude of the pBias was below 10% for all variables, while it reached up to 600% with the prior distributions ( Figure S3). The correlations between observed and simulated values were high for all variables with r 2 ≥ .90 for P. abies and r 2 ≥ .87 for F. sylvatica. The RMSE for the change in stem dry biomass (ΔWS) of P. abies and F. sylvatica was 15 and 18 Mg/ha, respectively, while pBias was −7% and −9%, respectively (Figure 2). The cross-validation based on 50% of all plots confirmed the high accuracy of the model ( Figure S4).

| Simulations of net primary productivity at the country scale
Annual mean net primary productivity (NPP) simulated on the 1 × 1 km grid for the described hypothetical stands at the age of 41-70 years within the species distribution range across Switzerland during the 1991-2018 period was 5.4 ± 1.5 Mg C ha −1 year −1 (mean ± standard deviation of the mean) for P. abies and 5.3 ± 1.0 Mg C ha −1 year −1 for F. sylvatica (Figure 3). There was a strong negative correlation between annual NPP and elevation (p < .001), with an average decrease of 2.86 ± 0.006 Mg C ha −1 year −1 km −1 for P. abies and 0.93 ± 0.010 Mg C ha −1 year −1 km −1 for F. sylvatica. On average, P. abies showed higher NPP (5.9 ± 4.1%) during the recent warmer period  compared to the reference period , while for F. sylvatica, the change was not significant. There was strong agreement in terms of the trend and the magnitude of NPP simulated by the 3-PG model and NPP derived from MODIS ( Figure S5) for P. abies, but less so for F. sylvatica.
The calibrated 3-PG model indicates that annual NPP of P. abies and F. sylvatica was considerably reduced during extreme years ( Figure 4). During the warm-dry year of 2018 ( Figure S1), NPP was strongly reduced (anomaly below −25%) for one-fifth of the potential habitat area in Switzerland (P. abies: 21%; F. sylvatica: 15%; Figure 5). Interestingly, the predictions for P. abies showed F I G U R E 3 Trajectory of net primary productivity (NPP) along the elevational gradient simulated by the 3-PG model for Picea abies (green) and Fagus sylvatica (orange) potential distribution ranges. The respective solid lines represent the average for the 1991-2018 period, the dashed and dotted lines represent 50% and 95% confidence interval, respectively. The shaded areas represent the density distribution of the potential species habitat along the elevational gradient F I G U R E 4 Trajectory of simulated net primary productivity (NPP) anomalies (percentage deviation) in selected extreme years ( Figure S1) relative to the 1961-1990 reference period for Picea abies (green) and Fagus sylvatica (orange). The respective solid lines represent the average for the 1991-2018 period, the dashed and dotted lines represent 50% and 95% confidence interval, respectively. NPP anomalies were calculated for each grid cell of the potential species distribution ranges

| D ISCUSS I ON
Estimating the impact of climate extremes on forest ecosystem productivity is essential for understanding their role for regulating the regional carbon cycle and its drivers. Many previous studies have examined the effect of climate extremes on forests focus on extremely warm-dry years. We stress here the importance to also account for cold extremes, even though these might become less likely under climate change. By assimilating observations from 271 permanent long-term forest monitoring plots with the 3-PG forest ecosystem model, we were able to quantify the spatiotemporal changes in forest ecosystem productivity in response to climate extremes at the country scale, highlighting that not only extremely warm and dry years, but also extremely cold and/or wet years significantly impact forest NPP. Our results further indicated a high altitudinal and spatial variation in forest productivity response to climate extremes, which provides important information on forest vulnerability across the species' range.

| Parameter estimation
So far, only few studies have assimilated extensive forest monitoring datasets into a DVM through techniques of parameter estimation (but see Cailleret et al., 2019;Fer et al., 2018;Minunno, Peltoniemi, et al., 2019;Thomas et al., 2017), even though recommended by several authors to improve large-scale model projections Hartig et al., 2012). Our study demonstrates that it is possible to integrate monitoring data from multiple networks across a wide bioclimatic gradient into a process-based forest ecosystem model 3-PG. The resulting uncertainty in the parameter estimates was relatively low (Tables S3 and S4) and comparable to other studies that calibrated the 3-PG model (Augustynczik et al., 2017;Thomas et al., 2017). Not surprisingly, the monitoring data were most informative for constraining parameters that are directly related to stand structure. However, the calibration also reduced parametric uncertainty in parameters not directly related to stand structure. For example, the maximum a posteriori estimates for the parameter LAIgcx (the LAI at which leaf area was not limiting transpiration) were reduced by 19% for P. abies and 5% for F. sylvatica, toward values consistent with empirical observations (Schulze, Kelliher, Korner, Lloyd, & Leuning, 1994). We conjecture that the lower reduction of parametric uncertainty in this parameter is both due to its lower influence on stand structure (to which the model was calibrated), but possibly also due to higher intraspecific variability, that is, the parameter values are not identical for all sites. To test this in future studies, spatially variable parameterization could be considered (cf. Vanderwel, Rozendaal, & Evans, 2017).
The relatively low predictive error of the calibrated model (e.g., pBias ≤ 9%) supports the use of Bayesian inference to estimate model parameters and their uncertainty. The observed reduction of the predictive error is comparable to that achieved in other studies in European forest landscapes (Minunno, Peltoniemi, et al., 2019;Van Oijen et al., 2013). The lower predictive error for P. abies compared to F. sylvatica (Figure 2) is likely due to the larger number of P. abies (N = 237) monitoring plots compared to F. sylvatica (N = 34) plots, especially due to the higher number of EFM plots (P. abies: 89; F. sylvatica: 5). In previous studies, data from permanent monitoring plots (like EFM) were shown to be more useful for model calibration than data from forest inventories (Minunno, Peltoniemi, et al., 2019;Van Oijen et al., 2013). Minunno, Peltoniemi, et al. (2019) argued that a main problem of the NFI data in their study was its shorter time span compared to their EFM data. The NFI and EFM data used by us have comparable time span, and might equally contribute to model calibration.
Even though our approach strongly improved model performance, we recognize and acknowledge some limitations. The reliability of initial conditions, climatic forcing data, and monitoring data are important for an ecologically meaningful data assimilation process (Van Oijen et al., 2013). Site nutrient status and available soil water are key variables in the 3-PG model, but accurate in situ measurements are rare for a large number of plots. The soil suitability map that we used is rather generic and actual values for specific monitoring plots may differ substantially from the mapped data.
Moreover, as for the vast majority of other DGVMs, we assumed that species-specific parameters are identical across their range, despite ample evidence for intraspecific variability of functional traits within species (Moran, Hartig, & Bell, 2016). Such a strong assumption simplifies the calibration process, but may also lead to inaccurate predictions about climate responses and forest resilience (see also Berzaghi et al., 2019). Thus, it will be beneficial to apply spatially variable parameterizations, despite the substantial computational cost. Thomas et al. (2017) successfully applied such an approach (Vanderwel et al., 2017) to constrain soil nutrient status in the 3-PG model based on site index and mean annual temperature.
Finally, the question of model data assimilation is closely connected to structural model error. For example, in 3-PG, the ratio between NPP and GPP is a constant, irrespective of environmental conditions (Amthor, 2000;DeLucia, Drake, Thomas, & Gonzalez-Meler, 2007). Empirical studies, however, show that the GPP/NPP ratio can vary considerably (Collalti & Prentice, 2019;Zhang, Xu, Chen, & Adams, 2009). Another example is that inter-annual variability in reproduction, not considered in 3-PG, can substantially impact carbon use and forest growth, especially for masting species such as F. sylvatica (Hacket-Pain et al., 2018). These and other limitations and uncertainties leave room for further improvement of 3-PG and DVMs in general.

| Simulations of NPP at the country scale
Simulated productivity for both species continuously decreased along the elevational gradient. This is consistent with results from other empirical studies in which forest productivity positively responds to temperature (Babst et al., 2013;Luyssaert et al., 2007).
We also found strong agreement in the NPP-elevation relationship between 3-PG and MODIS derived data for P. abies, but not for F. sylvatica ( Figure S5). We think the latter is most likely due to the fact that the presence of F. sylvatica is firstly less accurately mapped than for P. abies, and secondly, that F. sylvatica often occurs in mixtures, making comparisons to MODIS data less reliable.
Simulated NPP across species distribution ranges was primarily controlled by temperature, soil water, and VPD ( Figure S6).
Environmental constraints due to low temperatures increased nonlinearly with increasing elevation, reflecting the temperature control of photosynthesis. Environmental constraints due to reduced soil water availability and increased VPD decreased with increasing elevation, but at a lower magnitude. Accordingly, optimal conditions for both P. abies and F. sylvatica were found at the lowest elevations.
Surprisingly, the decrease in NPP along the elevational gradient for F. sylvatica is rather small compared to other studies (e.g., Zianis & Mencuccini, 2005). We hypothesize that this could be due to an incomplete coverage of the upper edge of the F. sylvatica distribution.
Our model simulated extreme years to cause substantial decreases in forest NPP along the Swiss elevational and bioclimatic gradient, which is in line with previous studies that assessed the impact of climate extremes on forest productivity during extremely warm-dry years (Kannenberg et al., 2019;Vitali, Büntgen, & Bauhus, 2017;Vitasse et al., 2019). Additionally, our study highlights the importance of accounting for both extreme cold and/or wet years. For both species, an extremely cold growing season can have a strong negative impact on NPP that is similar in magnitude to that from extreme warm-dry years ( Figure 5). Similarly, Vitasse et al. (2019) found that late frosts can impact F. sylvatica growth in a magnitude comparable to extreme drought. However, it is also important to mention that 3-PG and DGVMs in general have often overestimated the sensitivity of forest to drought, compared to observations (Klesse et al., 2018;Thomas et al., 2017).
Consistent with empirical studies, our analysis suggests that extreme years caused divergent forest growth responses along the Swiss elevational and spatial gradients (Hartl-Meier, Dittmar, Zang, & Rothe, 2014;Jolly et al., 2005;Vitali et al., 2017). An increase in temperature can enhance growth at higher elevations but lead to drought-induced growth decline at lower elevations, particularly for trees growing under high levels of competition Jolly et al., 2005;Primicia et al., 2015;Schurman et al., 2019).
As an example, the extremely warm-dry year of 2003 promoted better growing conditions at higher elevations for P. abies. At lower elevations, NPP decreased due to reduced soil water availability (caused by increasing evaporative demand) during the growing season ( Figure S6). The year 2003 had a different pattern in precipitation than 2018, with smaller negative precipitation anomalies, especially at higher elevations ( Figures S1 and S2). Thus, the abrupt decrease in NPP at lower elevations was compensated by an increase in NPP at higher elevations, which was not the case in 2018.
For F. sylvatica, the substantial reduction in NPP during both extremely warm-dry years (2003 and 2018) was consistent with the patterns along the elevational and spatial gradients (Figures 4 and   5). The difference in precipitation anomalies between 2003 and 2018 occurred only above 1,200 m a.s.l., which is above the simulated range of F. sylvatica.
Because of increasing frequency and intensity of warm-dry events due to climate change, our results suggest that P. abies and F. sylvatica will show a substantial reduction in NPP at the lower elevational band, up to 800 m a.s.l. This effect is exacerbated by the fact that the drought response along climatic gradients will likely be altered in a nonlinear way (Kannenberg et al., 2019). However, the impact of drought on tree performance and forest productivity strongly depends on its seasonal timing (Crimmins, Crimmins, Gerst, Rosemartin, & Weltzin, 2017;Vicente-Serrano et al., 2013), calling for more research on intra-annual tree growth and climate sensitivity. Still, significant reduction in NPP on a large area (up to 21% for P. abies and 15% for F. sylvatica) in drought years provides incentive to reconsider the current forest management strategies and favor more drought tolerant genotypes of present tree species (Fréjaville, Fady, Kremer, Ducousso, & Garzón, 2019), or alternative species at the lower elevations. Furthermore, reducing low-temperature constraints without necessarily reducing the probability of damaging late frosts (due to advanced phenology; Ma, Huang, Hänninen, & Berninger, 2019), our results suggest that F. sylvatica may experience a stronger reduction in NPP and potentially increased mortality in the future.

| CON CLUS IONS
We assimilated an extensive collection of data from 271 permanent monitoring sites into the 3-PG forest ecosystem model, and then simulated the climate sensitivity of two dominant European tree species across Switzerland. For the first time, it could be shown at a high spatial resolution that climate extremes impact forest productivity in more complex ways than simply shifting the response upwards in elevation. Our model suggests on the one hand that, during extremely warm-dry years, forests at lower elevations will suffer from soil water deficit and increased evaporative water demand, whereas forests at higher elevations, where trees are growing still below their temperature optimum, will benefit from warmer conditions. On the other hand, for trees growing on intermediate and less extreme sites, the model indicated highly differentiated year by year responses depending on the respective combination of weather forcing. The model-data fusion approach used in this study allowed us to model highly site-specific NPP from long-term monitoring data.
Such robust estimates of NPP are key for increasing our understanding of forests carbon dynamics under climate extremes. During the extremely cold and wet years, both species experienced strong reductions in NPP, which are comparable in magnitude to extremely warm and dry years. Neither of these broad effects, however, are linear or homogenous in space. The nonlinear shifts in NPP during extreme years along the elevational gradients indicate the value and necessity of spatially resolved analyses of the impacts of climate extremes and changes.

ACK N OWLED G EM ENTS
This study was funded by the SwissForestLab based on a proposal by W.E. and co-authors. V.T., W.E., M.S., M.F., and N.B. designed the research. V.T. and F.H. performed the analysis. V.T., F.H., M.C., F.B., D.F., and M.S. wrote the paper with substantial inputs from all coauthors. We are thankful to Dirk Schmatz for providing the gridded DAYMET data, and acknowledge MeteoSwiss for providing meteorological station data for the spatial interpolation. We acknowledge WSL and ETH and its scientists, field staff, laboratory personal, and database managers who designed, carried out, and maintained the measurements on the permanent monitoring plots used in this study. We also are grateful to Samuel Abiven, Susanne Burri, Frank The data that support the findings of this study are available on request from the Swiss NFI (https ://www.lfi.ch), Swiss Longterm Forest Ecosystem Research LWF (www.lwf.ch), and Swiss FluxNet (http://www.swiss fluxn et.ch). For our simulations, we used a reimplementation of the 3-PG model programmed in Fortran 90 (Minunno, Hartig, & Trotsiuk, 2019).