Nitrogen deposition shows no consistent negative nor positive effect on the response of forest productivity to drought across European FLUXNET forest sites

Atmospheric reactive nitrogen (N) deposition is an important driver of carbon (C) sequestration in forest ecosystems. Previous studies have focused on N-C interactions in various ecosystems; however, relatively little is known about the impact of N deposition on ecosystem C cycling during climate extremes such as droughts. With the occurrence and severity of droughts likely to be exacerbated by climate change, N deposition—drought interactions remain one of the key uncertainties in process-based models to date. This study aims to contribute to the understanding of N deposition-drought dynamics on gross primary production (GPP) in European forest ecosystems. To do so, different soil water availability indicators (Standardized Precipitation Evapotranspiration Index (SPEI), soil volumetric water) and GPP measurements from European FLUXNET forest sites were used to quantify the response of forest GPP to drought. The computed drought responses of the forest GPP to drought were linked to modelled N deposition estimates for varying edaphic, physiological, and climatic conditions. Our result showed a differential response of forest ecosystems to the drought indicators. Although all FLUXNET forest sites showed a coherent dependence of GPP on N deposition, no consistent or significant N deposition effect on the response of forest GPP to drought could be isolated. The mean response of forest GPP to drought could be predicted for forests with Pinus trees as dominant species (R 2 = 0.85, RMSE = 8.1). After extracting the influence of the most prominent parameters (mean annual temperature and precipitation, forest age), however, the variability remained too large to significantly substantiate hypothesized N deposition effects. These results suggest that, while N deposition clearly affects forest productivity, N deposition is not a major nor consistent driver of forest productivity responses to drought in European forest ecosystems.


Introduction
Terrestrial ecosystems have the potential to take up significant amounts of carbon dioxide from the atmosphere through photosynthesis and growth and thereby mitigate climate change. As the sequestration of carbon (C) in terrestrial ecosystems predominantly occurs in forest ecosystems (Pan et al 2011), forests largely drive the terrestrial C balance. The compound effect of many interacting drivers determines whether a forest is a net sink of C (i.e., taking up C from the atmosphere over multi-year timescales) or a net source of C (i.e., releasing C to the atmosphere over multi-year timescales). These drivers include edaphic and climatic factors such as soil nutrient and moisture conditions and air temperature and humidity. Drought, fires and outbreaks of insect herbivores and fungal pathogens can reduce forest productivity for many years and can also cause widespread forest mortality (e.g., Anderegg et al 2020). On the other hand, some factors can potentially increase ecosystem carbon storage, such as CO 2 fertilization, ozone exposure, nitrogen (N) deposition and forest management.
Earth system models (ESMs) can be used to quantify and predict ecosystem responses to a changing climate and the feedbacks involved. However, ESMs are known to involve large uncertainties in terrestrial C feedbacks (Friedlingstein et al 2014). These uncertainties partially result from a lack of knowledge of the physical and biogeochemical processes responsible for these C cycle feedbacks. Furthermore, we know relatively little about ecosystem responses to multiple, simultaneous stressors and their interactions, as most studies to date focus on the effect of one single stressor on plant growth. The impact of co-stressors is highly variable across ecosystems and, currently, ESMs are ill-equipped to model these interactions (Drewniak and Gonzalez-Meler 2017). One key example is the interplay between increased nitrogen (N) deposition and drought. With the frequency and intensity of droughts likely to increase globally as a result of climate change, and N deposition projected to increase to 88 Tg N yr −1 by 2100 under representative concentration pathway (RCP) 8.5 scenario (Lamarque et al 2013), understanding the effects of N deposition-drought interactions on forest productivity becomes increasingly important. N deposition-drought interactions are found to be interdependent and non-additive and remain one of the least quantified processes that vary locally and therewith one of the major uncertainties in ESMs to date (Drewniak and Gonzalez-Meler 2017).
Our objective is to contribute to understanding the variability of combined N deposition-drought interactions on C dynamics in forest ecosystems by (1) quantifying the effect of different types of drought on the gross primary production (GPP) from FLUXNET observations at European forest sites, and (2) linking these responses to different levels of N deposition for different forest types, climate zones and soil types.
1.1. Forest response to N deposition Nitrogen (N) is an important nutrient in ecosystems and critical for driving photosynthesis and growth (Evans 1989, Oren et al 2001, Fernández-Martínez et al 2014. Forest growth in temperate and boreal forest ecosystems is generally limited by N availability, and N deposition is, therefore, an important driver of forest growth in these ecosystems. The three main effects of N deposition on ecosystems are changes in foliar or leaf N, changes in biomass partitioning and increases in biomass N (Bobbink et al 2010). Reactive nitrogen (N) emissions have substantially increased during the last century, causing enhanced atmospheric N deposition on forest ecosystems (Dentener et al 2006, Lamarque et al 2013. With most forest ecosystems being N limited, this has resulted in significantly increased net primary production (NPP) and subsequently carbon (C) sequestration in trees (e.g. (Vries et al 2009, Thomas et al 2010).
Several studies have found increases in foliar N and, as a consequence, decreases in C:N ratios in leaves under elevated N deposition (e.g. (Boggs et al 2005, Elvir et al 2005, Pregitzer et al 2008). Increases in foliar N are often associated with increases in aboveground net primary production (ANPP) (Pregitzer et al 2008) as foliar N is generally found to be linearly and positively related to leaf photosynthetic capacity. Under elevated N availability, the ANPP may also increase because of reductions in N use efficiency and C allocation shifts away from mycorrhizae, leaving more C available for growth (Talhelm et al 2011).
Gains in forest productivity from increased N deposition may lead to an increase in litter production, which may ultimately lead to an increase of N mineralization in soils. Moreover, increased productivity in trees may result in higher C sequestration in soil due to higher soil C inputs by litterfall and reduction of organic matter decomposition (Lu et al, 2011, Janssens et al 2010. The increase in forest productivity is often correlated with an increase in aboveground biomass (e.g. (Pregitzer et al 2008, Quinn Thomas et al, 2010), resulting from changes in C partitioning in trees. Most of this increase in aboveground biomass is allocated to stems (Pregitzer et al 2008;Vries et al 2014), and generally resulting in faster biomass accumulation and taller, thinner trees. However, the increase in growth varies largely among different tree species. Needle-leaved boreal forests, for instance, have higher C/N ratios in all compartments and show a lower C-N response and N use efficiency than deciduous and evergreen broadleaved forests in temperate regions (Vries et al 2014). The C-N responses of forest ecosystems are non-linear and vary with N input level and time. At high and chronic N deposition rates, forest ecosystems will become N saturated. N loss rates by leaching, runoff and denitrification will increase, leading to a decrease in N retention (Vries et al 2014). Forest growth is initially stimulated by low levels of N deposition, as N limitations for growth diminish. At higher N deposition levels, when the ecosystem starts to become N saturated, the stimulating effect may decline and even reverse due to soil acidification and nutrient imbalances (Aber et al 1998).
N deposition can also lead to an increase in forest biomass N. Eventually, this extra N will flow from the canopy into the litter pool, where it can cause faster decomposition (Zhu et al 2015) and enhanced N mineralization to the forest floor. At low external N inputs, the breakdown of organic C due to microbial activity is stimulated, increasing respiration. This stimulating effect diminishes completely, however, after ecosystems start to become N saturated, due to N-induced microbial community and decomposing enzyme shifts, resulting in a reduction of forest soil respiration (Janssens et al 2010).
Modelling studies estimated that approximately 24% of the global historical C sink (between 1900-2006) was driven by N deposition effects (Fleischer et al 2015) and that N deposition accounts for the additional storage of approximately 175 Pg C in forests since pre-industrial times (Bala et al 2013). It is estimated that N deposition currently increases the global forest C sink by around 276 to 448 Tg C yr −1 (Vries et al 2014). Quantification of stimulation of forest growth as the result of nitrogen deposition, is still under discussion (Ehtesham and Bengtson 2017), with estimates of net ecosystem production (NEP) being simulated at rates of 30-75 kg C per kg N down to 16 kg C per kg N (e.g. Butterbach-Bahl et al 2011). A recent meta-analysis showed the difference between old and younger forest with a factor of 5 for the stimulation of aboveground woody production, but also low productivity forests respond more strongly than high productivity forests, (Schulte-Uebbing and Vries 2018).
Excess amounts of N deposition, on the other hand, can also cause nitrate leaching, reductions in forest biodiversity and may ultimately lead to growth reductions by N saturation (Bobbink at al., 2010). N deposition thus significantly influences the response of forest ecosystems, and consequently plays a vital role in understanding the long-term response of forests to the effects of climate change (increased CO 2 levels, elevated temperatures and changes in water availability) (Vries et al 2009).

Forest response to drought
The response of forest ecosystems to drought depends on various factors, including the sensitivity of dominant tree species to drought, soil characteristics -especially those related to soil water retention and rooting depth-, the climatic zone and the severity of the particular drought (e.g. (Schwalm et al 2010;von Buttlar et al 2018)). Generally, two mechanisms are identified through which plants are negatively impacted by drought, carbon starvation and hydraulic failure. Carbon starvation and hydraulic failure can co-occur and both mechanisms generally result in lower C assimilation and may ultimately lead to tree mortality (Sevanto et al 2011). Carbon starvation can occur when leaf stomata close to constrain water losses, also impairing the diffusion of CO 2 into the leaf and thereby limiting C assimilation. Reduced C assimilation results in fewer carbohydrates available for growth and maintenance and may ultimately lead to serious tissue damage if existing C reserves are insufficient to sustain plant maintenance requirements. Hydraulic failure occurs when xylem functioning is partially or completely lost through xylem embolism, which inhibits water and nutrient transport from the roots to the leaves and leads to tissue desiccation (Mcdowell et al 2008, Adams et al 2017. The drought tolerance of trees depends on many morphological and physiological traits and drought response mechanisms, such as stomatal control, hydraulic redistribution, tissue desiccation tolerance or allometric plasticity (Baker et al 2008). To withstand droughts, trees may also reduce C demand for instance by leaf senescence or down-regulation of respiration (Sala et al 2010). The effects forest may experience under droughts include changes in C availability, mobilization and transport, increases in N limitation and changes in biomass partitioning (Drewniak and Gonzalez-Meler 2017).
Stomatal closure or leaf senescence during droughts may lead to lower C availability and a higher N limitation. Stomata closure reduces photosynthesis and may lead to lower litter production and N mineralization in the soil over long periods, which in turn may lead to additional N limitations (Schimel et al 2007). Moreover, low soil moisture levels during droughts reduce nutrient flow and diffusion in soils, resulting in additional nutrient limitations, which may reduce photosynthesis levels even further. As a response to the early stages of drought, trees have shown to increase the root-to-shoot ratio to maintain transpiration. Several studies have found changes in C partitioning in trees following droughts, such as translocation of carbohydrates to roots and increases in root-to-leaf biomass production (e.g. (Hanson et al 2007, Hertel et al 2013). Trees may also alter morphological traits of their roots to help fulfil water demands in response to drought (Meier and Leuschner 2008). Furthermore, certain trees can extend root systems to deep soil layers or to redistribute water through the soil column via hydraulic redistribution (Hanson et al 2007). Hydraulic redistribution can effectively transfer water upwards into dry soil layers, and even move water deeper into the soil to be protected from evaporation or competition. These mechanisms may prevent hydraulic failure in trees by maintaining water potential above the failure limits.
The frequency, severity and timing of a drought plays an important role in the magnitude of a forest productivity response to drought. Over shorter time scales, stomata regulate water loss and can result in a decline of photosynthesis. Over longer time scales, frequent but less severe droughts may on the other hand increase forest drought tolerance as changes in tree physiology may occur that optimize hydraulic conductance (Martin-Stpaul et al 2013). The time scale at which forest growth responds to drought reflects its ability to cope with water deficits and is a proxy for drought vulnerability. The period between water shortage and impact on growth differ among different forest types and climate zones. For example, forests located in semi-arid and sub-humid areas tend to respond over longer time-scales than forests located in humid areas. Some forests may not respond to drought at all, for example, forests located in very cold and humid areas (Vicente-Serrano et al 2014). The timing of drought also plays a key role in forest response. Droughts coinciding with peak growth periods will for instance likely result in a stronger response of forest C uptake and higher tree mortality.
The response to drought varies strongly among different tree species. For example, in angiosperms, European beech (Fagus sylvatica) is generally found to be more vulnerable to drought-induced growth reductions compared to European oak (Quercus robur) (van der Werf et al 2007, Scharnweber et al, 2011). In gymnosperms, Norway spruce (Picea abies) is found to be more drought vulnerable in terms of radial growth compared to black pine (Pinus nigra) and Douglas fir (Pseudotsuga menziesii) due to its relatively shallow rooting depth (Lévesque et al 2014). Pine species (Pinus spp.), for example, Scots pine (Pinus sylvestris) and the Maritime pine (Pinus pinaster) are generally considered more drought-tolerant, although the Maritime pine shows a strong stomatal response (decline in stomatal conductance and photosynthesis) to drought (Picon et al 1996). In addition to species-specific drought tolerance, angiosperms and gymnosperms also have fundamental differences in their drought response strategies. A recent study, for instance, found that angiosperms, initially, have lower resistance to droughts, while gymnosperms generally show reduced recovery after droughts (Desoto et al 2020).

N deposition and drought interactions
The interactions between N deposition and drought are not always straightforward and vary with climate zone, forest type, as well as drought severity and duration. Furthermore, N deposition and drought may have counteracting effects on forest growth on different time scales. Where N deposition tends to increase photosynthetic capacity, photosynthesis is generally limited by drought. As a result, drought could negate increases in forest productivity resulting from increased N deposition (Wang et al 2012, Lui et al 2013. However, experiments have also shown that N addition may also partially alleviate drought impacts on growth (Wang et al 2012). N deposition and droughts both affect nutrient availability in forest ecosystems. Reduced soil moisture content during drought may lead to a decrease in organic matter decomposition, and immobilize nutrients in the soil. However, extra N available from N deposition may help reduce N limitations (Hanson and Weltzin 2000). Very low N deposition rates may also lead to a lack of N reserves, impairing the trees ability to sustain drought stress (Gessler et al 2017). Excessive N deposition, on the other hand, may result in nutrients imbalance within the soil and plants tissue, and cause a reduction of available cations (Mg, Ca and K), that play an important role in physiological drought defense mechanisms (Schulze 1989). Furthermore, excess nitrogen leads to an increase in the shoot-root ratio and a shallower rooting system and decreasing fine root biomass increasing the risk of drought (e.g., Li et al 2015). Therefore, we would expect an optimum value of N deposition where N deposition alleviates drought stress.
The anatomical and physiological traits of trees are also influenced by N deposition and drought. Chronic, elevated N deposition levels may lead to lower root to shoot ratios, a shallower rooting system and decreasing fine root biomass increasing the risk of drought (e.g., Li et al 2015). Reduction in root biomass under combined N deposition and drought appear to be more severe for younger trees than for older trees (Palátová et al 2002), as trees tend to allocate more biomass to roots as they age (Meyer-Grünefeldt et al 2015). Trees that allocated more C to stems under elevated N deposition may, on the other hand, initially experience less water stress during droughts due to extra water storage (Albuquerque et al 2013).

Paper setup
A global network of eddy-covariance measurements allows us to look at CO 2 fluxes between the atmosphere and biosphere on an ecosystem level. The FLUXNET2015 dataset provides an estimate of the net CO 2 balance of an ecosystem, as well as partitioning into upward-and downward CO 2 fluxes. Here we focus on the European forest sites because of the availability of more detailed nitrogen deposition data. One important component of the ecosystem C cycle is the net ecosystem exchange (NEE), which is the net CO 2 flux from the ecosystem to the atmosphere (Chapin et al 2006). The NEE corresponds to the difference of photosynthetic C uptake, or the gross primary production (GPP), and the total ecosystem respiration (R eco ), which includes both autotrophic and heterotrophic respiration (Papale et al 2006). So far, accurate estimates of nitrogen deposition are lacking for these sites. Flechard et al (2020) used the inferential method and the EMEP model with a coarse grid to provide estimates of nitrogen deposition. Here we used the recent updated LOTOS-EUROS model to provide sitespecific estimates using more detailed local estimates. Furthermore, we use multiple methods to quantify the response of GPP at European forest eddy-covariance sites to drought. For this, several drought indices and soil water availability products (Standardized Precipitation Evapotranspiration Index, soil volumetric water layer) were used, indicating different types and durations of droughts. The computed responses of forest GPP to drought were then linked to the modelled N deposition estimates to assess potential N-drought interactions. In doing so, the across-site variability in circumstances, such as different dominant forest type, soil types and climate zones is taken into account.
We hypothesize that the response of forest GPP to droughts is smallest at an intermediate level of N deposition of max 10 kg ha −1 yr −1 , which is generally accepted as the critical load for nitrogen deposition (Bobbink et al 2010). N deposition during drought may not always be detrimental and alleviate some of the impacts of drought. Due to a potential shortage of available nutrients, we expect the response of forest GPP to drought to be more severe in forests that experience very low levels of N deposition than in forest with intermediate levels of N deposition. At high N deposition levels, on the other hand, we expect that biomass partitioning and the anatomical and physiological development in trees plays a more dominant role in drought impact. As forest ecosystems under chronic, elevated N deposition levels may have reduced amounts of root biomass, we expect a more severe response of forest GPP to drought, especially for longer-lasting droughts. In the next section of the paper, a description of the used datasets and model is given. The methodology is described in section 3 and the results are presented in section 4. Finally, the results are discussed in section 5.

Datasets and model 2.1. FLUXNET2015 data
The FLUXNET2015 dataset provides ecosystem-scale data on CO 2 , water, and energy exchange between the biosphere and the atmosphere, as well as other meteorological and biological measurements (Pastorello et al 2020, Gilberto et al 2020). The eddy covariance technique measures land-atmosphere energy and greenhouse gas fluxes at an ecosystem level at a high temporal resolution of 30 min. High frequency (10-20 Hz) measurements of the vertical wind velocity and scalar variables (e.g. CO 2 , temperature) are used to provide an estimate of the net exchange of that scalar variable over a footprint area around the point of measurement (Aubinet et al 2012). The measured NEE is partitioned into the GPP and R eco using the night-time method (Reichstein et al 2005).
Here, we use estimates of the daily and monthly GPP and R eco , as well as some other meteorological observations.
We selected forest sites in Europe with at least 6 years of observations left after the application of a preprocessing filter (see Methods section). An overview of the used sites including auxiliary information is given in table S1 and S2 (available online at stacks.iop.org/ERC/3/125003/mmedia) in the supplement. The locations of the selected sites are shown in figure 1.

Soil type
The physical properties of the soils on which the forests are growing are obtained from the European Soil Data Centre (ESDAC) Topsoil physical properties for Europe (based on LUCAS topsoil data) (Ballabio et al 2016). This dataset includes the percentage of clay, silt and sand in the topsoil layer (upper 20 cm), and is based on interpolation of around 20.000 survey points all over Europe from the LUCAS topsoil database using Multivariate Additive Regression Splines.

Drought characterization 2.3.1. Standardized Precipitation Evapotranspiration Index (SPEI)
The Standardized Precipitation Evapotranspiration Index (SPEI) is a drought index that takes into account both precipitation and potential evapotranspiration (Vicente-Serrano et al 2010). The SPEI has been widely used because it allows the comparison among sites with contrasting climates and accounts for droughts at different time scales. The SPEI index is a standardized variate, i.e. a z-value, and expresses deviations from the current climatic balance (precipitation minus evapotranspiration potential) in respect to the long-term balance. Positive and negative SPEI values correspond to relatively wet and dry conditions, respectively. The SPEI values are computed globally on a grid for different time scales, depending on the time window that is used to calculate the SPEI values of the previous n months. Here, monthly SPEI values were obtained from the SPEIbase v.2.5, with 1-, 3-, 6-and 12-month time windows and a spatial resolution of 0.5°.

Soil volumetric water layer (swvl)
The soil volumetric water layer is the volume of water present in the total volume of soil (m 3 /m 3 ), divided into different soil layers. Compared to the SPEI, which is a statistical metric, the soil volumetric water layer is a soil physical variable which is related to plant available soil water. The volumetric soil water layer is a product of the European Centre for Medium-Range Weather Forecasts (ECMWF) ERA5-Land reanalysis (Copernicus Climate Change Service 2019) and is associated with soil texture, soil depth and the underlying groundwater level. The model used by ECMWF has a four-layer representation of soil: the first layer extending from 0 to 7 cm depth, the second layer from 7 to 28 cm, the third layer from 28 to 100 cm and finally the fourth layer from 100 to 289 cm depth. Here, the same partitioning is kept.

N deposition fields
The N deposition fields used in this study are modelled using the LOTOS-EUROS chemical transport model (Manders et al 2017). LOTOS-EUROS is an Eulerian chemistry transport model that simulates air pollution in the lower atmosphere. For this study, we used the mixed layer approach version of the model. We used a fivelayer grid that extends up to 5 km vertically. The bottom layer is the surface layer that has a fixed thickness of 25 meters. This layer is followed by a mixing layer, which is in turn followed by two dynamic reservoir layers of equal thickness and finally a top layer that extends up to 5 km. The model performs hourly calculations using ECMWF meteorology (ECMWF 2016), and uses the TNO CBM-IV gas-phase chemistry scheme (Schaap et al 2008). The anthropogenic emissions are taken from CAMS-REG-AP (Copernicus Atmospheric Monitoring Services Regional Air Pollutants) emissions dataset v2.2 (Granier et al 2019). The wet deposition computation simulates both in-cloud and below-cloud scavenging and is based on the CAMx approach (Banzhaf et al 2012). The dry deposition flux is computed using the DEPAC3.11 (Deposition of Acidifying Compounds) module (Van Zanten et al 2010). This module follows the resistance approach, in which the exchange velocity is equal to the reciprocal sum of the aerodynamic, quasi-laminar and canopy resistance. A canopy compensation point for NH 3 is implemented in this routine, allowing bi-directional fluxes. The model uses the CORINE/Smiatek land use map to compute exchange velocities for different land use classes. In this study, we computed 12 years (2003-2014) of land use specific fluxes for deciduous and coniferous forest in Europe (35°N to 70°N, 15°W to 35°E). We used a horizontal resolution of 0.50°longitude by 0.25°latitude, corresponding to approximately 28 by 28 km 2 . We match the modelled deposition fluxed with the FLUXNET sites based on their plant functional type. For mixed forest, the average flux for deciduous and coniferous forest is computed.

Deviations from reference GPP values
To quantify the response of forest GPP to drought, the observed daily GPP values are first compared to two types of reference GPP cycles, which we assume are 'expected' or default GPP values at a specific day on the long term. The deviations from these reference cycles are then computed and matched with drought indices to quantify the response of forest GPP to a specific drought. Here, two types of reference cycles are defined: (1) the detrended daily mean and maximum GPP values, and (2) the modelled daily GPP values, estimated from the observed meteorological variables at the same location using multiple regression analysis. The daily deviations are computed and then accumulated over a month or season (spring and summer), so that they can be linked to the monthly drought indices. The two types of reference GPP cycles are explained in more detail in the sections below.

Mean and maximum daily GPP values
The first method directly uses the daily FLUXNET measurement of the GPP. For each site, the mean and maximum GPP values are computed for each day of the year and used as GPP reference cycles. For the computation of the daily mean GPP value, all measurements with a medium-to high-quality indication are used (NEE_VUT_REF_QC>0.5). For determining the daily maxima, only observations with a high-quality flag indication (NEE_VUT_REF_QC>0.9) were used to avoid the inclusion of outliers. The monthly accumulated deviations from the daily mean GPP values will be referred to as δGPP 1 .

Modelled daily GPP values
In the second method, the reference GPP values are computed using multiple regression analysis. Here, daily FLUXNET observation of the ambient air temperature (T air ), incoming short-wave radiation SW in ( ) and longwave radiation (LW in ), the sensible-(H) and latent (LE) heat flux and the vapor pressure deficit VPD ( ) are used. These parameters are used as predictors to estimate the expected daily GPP at a certain day of the year given the meteorological conditions. As most biochemical reactions follow a hyperbolic curve, the following multi polynomial regression model is fitted to the data in a least-squares sense: The relationship between these predictors and the GPP changes throughout the year and may be quite asymmetrical, especially for a deciduous forest. To avoid systematic over-or underestimations in specific months, we fit the model for two separate periods: a startup growing period and the remainder of the growing season. In the startup period, we commonly observe a delayed response of GPP to change in these predictors. After this period, the relationship stabilizes and can be modelled well with a single multi-regression fit (see Result section). For simplicity, the split is made at the beginning of May for sites with deciduous forests and at the beginning of June for sites with coniferous forests.
The monthly accumulated deviations from the daily modelled GPP values will be referred to as δGPP 2 .

Detecting extreme drought events
Extreme droughts are detected using the following criterium: SPEI less than the 10th percentile of the sitespecific SPEI distribution. Here, only droughts that occur during spring and summer (March-August) are taken into account. The distribution of the different SPEI values (i.e., the monthly, 3-monthly, 6-monthly and 12monthly) is computed at each forest site during the time the measurements took place and the months in the lower 10th percentile are selected as drought events. The deviations from the reference cycles δGPP 1 and δGPP 2 during the months indicated as 'drought event' are then matched and (1) the distribution of δGPP 1 and δGPP 2 values and (2) the most extreme (negative) δGPP 1 and δGPP 2 values in the corresponding months are computed.

Comparison to N deposition, climate zone, soil type and forest age
After the deviations from the reference cycles δGPP 1 and δGPP 2 were computed and matched with the SPEI indices, they were compared to climatic and edaphic variables and N deposition estimates at each site. As a way of standardizing across forest sites, the δGPP 1 and δGPP 2 values are expressed as a percentage of the monthly mean GPP. As forest located in different climate zones respond differently to droughts of different durations, we determine which SPEI time window (1-, 3-, 6-or 12-months) shows the largest negative median δGPP 1 and δGPP 2 values at each site. The median and the largest negative, i.e., most severe, δGPP 1, 2 values corresponding to the SPEI with this time window are compared to the average amount of N deposition, the mean annual temperature (MAT), the mean annual precipitation (MAP), the age of the forest, and the percentage of clay and sand in the soil.

Detecting low soil moisture content
The soil volumetric water layer is used to look at the direct effect of soil water deficits on forest GPP during spring and summer. The analysis is as follows: first, the weighted averages of the soil volumetric water layer in the upper three layers (up to 1 m depth) at each site are computed. The weighted values are normalized so that 1 corresponds to the wettest-and 0 to the driest soil conditions occurring at that site, respectively. Subsequently, for each year the average weighted soil volumetric water layer during the spring (March-May) and summer (June-August) is computed, which we call the 'svwl WA '. Here, the daily maximum GPP values are used as reference values. For each year, the observed, total GPP during spring and summer is divided by the maximum GPP values during the same season, which we will refer to as 'fGPP max '. The svwl WA during spring and summer are then plotted against the corresponding fGPP max values. A simple linear regression model is fitted to all points and the slope and correlation coefficient are calculated. The computed slope and correlation coefficients are then compared to the average amount of N deposition, the mean annual temperature (MAT), the mean annual precipitation (MAP), the age of the forest, and the percentage of clay and sand in the soil for all sites simultaneously.

Relationship between N deposition and GPP and R eco fluxes
First, the yearly N deposition at each forest site is directly compared to the corresponding GPP and R eco values. The annual GPP values show a clear dependence on the annual N deposition values, increasing linearly at first, and peaking at approximately 10-15 kg N ha −1 yr −1 and 2000 gC m −2 yr −1 (figure 2). For higher N deposition values (>15 kg N ha −1 yr −1 ) a decrease in annual GPP is observed. This is in correspondence with the previously reported growth optimum related to N deposition. The relationship between R eco and N deposition shows more inter-site variability but follows a similar pattern. The optimum again lies around 10-15 kg N ha −1 yr −1 and around 1500 gC m −2 yr −1 . The modelled N deposition is then split up into NO y and NH x deposition. The breakdown of the N deposition into NO y versus NH x components at each FLUXNET site is shown in figure S3. The yearly NO y and NH x deposition are compared to the corresponding GPP values in figure 3. For the individual components, a similar relationship to GPP is observed. The optima for each component, however, lie at different values. The optimum GPP for NO y occurs around 8 kg N ha −1 yr −1 , while the optimum GPP for NH x lies around 12 kg N ha −1 yr −1 . Beyond these optima, the GPP decreases at different rates. Here, the decrease is particularly strong for NO y deposition values above 8 kg N ha −1 yr −1 , and relatively steeper than for NH x .

Changes in GPP during extreme droughts
The mean and maximum GPP cycles at each of the FLUXNET forest sites, that are used as reference values, are shown in figure 4. The GPP cycle of the forest sites with Picea abies as dominant species is the most symmetrical throughout the year and the least variable across sites. Forest sites with predominant deciduous broadleaved species, e.g., with Quercus and Fagus sylvatica as dominant species, show the most asymmetrical GPP cycles, as was expected. The GPP cycles of most forest sites with Quercus subspecies as dominant forest type are twofold, with a large initial peak during spring-and summertime followed by a smaller peak during autumn. Moreover, a clear dependence on climate is observed. Forest located in colder climates, for instance in Finland (Fi-Sod, Fi-Hyy), have a relatively short growing season. For forests in Mediterranean climates (e.g., IT-Cpz), on the other hand, an early onset of the growing season, as well as relatively high GPP values in wintertime ( figure 4(c)) are observed.
The distribution of the computed δGPP 1 and δGPP 2 values during months indicates as droughts are shown in figure S2. Table S1 shows an evaluation of the polynomial regression model per site, which includes the rootmean-squared deviation (RMSD), the mean absolute deviation (MAD), Pearson's correlation coefficient r and the slope between the measured and modelled daily GPP values. First of all, the spread of the δGPP 1 and δGPP 2 values plotted in figure S2 varies quite a lot across sites. There seems to be significant variability in the GPP responses to different drought durations as well. For example, there are more frequent, negative GPP responses corresponding to longer-lasting droughts in sites with Pinus spp. as dominant tree species (e.g. FI-Hyy, FI-Sod, FR-LBr), as well as sites dominated by Quercus spp. (FR-Pue, IT-Cpz, IT-Ro2). The difference between δGPP 1 and δGPP 2 values indicates to what extent the response of forest GPP to drought can be explained by the weather conditions. Some forest sites for instance strongly respond to high ambient air temperatures in combination with high vapor pressure deficit values with a consistent decrease in GPP, which could indicate direct stomatal closure at these conditions. This is visible for forest sites with Quercus species as dominant forest type in particular (i.e. FR-Fon, FR-Pue, IT-Cpz, IT-Ro1, IT-Ro2).
In the following section, only months identified as 'extreme droughts' that correspond to negative impacts on GPP are considered. Figures 5(a) and (b) show the most negative response in δGPP 1 and δGPP 2 at each forest site plotted against the mean N deposition, the MAT and the MAP, the forest age and the sand and clay content of the soil. Figure S2 shows the same information but grouped by dominant tree species. The largest negative δGPP 1 and δGPP 2 values differ the most for relatively young forest (< 80 years), young forest having a relatively larger negative response in δGPP 1 . This indicates that young forests tend to respond more rapidly and consistently to changes in weather conditions. At first sight, there is a decreasing response in forest GPP to drought with increasing N deposition values, which levels off around 15 kg N ha −1 yr −1 . However, the variation is large and the R 2 values of the best fitting curves are relatively small (R 2 =0.15 and R 2 =0.39). In general, no clear pattern for all sites in relation to either MAT or MAP, or sand or clay content is observed. Here, too, the relationships are rather weak.   Figures 6(a) and (b) show the mean negative response in δGPP 1 and δGPP 2 per drought event at each forest site plotted against the mean N deposition, the MAT and the MAP, the forest age and the sand and clay content of the soil. Figure S3 shows the same information but grouped by dominant tree species. In these figures, the  (δGPP1 and δGPP2), plotted against the mean N deposition (Ndep), mean annual temperature (MAT) and mean annual precipitation (MAP). The δGPP1 and δGPP2 values are plotted as percentage of the monthly mean GPP at each site, and represents the most negative δGPP1 or δGPP2 value out of the all values computed with 1-, 3-, 6-and 12-monthly SPEI values. The symbols represent the dominant forest types at each site. The gray lines represent the best fitting polynomial function using least-squares optimization. (b): Most severe, negative response in forest GPP to drought (δGPP1 and δGPP2), plotted against the forest age and the soil sand and clay content. The δGPP1 and δGPP2 values are plotted as percentage of the monthly mean GPP at each site, and represents the most negative δGPP1 or δGPP2 value out of the all values computed with 1-, 3-, 6-and 12-monthly SPEI values. The symbols represent the dominant forest types at each site. The gray lines represent the best fitting polynomial function using least-squares optimization.
values are plotted in g C m −2 per drought event, which corresponds to the mean decrease in GPP in one month identified as extreme drought (<10th percentile) by the used SPEI indices. First of all, all plotted relationships are rather weak. The mean negative response in forest GPP to drought increases with increasing MAT values, Figure 6. (a): The mean negative response in forest GPP to drought (δGPP1 and δGPP2) per drought event, plotted versus the amount of N deposition (Ndep), the mean annual temperature (MAT) and the mean annual precipitation (MAP). The GPP response per drought event is computed by summing all negative δGPP1 and δGPP2 values and then dividing that by the number of drought events per SPEI index (1-, 3-, 6-and 12-monthly). The gray lines indicate the spread in outcomes using different SPEI indices. The symbols represent the dominant forest types at each site. (b): The mean negative response in forest GPP to drought (δGPP1 and δGPP2), plotted versus the forest age and the soil clay and sand content. The GPP response per drought event is computed by summing all negative δGPP1 and δGPP2 values and then dividing that by the number of drought events per SPEI index (1-, 3-, 6-and 12-monthly). The gray lines indicate the spread in outcomes using different SPEI indices. The symbols represent the dominant forest types at each site. which indicates that forest sites located in warmer climates experience relatively strong GPP reductions during extreme droughts. Despite significant uncertainties, there appear to be some detectable relationships in forest sites with the same dominant tree taxa. For example, at forest sites with Pinus spp. as dominant species there seems to be a linear relationship between the mean negative δGPP 1 and the mean annual temperature and precipitation. Furthermore, in oak forests (Quercus ssp.), and to a lesser extent in spruce forests (Picea ssp.), linear relationships between the mean negative δGPP 1 and the soil sand and clay content can be observed at each site. In these forests, the response of forest GPP to drought seems to become more severe with decreasing soil sand content and increasing clay content. This suggests that Quercus ssp. are in general more sensitive to changes in soil water content and that Quercus ssp. forests on sandy soils are generally less susceptible to droughts than those on clayey soils. Figures 7(a) and (b) show the slope and the correlation between the svwl WA and fGPP max during spring and summer at each forest site. Positive slopes and correlation coefficients represent a positive relationship between the svwl WA and fGPP max , e.g., increasing GPP with increasing soil volumetric water layer. In general, GPP in most forest sites benefits from relatively low soil moisture levels during spring, and relatively high soil moisture levels during summer.

Effect of soil volumetric water layer on seasonal GPP
The highest correlations between svwl WA and fGPP max in spring are found for sites with Fagus sylvatica as dominant species (diamond). The growth during the beginning of the growing season is thus highly dependent on the amount of available soil water at sites with Fagus sylvatica as the dominant tree species. Generally, the soil volumetric water content in the topmost layer (>1 meter below the surface) decreases gradually from wintertime towards the summer as accumulated evaporation becomes higher than accumulated precipitation. As a result, trees are thus less likely to experience serious water shortages during spring. Moreover, as the start of the growing season seems to be fairly consistent at these sites (see figure 4), the apparent positive effect of low soil moisture conditions on GPP in spring is likely caused by favorable meteorological conditions, such as higher air temperatures and higher amounts of incoming short-wave radiation, and not low soil moisture conditions. During summertime, the relationship between svwl WA and fGPP max is mainly dependent on the mean annual temperature (MAT) and the forest age. The slope and correlation coefficient increase with increasing MAT values, indicating that forest sites located in warm climates are relatively sensitive to changes in soil moisture. The slope and the correlation coefficients decrease with increasing forest age, suggesting that the GPP in younger forests is more sensitive to changes in soil moisture compared to older forests. This is in line with the literature, as young forests have less developed root systems. In general, there are no clear pattern in relation to N deposition.

Case study: N deposition and drought interactions in pine forest (Pinus ssp.)
The results from the previous sections suggest that there is a differential response of forest GPP to drought across FLUXNET sites, seemingly independent from N deposition. However, by grouping sites with the same dominant tree taxa, some patterns could be distinguished. In this section, the group of sites located in pine forests (Pinus spp.) is studied in more detail. These pine forest sites (FI-Hyy, FI-Sod, IT-SRo, FR-LBr, NL-Loo, BE-Bra) showed the strongest relationship between the response in forest GPP to drought and the MAT, MAP and forest age. Moreover, the largest variation in N deposition levels is observed across these sites. Figure S5 shows the time series of the daily and monthly summed GPP deviations (δGPP 1 ). At some sites, there are consistently higher or lower GPP values throughout the entire year (IT-SRo, FR-LBr, FI-Sod). Other sites, on the other hand, show much more intra-annual variability (NL-Loo, FI-Hyy). Figure 8 shows the mean negative response in forest GPP (δGPP 1 ) per drought event identified by different drought indices. Here, the normalized soil volumetric water layer (svwl) values at each forest site is used as an additional drought indicator. Most sites respond more severely to longer-lasting droughts, except for BE-Bra. The addition of the svwl did not lead to a more severe response in forest GPP to drought at all sites. BE-Bra, NL-Loo and FI-Sod are insensitive to svwl as an additional drought indicator. Figure 9 shows the mean negative response in forest GPP (δGPP 1 ) per drought event, in relation to the MAT, MAP, forest age and the amount of N deposition. Overall, the response in GPP increases with increasing MAT and MAP values and decreasing with increasing forest age. To extract the influence of the MAT, MAP and forest age, a multi-linear regression was fitted. The resulting, modelled δGPP 1 are plotted in white in relation to N deposition, alongside the initial δGPP 1 values. From our hypothesis, we expect the response in forest GPP to drought to be more severe at forest sites with very low levels of N deposition (FI-Hyy and FI-Sod), and with very high levels of N deposition (BE-Bra and NL-Loo). At forest sites with intermediate levels of N deposition (IT-SRo and FR-LBr), N deposition is expected to alleviate the drought response in GPP. The fitted δGPP 1 values, however, at different levels of N deposition, are both higher and lower than the observed δGPP 1 values.
Moreover, the response of forest GPP to drought was modelled with another multi-linear regression, this time including the amount of N deposition. The resulting fitted model (δGPP 1 =5.26 * MAT − 0.38 * MAP − 0.41 * forest age +1.13 * N dep +225.36) improved slightly, but not significantly (R 2 =0.85, RMSE=8.1) and shows that N deposition only explains another ∼10% of the variation in the data. These results suggest that there is no consistent nor significant effect of N deposition on the response of forest GPP to drought in these pine forest sites and therefore do not support our hypothesis. The slope and the correlation coefficient between the weighted average normalized soil volumetric water layer (swvlWA) and the fraction of the maximum GPP (fGPPmax) during spring and summer at all forest sites, plotted against the amount of N deposition (Ndep), the mean annual temperature (MAT) and the mean annual precipitation (MAP). Negative slopes and r values (red) indicate a preference for drier soil conditions, whereas positive slopes and r values (blue) indicate a preference for wetter soil conditions. The symbols represent the dominant forest types at each site. The slope and correlation coefficients are obtained by fitting a simple linear regression using least-squares optimization. (b): The slope and the correlation coefficient between the weighted average normalized soil volumetric water layer (swvlWA) and the fraction of the maximum GPP (fGPPmax) during spring and summer at all forest sites, plotted against the forest age and the soil sand and clay content. Negative slopes and r values (red) indicate a preference for drier soil conditions, whereas positive slopes and r values (blue) indicate a preference for wetter soil conditions. The symbols represent the dominant forest types at each site. The slope and correlation coefficients are obtained by fitting a simple linear regression using least-squares optimization.

Discussion
This study discusses various pathways in which N deposition and drought can interact as co-stressors on forest productivity in European forest ecosystems. Based on the literature, our initial hypothesis was that the response of forest GPP to drought is relatively small in forest ecosystems that experience intermediate levels of N deposition (5-10 kg ha −1 yr −1 ). Productivity in forest ecosystems with either very low N deposition levels (<5 kg ha −1 yr −1 ) or higher N deposition levels approaching N saturation (>10 kg ha −1 yr −1 ) was expected to show a relatively large magnitude response to drought. Due to a differential response to the used drought indices, we were unable to isolate a clear, overall dependence of the response of forest GPP to drought on N deposition. This, in part, results from the large inter-site variability in dependencies of the response of forest GPP to drought on varying physiological, climatic and edaphic conditions. Grouping sites with the same dominant tree species enabled us to extract some of these dependencies, such as the sensitivity of drought response to soil texture in oak (Quercus ssp.) forests. Furthermore, the GPP response to drought within the group with the largest variation in N deposition (Pinus ssp.) could be predicted using the site mean annual temperature and precipitation and the age of the forest (R 2 =0.76, RMSE=10.0). However, after extracting the contribution of these most prominent parameters, the variability in drought responses remained too large to substantiate any hypothesized N deposition effects.
We were unable to draw any generalized conclusions regarding the impact of N deposition on the response of GPP to drought in European forest ecosystems. In addition to the large inter-site variability in physiological, climatic and edaphic conditions in general, another limiting factor in this study setup was the low variation in N deposition levels within some groups of dominant forest species (e.g. Quercus ssp. or Fagus sylvatica). Moreover, even though the FLUXNET2015 dataset contains the most extensive network of eddy-covariance measurements to date, this study is limited by the relatively short time series at some forest sites. As a result, some forest sites for instance only experienced a handful of months indicated as severe droughts. As we established that for example both frequency and severity of a particular drought event play an important role in the forest response, it is especially hard to generalize the response with only a few droughts occurring. Moreover, a relatively short time series of GPP observations could result in a bias in the mean GPP cycle and therefore the drought-induced anomalies in GPP. For example, if the majority of the years with measurements at a particular site were relatively warm and dry, the difference in GPP with a classified drought period would be smaller than if these years were relatively cold and wet. Also, longer time series would allow us to look at the impact of re-occurring or longerlasting droughts, and the impact of severe drought on the long-term C cycling in forest ecosystems.
A recent review by Speich et al 2019 found that, for temperate forests, drought indices accounting for evaporative demand performed better than indices based on precipitation alone. The SPEI index used here does account for both precipitation and evaporative demand, but does not include any additional levels of information. The results showed that not all month indicated as extreme droughts by the SPEI-indices correspond to negative responses in GPP ( Figure S2). This is likely because the SPEI index does not account for soil moisture storage, the ease of soil water extraction by plant roots, or water table depth. Therefore, the first couple of months of drought, with possibly sufficient water availability coming from deeper soil layers, were included in the analysis. As long as tree root systems have access to sufficient amounts of soil water, favorable meteorological conditions during drought, such as higher incoming solar radiation, could even result in a relative increase of GPP. This illustrates that using the 10th percentile of the SPEI index as a definition of extreme drought conditions, as was done in this study, does not necessarily only include impactful droughts. Using a more applied drought index that uses more site-specific variables (e.g., soil moisture storage, stand properties, physiological thresholds) could potentially improve the selection of droughts.
Other factors contributing to the uncertainty in our analyses are related to the modelling of N deposition. The modelled N deposition fields used in this study have a relatively large uncertainty. N deposition Figure 8. The mean negative response in forest GPP (δGPP 1 ) per drought event identified by different drought indices. A specific month is classified as drought if the used drought index is below its 10th percentile value. The soil volumetric water layer (svwl) values are standardized. measurements, both the dry and wet component of the deposition fluxes, at the same locations as the FLUXNET sites could help to improve on this and would allow us to represent the current N deposition levels per site more accurately. At the same time, most of these forests are around a 100 years old, and better estimates or measurements of current N deposition levels are, therefore, not necessarily correlated with historic N deposition levels. Here, we used a relatively long time-series (12 years) of modelled N deposition and assumed that these N deposition levels represent the N availability in these forests as a whole relatively well. N retention and allocation in forest ecosystems is, however, not only influenced by the current level of N deposition, but also by the accumulated N deposition that a forest ecosystem has received (Vries et al 2014). Using additional measurements of N in plant and soil pools could help us to get a better grip of the accumulated N in different forest compartments and herewith the historic N deposition levels.
Nutrient availability in general is a key regulator of the forest carbon balance. Not only N limitation plays a role in forest productivity, but also the availability of other nutrients (e.g., potassium, phosphorus, calcium, magnesium) (Fernandez-Martinez et al 2014). For a more complete analysis, not only N availability but also other potentially limiting nutrients such as phosphorus (P) or other cations (Mg, Ca and K) could be considered. Moreover, surface ozone (O 3 ) also plays a critical role in forest growth and drought response (e.g. (Karlsson et al 2004, Kronfuß et al 1998). In the current study setup, however, adding more variables (for instance O 3 surface concentrations) would likely result in an even greater divergence in conditions between the FLUXNET forest sites. To effectively isolate N deposition-drought interactions, an extensive measurement network in forests with similar edaphic, physiological and climatic conditions is needed. To get more insight into N deposition-drought interactions for specific tree species, one could for instance set up N addition experiments in a managed forest. Another option would be to do measurements for the same type of trees in regions with strong, local N deposition gradients. This is the first study to examine both the effects of modelled N deposition on forest productivity across the European FLUXNET sites and the possible effects of N deposition on the magnitude of a productivity response to drought. We find good agreement across the FLUXNET sites in the dependence of GPP to N deposition and its component (NO y and NH x ), with an initial strong increase in GPP at low N deposition levels followed by a slow decline, and even a decrease at high N deposition levels. However, the effect of N deposition on the magnitude of the forest productivity response to drought could not be isolated. These results suggest that, while N deposition might play a critical role in the response of forest productivity to drought within specific forest ecosystems, N deposition does not seem to be a major or consistent driver of the magnitude of the GPP response to drought across a diverse set of FLUXNET forest sites in Europe. Figure 9. The mean negative response in GPP (δGPP 1 ) per drought month in relation to the mean annual temperature (MAT), the mean annual precipitation (MAP), the forest age and the N deposition. Each triangle stands for one forest site, and the grey lines indicate the spread in outcomes using different SPEI values. The triangles with the white outline are the modelled δGPP 1 values, using a multi-linear regression fitted to the MAT, MAP and forest age (δGPP 1 =6.29 * MAT − 0.35 * MAP − 0.02 * forest age+183.43, R 2 =0.76, rmse=10.0).