How Nitrogen and Phosphorus Availability Change Water Use Efficiency in a Mediterranean Savanna Ecosystem

Nutrient availability, especially of nitrogen (N) and phosphorus (P), is of major importance for every organism and at a larger scale for ecosystem functioning and productivity. Changes in nutrient availability and potential stoichiometric imbalance due to anthropogenic nitrogen deposition might lead to nutrient deficiency or alter ecosystem functioning in various ways. In this study, we present 6 years (2014–2020) of flux‐, plant‐, and remote sensing data from a large‐scale nutrient manipulation experiment conducted in a Mediterranean savanna‐type ecosystem with an emphasis on the effects of N and P treatments on ecosystem‐scale water‐use efficiency (WUE) and related mechanisms. Two plots were fertilized with N (NT, 16.9 Ha) and N + P (NPT, 21.5 Ha), and a third unfertilized plot served as a control (CT). Fertilization had a strong impact on leaf nutrient stoichiometry only within the herbaceous layer with increased leaf N in both fertilized treatments and increased leaf P in NPT. Following fertilization, WUE in NT and NPT increased during the peak of growing season. While gross primary productivity similarly increased in NT and NPT, transpiration and surface conductance increased more in NT than in NPT. The results show that the NPT plot with higher nutrient availability, but more balanced N:P leaf stoichiometry had the highest WUE. On average, higher N availability resulted in a 40% increased leaf area index (LAI) in both fertilized treatments in the spring. Increased LAI reduced aerodynamic conductance and thus evaporation at both fertilized plots in the spring. Despite reduced evaporation, annual evapotranspiration increased by 10% (48.6 ± 28.3 kg H2O m−2), in the NT plot, while NPT remained similar to CT (−1%, −6.7 ± 12.2 kgH2O m−2). Potential causes for increased transpiration at NT could be increased root biomass and thus higher water uptake or rhizosphere priming to increase P‐mobilization through microbes. The annual net ecosystem exchange shifted from a carbon source in CT (75.0 ± 20.6 gC m−2) to carbon‐neutral in both fertilized treatments [−7.0 ± 18.5 gC m−2 (NT) 0.4 ± 22.6 gC m−2 (NPT)]. Our results show, that the N:P stoichiometric imbalance, resulting from N addition (without P), increases the WUE less than the addition of N + P, due to the strong increase in transpiration at NT, which indicates the importance of a balanced N and P content for WUE.

Current and projected levels of anthropogenic N deposition are expected to cause stoichiometric imbalances in plant-available N and P in several terrestrial ecosystems (Du et al., 2020;Peñuelas et al., 2010Peñuelas et al., , 2013. In such a scenario, plants must adapt to altered and potentially imbalanced nutrient availability (Elser et al., 2010;Güsewell, 2004;Oldroyd & Leyser, 2020;Zhu et al., 2016). In terms of carbon and water fluxes, ecosystem-level responses due to nutrient availability and stoichiometric N:P imbalances are poorly understood and may differ from those at the leaf, plant, or community level, due to interactions and compensatory effects among and between plants and soil. At the ecosystem level, it is especially important to understand and characterize how resource-use efficiencies such as WUE are changing with the stoichiometric imbalance and climate change. Especially in semi-arid ecosystems, increases in WUE could have positive effects on carbon sequestration potential as more carbon can be fixed with the same limited amount of available water (Grunzweig et al., 2003). In addition, increased WUE allows for maintaining photosynthesis while losing less water through transpiration. This leads to less soil moisture depletion and avoids the earlier onset of water stress and senescence in summer (Luo et al., 2020). There is a physiological trade-off between nutrient-use efficiency and WUE (Field et al., 1983;Han et al., 2016) arising from the importance of stomatal openings in enhancing the former while decreasing the latter. However, it is unclear to what degree an imbalance of the N:P ratio at ecosystem level will influence transpiration and WUE or how it is regulated by climate variability (Fernández-Martínez et al., 2014;Jiang et al., 2019;Luyssaert et al., 2014;Peñuelas et al., 2010Peñuelas et al., , 2013. This knowledge gap owes to the small spatial scales of experiments that manipulate nutrient availability, which have typically ranged from individual leaves and plants in mesocosm experiments to the order of a few tens of meters in small plot studies (Fay et al., 2015;Ford et al., 2016;Wicklein et al., 2012). While these experiments provide answers for nutrient manipulations for specific ecosystem compartments, such as soil and plants, they do not allow for the quantification of ecosystem-scale responses by carbon and water fluxes or water-use efficiency (WUE) in an integrated system at high temporal resolutions. The reliability of previous ecosystem-scale studies of surface properties and N effects on ecosystem-atmosphere exchange has been widely debated. Many studies were based on natural spatial variability in N content and may be influenced by confounding factors, e.g., temperature, species, and stand density (de Vries et al., 2008;Luyssaert et al., 2014;Magnani et al., 2007Magnani et al., , 2008. While studies concerning the relationship of leaf N and ecosystem functioning exist and can be improved upon, systematic studies of the role of P and N:P ratio imbalances at the ecosystem level are lacking (Du et al., 2020;Penuelas et al., 2020). In general, investigations of fluxes or resource-use efficiency at the ecosystem level, under controlled N and P nutrient treatments are missing.
This study aims to address these knowledge gaps and comprehensively quantify the effects of concomitant increasing N availability and increasing N:P imbalance on ecosystem functioning and on WUE in a semi-arid tree-grass ecosystem.
EL-MADANY ET AL. 10.1029/2020JG006005 2 of 21 N + P (21.5 ha). A third plot served as the control treatment. While the N-only treatment created an imbalance between the available N and P, this imbalance was relieved in the N + P treatment where both N and P were provided. Our measurements showed that both fertilized treatments increased their carbon uptake and turned the ecosystem from a carbon source to carbon neutral. One of the main differences between the fertilized treatments which is associated with the imbalance of available N and P is the loss of water through the vegetation (transpiration). This increase in transpiration was only observed in the N-only but not in the N + P treatment. Our results show, that the N:P stoichiometric imbalance, resulting from N-only addition, increases the water-use efficiency (i.e., the carbon gain per water loss) less than the addition of N + P, due to the strong increase in transpiration at the N-only treatment.
We present here, 6 years (2014-2020) of data from a large-scale N and P manipulation experiment, in which we combine eddy covariance data with satellite and plant trait data. We measured and compiled ecosystem responses of water fluxes [evapotranspiration (ET), transpiration (T), evaporation (E)], carbon fluxes [net ecosystem exchange (NEE), gross primary productivity (GPP), ecosystem respiration (Reco)], and different metrics of WUE (i.e., the carbon uptake per unit of water lost), with changes in canopy structure (leaf area index, LAI) and foliar nutrient content and isotopes, mainly during the winter and spring period.
We analyze the available data streams to evaluate if and how changing nutrient availability and stoichiometric N:P imbalance influences carbon and water fluxes, canopy structure and WUE at ecosystem scale.

Site Description
The study site, Majadas de Tiétar (Casals et al., 2009;El-Madany et al., 2018;Perez-Priego et al., 2017), is located in the center of the Iberian Peninsula (39°56′25″N 5°46′29″W). It is a tree-grass ecosystem, a typical "Iberic dehesa," with 20-25 Quercus ilex trees ha −1 . The average canopy height is 8.7 ± 1.25 m, and the fractional canopy cover is 23.0 ± 5.3% (Bogdanovich et al., 2021). The dehesa is managed and used for extensive cattle farming with a cow density of ≤ 0.3 ha −1 . Precipitation and temperature were measured at the site between 2004 and 2019. The mean annual precipitation was 636 mm, ranging from 440 to 965 mm. Nearly 85% of the annual precipitation fell between October and April. The mean annual temperature was 16.7°C with a mean annual minimum temperature of −4.7°C and a mean annual maximum temperature of 41.1°C. The mean LAI of the green vegetation (LAI green ) changes throughout the seasons, with a mean value between 0.25 ± 0.07 m 2 m −2 in summer and 1.75 ± 0.25 m 2 m −2 at the peak of the growing season in spring, for the herbaceous layer and between 1.5 and 2.0 m 2 m −2 for the trees. Due to their fractional tree canopy cover, this converts to 0.3-0.4 m 2 m −2 for the total area. The trees have a rather constant LAI throughout the year and seasonal variability of the ecosystem LAI is driven by the phenology of the herbaceous layer (Luo et al., 2018). The growing season begins with the re-greening of the herbaceous layer after the autumn rains start, which usually is in October (Figure 1). The growing season lasts usually until the end of April to the end of May when the depletion of soil moisture, high radiation and temperatures lead to the senescence of the herbaceous layer.
The soils are characterized as Abruptic Luvisol with a sandy upper layer . The prevailing wind directions are west-southwest and east-northeast.

Experimental Design and Fertilization
The study-site was divided into three plots. N was added to one plot (NT), N + P to another one (NPT), and the third plot served as the unfertilized control area (CT).
The fertilizers were applied with a tractor around the eddy covariance towers. The fertilized areas were 21.5 and 16.9 ha at NPT and NT, respectively ( Figure 2   Small black dots are half hourly data and the big colored points correspond to weekly averages. Ta is the air temperature in 2 m (a), VPD is vapor pressure deficit in 2 m (b), SWCn is the normalized soil moisture content for the top 20 cm (c), Rain are daily precipitation sums (d), H is the sensible heat flux (e), LE is the latent heat flux (f), NEE is the net ecosystem exchange of carbon dioxide (g) big green dots are daytime averages and big brown dots are night time averages. The vertical dashed lines correspond to the winter (blue) and spring period (red). NT received 100 kg N ha −1 in the form of ammonium nitrate (NH 4 NO 3 ) and calcium ammonium nitrate (5Ca(NO 3 ) 2 NH 4 NO 3 ), respectively. Calcium ammonium nitrate was used at NT to account for the calcium included in the triple superphosphate fertilizer used at NPT. A second dose of fertilizer was applied in March 2016, 10 kg P ha −1 and 20 kg N ha −1 , using the same fertilizers as before. No further additions of nutrients were made to follow the outcome of the added nutrients and the ecosystem response in time.

Eddy Covariance Setup
Three ecosystem eddy covariance (EC) towers are operated at the site, one in each plot. The long-term FLUXNET site Majadas de Tiétar (ES-LMa in FLUXNET, since 2003) serves as CT, and two additional towers are added for the fertilization experiment and are operated since March 2014. These EC towers are named ES-LM1 (NT) and ES-LM2 (NPT). Each tower is equipped with a R3-50 sonic anemometer (Gill Instruments Limited, Lymington, UK) to measure three-dimensional wind components and sonic temperature, and a LI-7200 infra-red gas analyzer (Licor Bioscience, Lincoln, Nebraska, USA) to measure CO 2 and H 2 O mixing ratios. The measurement height is at all towers 15 m above ground.

Data Treatment, Gap-Filling and Partitioning
Eddy covariance data and meteorological data were collected and calculated, as described in El-Madany et al. (2018). The processing of the raw high-frequency data was conducted using EddyPro v 6.2.0 (Fratini & Mauder, 2014). Calculated CO 2 fluxes were storage corrected from seven-point profile measurements.
EL-MADANY ET AL.  A friction velocity (u*) threshold [0.194 ± 0.03 ms −1 (CT), 0.200 ± 0.03 ms −1 (NT), 0.207 ± 0.02 ms −1 (NPT); mean ± sd] was determined for each year following Papale et al. (2006), and all data with u* values below the established threshold were removed from further analysis. The gap-filling of missing and bad quality data [QC values > 1 according to Mauder and Foken (2011)] was then performed (Reichstein et al., 2005). Subsequently, flux partitioning of NEE into GPP and Reco was conducted by applying the nighttime flux partitioning method (Reichstein et al., 2005). All post-processing was conducted using the R package REd-dyProc v 1.2.2 (Wutzler et al., 2018). The details of the full processing including statistics on the number of gaps for the individual fluxes are provided in the Supplementary Material and Figure S1.

Flux Footprint Calculations
To estimate which areas were sampled by the EC systems and to evaluate how similar these areas were the flux footprint (FP) model of Kljun et al. (2015) was used, as integrated into the Flux Footprint Prediction (FFP)-R-code version 1.4 (available at http://footprint.kljun.net/download.php). We calculated the FPs for all half-hourly (HH) periods between March 2014 and February 2020, when all input parameters (wind speed, wind direction, the standard deviation of lateral wind component, u*, Obukhov length and boundary layer height) for the model were available and friction velocity values were above the u*-threshold. The boundary layer height estimates were used from ERA5 at 0.25-degree spatial resolution (Hersbach et al., 2018) and linearly interpolated to fit the 30-min temporal resolution of the EC data. The flux footprint climatology (FFC) was subsequently calculated based on 64271 (CT), 66920 (NT), and 67862 (NPT) individual HH footprints, by integrating the flux contribution from the peak location into all directions until that is, 80% of the total contribution was reached. The 80% contribution was selected because it was the largest at which no overlap between the treatments was observed ( Figure 2). The size of the FFC for the three plots were 26.0 ha (CT), 23.2 ha (NT), and 24.3 ha (NPT). At NT and NPT the fertilized areas covered 16 and 19.1 ha centered on the respective FFC.
The FFC was intersected with Landsat 7 Normalized Difference Vegetation Index (NDVI) maps (see 2.4.1) and classification maps from airborne hyperspectral measurements (see 2.4.2 and Pacheco-Labrador et al. (2017;2020)) to calculate, for example, average NDVI values, tree canopy cover and FFC size. The medians of the distribution of tree and grass fraction from the individual HH FPs of all treatments agree within 2% with the fractional cover within the FFC. Between treatment differences of the fractional covers were below 6% ( Figure S2) for the FFC as well as the medians of the HH flux footprints.

Landsat Data
The Landsat 7 NDVI at 30 m spatial resolution (Masek et al., 2006) data were obtained from Google Earth Engine Data Catalog (https://developers.google.com/earth-engine/datasets/catalog/LANDSAT_LE07_ C01_T2_SR#description) and were used to characterize pre-fertilization (1999-2014) and post-fertilization (2015-2019) differences within the 80% isolines of the FFC at each tower ( Figure S3). First, a subset including the outer margins of the 80% FFCs was taken from each scene. Depending on their size individual FFC included 264 to 294 30 × 30 m Landsat 7 pixels. After May 31, 2003, the Scan Line Corrector (SLC) of Landsat 7 broke, which resulted in a reduction of available pixels per scene. The data gaps form alternating wedges that increase in width from the center to the edge of a scene and change position in the different overpasses. Additionally, the following quality control scheme was used to ensure that pixels within the subsetted scenes were cloud-free and of good quality: Atmospheric opacity ≤0.1, cloud cover ≤2, pixel quality ≤66 and pixel saturation = 0. In case that more than 20% of the pixels were flagged as bad quality the whole scene was removed from the analysis. The lowest number of good quality pixels within the FFCs for an individual scene ranged between 60% and 68% of all available pixels (including data gaps due to SLC error). Average NDVI values were calculated for each treatment's FFC when the quality of the respective scene was good. Overall, 170 and 92 scenes fulfilled the quality criteria for the pre-and post-fertilization periods, respectively.

Airborne Spectral Measurements
Airborne hyperspectral measurements were conducted using a Compact Airborne Spectrographic Imager CASI-1500i (CASI) (Itres Research Ltd., Calgary, AB, Canada) with a spatial resolution of 1.25 m and 140 spectral bands ranging from 380 to 1,050 nm (Pacheco-Labrador et al., 2017, 2020. The CASI was on board a C-212-200 RS airplane, which was operated by the Spanish Institute for Aerospace Technology. Hyperspectral flight campaigns were conducted on May 5, 2011 (before the experiment), April 8, 2014 (during the experiment but before the fertilization), and April 23, 2015 (after the fertilization) around solar noon and using similar flight patterns. The Envisat MEdium Resolution Imaging Spectrometer (MERIS) Terrestrial Chlorophyll Index (MTCI) was calculated based on the CASI data, to have an estimate of canopy chlorophyll content in the different treatments (Dash & Curran, 2007).
Here, ρ is the hemispherical-directional reflectance factor for a specific wavelength in nm. The numbers denote the wavelengths. Similarly, NDVI was calculated from the CASI data following Rouse et al. (1974).
For details about the processing of the hyperspectral images we refer to Pacheco-Labrador et al. (2017, 2020. The number of CASI pixels within each FFC ranged, depending on its size, between 148763 -166350. MTCI and NDVI values from pixels classified as water bodies were excluded from the analysis. Differences between distributions were then used to estimate the effect of the fertilization ( Figure S4).

Metrics of Photosynthetic Capacity
The maximum net CO 2 uptake rate of the canopy at infinite global radiation (A max ) was calculated as part of the daytime flux partitioning (Lasslop et al., 2010) where Rg is the global radiation, α the initial slope of the LRC and γ the Lloyd and Taylor respiration model (Lloyd & Taylor, 1994). We want to emphasize, that A max includes a VPD dependency based on an exponentially decreasing function at VPD values > 10 hPa).
To calculate the GPP value at light saturation that is, R g = 1,000 W m −2 (GPP sat ), the "light.response" function of the "bigleaf" R-package (Knauer et al., 2018 Equation 29) was used as shown in Equation 4.
where α is the initial slope of the LRC, R gRef the reference R g of 1,000 W m −2 , and R eco the ecosystem respiration. The parameters α and GPP sat are derived using the "nls" function in R (R Development Core Team, 2015). GPP sat was calculated for every day using a three-day moving window.

Water-Use Efficiency Metrics
We used the stomatal slope parameter G 1 as a physiologically significant WUE metric from the eddy covariance data . It is calculated based on the stomatal model of Medlyn et al. (2017) and is inverse to the WUE.
s G is the surface conductance for water vapor, 0 G is the minimum canopy conductance, 1 G is the stomatal slope parameter, VPD the vapor pressure deficit, Ca the CO 2 concentration and GPP is the gross primary productivity.
For the calculation of G 1 , the R package "bigleaf" was used . To ensure meaningful estimates of G 1 , data were first filtered for the following reasons: 1) to only include measured, flux and meteorological input parameters (no gap-filled data); 2) to remove 48 h after precipitation; and to ensure that values were within the following ranges 3) global radiation (Rg) > 100 W m −2 ; 4) GPP > 1 µmol m −2 s −1 ; 5) VPD > 0.1 kPa; 6) LE and net radiation (Rn) > 0 W m −2 ; and 7) relative humidity (rH) < 95%. After filtering, 18% (CT), 20% (NT), and 20% (NPT) of the site data remained. The difference of 2% in data availability was due to a longer data gap at the CT plot, which was caused by lightning in the summer of 2014. G s was calculated based on the inverted Penman-Monteith equation (Knauer et al., 2018 Equation 15), where aerodynamic conductance was estimated, according to Thom et al. (1972). G 1 was then calculated for the spring of each year, individually, to compare the changes in G 1 between sites and pre-and post-fertilization. The uncertainty of G 1 was estimated by randomly resampling 10% of the data used to calculate G 1 . This procedure was repeated 100 times for each spring season of each year (2014-2019).
To estimate the fertilization effect on G 1 the pre-fertilization G 1 estimate for spring 2014 was subtracted from each post-fertilization G 1 estimate for the respective spring periods (2015-2019). Thus, the differences between G 1 at the CT plot and fertilized treatment plots (NT, NPT) correspond to the fertilization effect. To compare whether different WUE metrics detected the same patterns we also calculated inherent WUE (iWUE, Beer et al., 2009) andunderlying WUE (uWUE, Zhou et al., 2014) for the same time periods as G 1 . The descriptions of the different WUE formulations are described in the SI.

Transpiration Estimates
The transpiration estimation algorithm (TEA) was used to partition measured ET into T and E (Nelson et al., 2018). TEA operates by modeling ecosystem WUE T (as GPP/T) from GPP and ET data after controlling for E contamination in the ET signal. First, the data is filtered for periods when plants are active (i.e., GPP > 0.5 gC m −2 d −1 ) and when ecosystem surfaces are probably dry to minimize the influence of abiotic evaporation. A machine learning model [random forest (Breiman, 2001)] was then trained on the WUE ET (as GPP/ET) within this filtered subset using metrics derived from the eddy covariance fluxes and associated meteorological data as predictors. The trained model then predicts WUE T at every time step using a 75% quantile prediction (Meinshausen, 2006) to further account for the effect of residual E in the training data set. The 75% quantile was chosen based on previous experiments using the model output as synthetic eddy covariance data when T was known (Nelson et al., 2018). T is then ultimately calculated as GPP/WUE T . See Nelson (2020) for the related code and tutorials describing the use of TEA.

Vegetation Data
Every year, the herbaceous layer was sampled in spring to measure biomass, green LAI and to analyze nutrient content. At each sampling date (Table 1), the aboveground biomass was harvested in multiple 25 × 25 cm squares of the herbaceous layer within the 25 × 25 m plots in each FFC (Figure 2). A subsample of the total biomass was selected and green and senescent fractions were separated by hand in the laboratory. The green and senesced fractions were weighted, scanned, dried and weighted again. The LAI green was then calculated as With W DT the dry weight of the total sample, W DSG the the dry weight of the green fraction of the subsample, W DS the dry weight of the subsample, LAI greenSub the leaf area of the green fraction of the subsample and A T the area of the sample plot. Samples were oven-dried at 60°C for 48 h. Plant material was ground using a ball mill (Co. Resch) at a frequency of 20 Hz for less than 4 min. Subsequently, the material was stored for further chemical analysis.
N concentration (mg/g) was analyzed by the dry combustion method (Vario EL; Elemantar GmbH, Hanau, Germany) using 25 mg of the ground sample material. P concentrations were analyzed using an inductive coupled plasma optical emission spectrometer (ICP OES, Optima 3300, Perkin Elmer). First 0.1 g of the dried and ground plant material was weighed and 3 ml of HNO 3 were added. The sample was then heated for 20 min with 1,000 W and then cooled down for 10 min. Afterward, the sample was filtered and diluted with deionized water to get a 50 ml sample. The prepared sample was then analyzed in the ICP OES.
Isotopic signatures are based on the determination of the carbon isotopic composition of dried grass and tree foliar samples. They were analyzed with a DeltaPlus isotope ratio mass spectrometer (Thermo Fisher, Bremen, Germany) coupled via a ConFlowIII open-split to an elemental analyzer (Carlo Erba 1100 CE analyzer; Thermo Fisher Scientific, Rodano, Italy). The measurement protocol is based on Werner et al. (1999;2001) and Brooks et al. (2003). Carbon isotope signatures of tree and grass samples (δ 13 C) were calculated using Equation 7.
Here 13 Rsample is the 13 C/ 12 C ratio of the respective sample, and 13 RStandard is the 13 C/ 12 C ratio of the standard. All values are in per mil (‰) by multiplying the δ 13 C values by a factor of 1,000 (Brand, 2013;Coplen, 2011). The δ 13 C values are based on the δ 13 CIAEA-603 -LSVEC scale. Samples were analyzed against a calibrated in-house-standard (Acetanilide: −30.06 ± 0.05‰). Additionally, a quality control standard (Caffeine: −40.46‰) was interspersed between samples. The precision of the sequences was smaller or equal to 0.1‰.
The leaf samples of the Quercus ilex trees were analyzed in the same way as the samples of the herbaceous layer. The only difference were the sampling dates (Table 2). For the Quercus ilex, the leaf sampling was done during the winter period when the foliar nutrient concentrations are more stable. Sampling times and numbers are shown in Table 2. For nutrient analysis branches from the upper third of the crown (both North and South orientations) were detached and 100 leaves per tree were sampled. Between 6 and 8 trees (except for 2013 when 24 trees were sampled) were sampled per treatment (

. Estimation of Fertilization Effects and Uncertainties
To quantify the changes in carbon and water fluxes after the fertilization the pre-fertilization differences and post-fertilization differences between the treatments were assessed. We identified two sources of uncertainty for this comparison that relate to (1) within-FFC heterogeneity and (2) differences in seasonality and interannual variability (IAV) between treatments ( Figure S5). These independent sources of differences and uncertainty were considered by following the subsequent steps: 1) Within-FFC heterogeneity was accounted for by calculating between-treatment flux differences at a half-hourly temporal scale. The variability in meteorological conditions and wind direction modifies the location and size of each footprint. The vegetation properties of the different footprints vary according to microtopography, tree density, and species composition of the herbaceous layer within each footprint. Since these factors are not spatially correlated between FFCs, this variability must be accounted for when assessing the overall signals of the footprint climatologies. Before fertilization, flux differences carry information about pre-fertilization differences, spatial variability, and flux uncertainties . After fertilization, the effect of nutrients could be reflected in both the average difference between fluxes and the variability within the FFCs.
To characterize the differences between treatments and the uncertainty that the intra-FFC variability causes in the determination of such differences, the half-hourly fluxes in the unfertilized, or CT were subtracted from those of the fertilized treatments (NT and NPT). For half-hourly (HH) differences for a flux (F) between NT and CT, this would be: where Δ indicates the difference between a treatment and the control, F is any flux of interest (e.g., GPP or ET), and i stands for any individual half-hour. Due to the proximity of the treatments, the meteorological conditions (e.g., atmospheric pressure, precipitation, shortwave-and longwave incoming radiation, wind speed, wind direction, air temperature) are virtually identical for all treatments (data not shown). Therefore, the distribution of Δ i F within a particular time window (e.g., spring (Δ i f F)) holds information about the mean (  ΔF i ) flux differences between the control and fertilized treatments and the variability within the FFCs. The associated uncertainty ( ) is determined as the effective standard error ( ΔF i ): with  ΔF i being the standard deviation of   Δ i f F and eff n the effective number of observations. eff n accounts for the autocorrelation of the time series and was estimated using the first coefficients of the empirical autocorrelation function, as described by Wutzler et al. (2019). The empirical autocorrelation function was estimated using continuous and, thus, gap-filled time series. The gap-filled time series were smoother than the time series that excluded gaps; thus, both the autocorrelation and the resulting uncertainty were probably conservatively overestimated.
2) In our study, the fluxes measured in 2014 served as the pre-fertilization reference where all treatments were measured in parallel and, thus, allowed us to estimate the between-treatment differences for each season. The use of only 1 year as the baseline to determine the treatment effect might be hampered by the fact that flux magnitude varies between different years and seasons according to meteorological conditions, and proper characterization of the differences in the interannual variability of fluxes before fertilization is not possible.
To overcome such a limitation, we used the Landsat NDVI time series for the period 1999-2019 calculated for each FFC. Optical remote sensing is sensitive to vegetation, structural, and biochemical properties, including LAI and chlorophyll content (Baret & Buis, 2008;Homolová et al., 2013), and is, therefore, suitable to quantify the differences and variability of vegetation properties in different FFCs. We estimated the linear relationship from NDVI between the control and fertilized plots. The long-term NDVI relationship between the FFC of the treatments was considered as the average longterm pre-fertilization relationship between the treatments. Since NDVI values within the FFC are expected to have similar variance and errors, reduced major axis regression was performed using the "lmodel2" package version 1.7-3 (Legendre, 2018;Legendre & Legendre, 1998) in R (R Development Core Team, 2015).
Because NDVI is a good proxy of greenness and relates to ecosystem productivity, the relative deviation between the 2014 NDVI values and the long-term linear regression, can be interpreted as the relative pre-fertilization uncertainty originating from using only 2014 as a reference year instead of using a long-term average. The resulting uncertainties were 4.6% and 3.3% for CT/NT and CT/NPT, respectively. Using the assumption that the relative error of a flux F was similar to the relative error of NDVI, the corresponding pre-fertilization uncertainty for NT was calculated as: As the relative pre-fertilization uncertainty is derived from a linear relationship, it depends on the value of the mean differences (i.e.,  ΔF Pre i ). Following Equation 8-10 and assuming the error components to be independent, the total fertilization effect (FE) of F in a particular post-fertilization year and season and its associated uncertainty ( unc FE ) is: The same procedure was followed for all spring and winter periods, years, and flux variables (Figure 4, S5). The individual seasons were selected to cover their key features. The spring periods last for 56 days each and capture the peak of the growing season. The winter periods last for 61 days and cover the lowest seasonal temperatures. Table 3 shows the exact starting and ending dates for each season used in the analysis and in Figure 1 the respective meteorological conditions and flux magnitudes can be seen. Due to the variability in meteorological conditions the starting and end dates of each season change from year to year but they always have the same length.

Seasonal and Annual Sums
To estimate the impact of fertilization on carbon and water fluxes and the WUE of the ecosystem, we calculated annual and seasonal cumulated sums (Table 4). The differences between the fertilized and control treatments were then used to calculate the fertilization effect on the annual and seasonal sums. We emphasize that differences between the fertilization effects from the high-quality half-hourly data and the gapfilled cumulated sums are expected. The main reason is the larger fraction of good quality EC data during daytime hours compared to nighttime data.

Functional Relations and Regression Analyses
To characterize the changes of functional properties due to fertilization, we analyzed the relationships between LAI and LAI * [N] with A max , GPP sat , and G s for each treatment, individually. The goodness of fit allows us to evaluate whether leaf [N] partially explains the respective relationships or not. LAI and leaf [N] were based on spatial sampling within the FFC of the treatments, and the values A max and GPP sat were time averages for periods of ±7 days around the sampling date as individual values were only calculated every second day. The G s values were averaged for the same period but were based on midday values between 12:00 and 14:00 UTC. Linear, nonlinear, and segmented regressions were fitted to the data.
Linear and nonlinear least squares regressions were performed in R using the "lm" and "nls" functions (R Development Core Team, 2015), respectively. The nonlinear regressions used the following formula to EL-MADANY ET AL.  The averages are based on the 5 years of the post-fertilization period between March 2015 and February 2020. The displayed fluxes are net ecosystem exchange (NEE), gross primary productivity (GPP), ecosystem respiration (Reco), evapotranspiration (ET), transpiration (T), and evaporation (E) including their respective uncertainties (± unc). All carbon fluxes include the mean uncertainty based on u* filtering and the uncertainty of marginal distribution sampling for NEE. For ET, the flux uncertainty was estimated based on marginal distribution sampling. All uncertainties were calculated for each season and year and then propagated based on standard error propagation. The fertilization effects (Δ) were calculated as the difference between the fertilized and control treatment including the respective standard deviations (± sd) to show the year to year variability of the fertilization effects. Additionally, the relative change compared to the control is given as percent. Table 4 Average

Seasonal (Spring, Winter) and Annual Carbon (C) and Water (H 2 O) Fluxes for the Control (CT), Nitrogen (NT), and Nitrogen + Phosphorus Treatments (NPT)
characterize the saturating behavior of a relationship in which the parameters a, b, and c were optimized to minimize the residual sum of squares (RSS) between model and observations: Here, y is the property of interest (e.g., GPP sat ), and x is either LAI or LAI * [N]. The starting values for a, b, and c were the same for all treatments (i.e., LAI * [N] vs. GPP sat ) in the CT, NT, and NPT plots.
In addition to the nonlinear regressions, segmented regressions were performed using the "segmented" package (Muggeo, 2003(Muggeo, , 2017 in R to estimate the breakpoint between the increase and plateau of the nonlinear relationship. No prior value was provided to the breakpoint estimation algorithm. When comparing the means of two or more distributions a two-sample t-test was used. The "t.test" function in R (R Development Core Team, 2015) was used. Each distribution was tested against the others and the significant differences between means were encoded in letters to show which distributions means were significantly different from each other.
Analysis of Covariance (ANCOVA) was used to evaluate the treatment effect on changes in the linear relationship between variables. The "aov" function in R was used.

Leaf Level Nutrients
The increase in [N] in the herbaceous layer of both fertilized treatments, NT and NPT, was most significant in the first 2 years following fertilization (2015 and 2016). In 2017, it was higher for only a few sampling periods ( Figure 3). The differences in [N] between the fertilized and control treatments became smaller from 2015 to 2019.
[P] increased significantly after fertilization for all years in the NPT plot, except during the first sample in 2015. While both nutrients were taken up by the herbaceous vegetation, the resulting increased [N] reduced faster than [P]. The comparison of N:P ratios between CT, NT, and NPT showed that fertilization modified the nutrient stoichiometry of the herbaceous vegetation as expected, with a higher average N:P ratio (stoichiometric imbalance) in the NT plot (9.3 ± 2.5; mean ± s.d.) than in the CT (7.3 ± 0.8) and NPT plots (7.5 ± 2.4) (Figure 3). In contrast, the trees did not consistently vary in leaf nutrient content during the post-fertilization period (Figure 3).

Carbon Fluxes
The impact of fertilization on carbon fluxes was largest for the spring and winter periods of 2015 and 2016 (Table 4 and Figure 4). GPP generally increased, while Reco increased much less or even reduced, which caused NEE to become more negative. For the spring periods from 2015 to 2019, the average HH daytime ΔNEE was more negative in the fertilized treatments (−0.0156 ± 0.0076 gC m −2 (mean ± s.d.) at NT and −0.0161 ± 0.0077 gC m −2 at NPT). This resulted in 57% and 48% more negative cumulated spring NEE in the NT and NPT, compared EL-MADANY ET AL.  to CT (−46.8 ± 6.0 gC m −2 ; Table 4). Spring ΔGPP was higher in the NPT plot, +43.6 ± 18.7 gC m −2 (14%) and +39.4 ± 7.0 gC m −2 (13%) in the NT, while ΔReco was +12.7 ± 11.0 gC m −2 (5%) and +21.0 ± 19.7 gC m −2 (8%) relative to the CT (Table 4).
The average cumulated winter NEE became 82% (NT) and 110% (NPT) more negative in the fertilized treatments compared to the CT (−12.4 ± 13.6 gC m −2 ). Again, this increase in carbon uptake was driven by a substantial increase in GPP in both fertilized treatments (i.e., ΔGPP +19.0 ± 16.3 gC m −2 (13%) in the NT and +19.6 ± 18.9 gC m −2 (14%) in the NPT) compared to the increase in Reco.
When comparing drier (2015, 2017, 2019) and wetter (2016, 2018) spring periods, the fertilization effects vary in their strength and between treatments. The fertilization effects in dry and wet spring periods ( Figure S6) were in opposite directions in the NT and NPT plots and point toward stronger carbon uptake in the NPT during drier conditions and in the NT during wetter conditions. At the annual scale, ΔGPP was positive in the fertilized treatments +131.6 ± 31.2 gC m −2 (12%) in the NT and +143.2 ± 46.7 gC m −2 (13%) in the NPT, but ΔReco was less high +49.6 ± 45.4 gC m −2 (NT) and +68.5 ± 54.4 gC m −2 (NPT), respectively. As a result, the annual NEE changed from an average carbon source in the CT (+75.0 ± 20.6 gC m −2 ) to carbon-neutral in the NT (−7.0 ± 18.5 gC m −2 ) and NPT (+0.4 ± 22.6 gC m −2 ) plots.

Water Fluxes
ET was highest during the spring and lowest during winter with cumulated values at CT of 122.7 ± 2.2 kgH 2 O m −2 and 29.4 ± 1.1 kgH 2 O m −2 , respectively. The fertilization effects on water fluxes differed between the NT and NPT. While changes in ΔET at NPT during spring and winter were on average negative, they were positive in the NT (Table 4, Figure 4). During spring and winter, the cumulated seasonal sums of T increased for both fertilized treatments, and they decreased for E. ΔT was highest in the NT (+12.4 ± 6.3 kgH 2 O m −2 ; 13%) while in the NPT, it was only +4.4 ± 3.8 kgH 2 O m −2 (5%). In addition, E was substantially reduced in the NPT plot compared to the NT plot, with an overall decrease in ET in the NPT plot (Table 4).

Water-Use Efficiency
Changes in carbon and water fluxes due to fertilization display increased uptake of carbon during the growing season for both treatments and a reduction in ET in the NPT in spring, while ET and T increased in the NT for the same period (Table 4). A comparison across three measures of WUE (G 1 , IWUE, and uWUE) showed that they all follow the same main patterns as a result of fertilization ( Figure S7). Changes in G 1 , as well as G s , T, and GPP during the post-fertilization spring periods, are reported in Figure 5. Overall, strong interannual variability with consistent treatment effects was observed.
On average, NPT had the lowest G 1 values, which means it had the highest WUE, followed by NT. Both fertilized treatments showed increased G s and T compared to the CT. However, for both, G s and T values of the NPT were significantly lower in most springs compared to the NT (p < 0.01; Figure 5). As seen in Table 4 there was less evapotranspiration and transpiration in the NPT relative to the NT plot. The values of GPP were always (except for NPT in Spring 2016) significantly higher in the NT and NPT plots (p < 0.01) than in the CT plot ( Figure 5). Overall, the increase in WUE of the fertilized treatments, compared to the CT, was driven by the increase in GPP due to the fertilization. However, the loss of water shown by T mediated the magnitude and differences between the NT and NPT. The different behavior of GPP and G s between the fertilized treatments resulted in different WUEs. In addition, in three out of five and four out of five post-fertilization years, G s and T were significantly (p < 0.01) lower in the NPT than in the NT, respectively, indicating a better water-saving strategy in the NPT plot.
The δ 13 C isotopic signatures from the plant samples of the herbaceous layer show similar patterns as the eddy covariance based WUE estimates. While we did not find significant differences between the treatments in the control year (2014), we observed significantly higher δ 13 C values of the herbaceous layer in the NPT (p < 0.05; Figure 6). In contrast, the isotopic signatures of the herbaceous layer in the NT plot were significantly higher in the first year than in the CT plot and, thus, underline the more substantial increase in WUE in the NPT compared to the NT plot.
The δ 13 C values agree with the average post-fertilization WUE T estimates. WUE T spring values were 3.327 (CT), 3.307 (NT), and 3.624 gCkgH 2 O −1 (NPT) and show the highest WUE in the NPT and similar values for the NT and CT plots.

Functional Relationships
Strong linear relationships between LAI green and A max [R 2 = 0.80 (CT), 0.93 (NT), 0.92 (NPT)] across all treatments were observed (Figure 7a). The offsets in the linear fits ranged between 10.9 and 14.3 µmol m −2 s −1 . Also, both fertilized treatments showed higher A max values than the CT after fertilization.
GPP sat values collapsed into a small range between 5 and 7 µmol m −2 s −1 for the sampling points of July when LAI green values across treatments were below 0.5 m 2 m −2 and water stress was high (Figures 7b and 1). The fit of LAI green and GPP sat was not linear and GPP sat flattened at higher LAI green values (Figure 7b). The change point of the segmented regression analysis showed that the plateau started at a LAI green of about 1.8 ± 0.15 m 2 m −2 . Similar to A max , GPP sat values in the CT were lowest, showing a lower photosynthetic capacity in the CT compared to the fertilized treatments. Figure 5. Between-tower differences for water-use efficiency-related parameters for the spring periods from 2015 to 2019. Boxplots for each treatment [control (CT, purple), nitrogen (NT, dark blue), and nitrogen + phosphorus treatments (NPT, light blue)] and year are the differences between the respective post-fertilization years (2015-2019) and the pre-fertilization period of 2014. The panels display the stomatal slope parameter (G 1 , (a)), surface conductance (G s , (b)), canopy transpiration based on the TEA algorithm (c), and gross primary productivity (GPP, (d)). Letters above the boxplots indicate significant differences in the means based on two-sample t-tests at the p < 0.01 level.
The G s were nonlinearly related to LAI green similar to GPP sat , and the treatments reached their plateau around 6 mm s −1 (Figure 7c). The highest values of G s were detected in the NT, while NPT and CT had lower and more similar values. The largest differences between NT and the other fitted curves were detected for LAI green values below 1.1 m 2 m −2 . The G s values converged with increasing LAI green , especially for values after the change point of the segmented regression.
When analyzing the functional relationships above with LAI green * [N] the R 2 values for LAI green  Based on linear relationships between G s and GPP sat , the fertilized treatments have higher GPP sat for the same G s values that is, a significant difference in offset but no significant difference in slope. On the other hand, for the same increase in GPP sat , the increase in G s is highest in the NT plot but not significantly different from the NPT and CT plots ( Figure S8    Functional relationships between maximum net carbon uptake (A max ), photosynthetic capacity (GPP sat ), and surface conductance (G s ) with green leaf area index (LAI green , top row) and green leaf area index * nitrogen concentration (LAI green * [N], bottom row) of the herbaceous layer for the control (CT, purple), nitrogen (NT, blue), and nitrogen + phosphorus treatments (NPT, light blue). All figures include the respective regression formulas and fits (solid lines). In addition, the nonlinear fits contain a segmented regression (dashed black line) with the respective breakpoint of the two linear fits and its uncertainty (black dot with dashed line at the bottom). Horizontal error bars display the standard deviation from the LAI spatial sampling, and vertical error bars display the temporal variability of A max , GPP sat , and midday G s for ±7 days around the sampling day. For the nonlinear fits, the residual sum of squares (RSS) was used to show the goodness of fit.

Discussion
We will discuss the impacts of nutrient availability and stoichiometric imbalance on carbon and water fluxes, ecosystem-scale WUE, functional relationships, canopy structure, and nutrient content of the vegetation.
The observed changes in leaf nutrients within the herbaceous layer, but not in tree leaves, suggest that the ecosystem response to the fertilization was driven by the herbaceous layer rather than by the trees. The higher response of the herbaceous layer could be related to the common competition between trees and grasses for resources in savanna ecosystems (Riginos, 2009;Rivest et al., 2011;Sankaran et al., 2004;Scholes & Archer, 1997). The marked differences in root system profiles between the herbaceous layer and trees in the study system (Rolo & Moreno, 2012) would have allowed a higher uptake of nutrients by herbaceous species in the uppermost soil layers and hampered any response in trees. This agrees with the findings of Scalon et al. (2017), who analyzed the leaf nutrient concentrations of woody species in a Brazilian savanna in a long-term nutrient manipulation experiment and found no influence of N and P additions to the leaf nutrient concentrations. In contrast, a nutrient manipulation experiment conducted in the Harvard experimental forest for 15 years showed a constant increase in leaf nutrient content for oak and pine trees (Magill et al., 2004). The absence of increased leaf nutrients in the trees with simultaneously increased leaf nutrients of the herbaceous layer allows us to suggest that the ecosystem response to fertilization was driven by the herbaceous layer rather than by the trees. With only leaf data available, we cannot rule out that the trees have responded to the higher N availability different than increased leaf N (e.g., by increased LAI or storing N in other organs).
The strong response in carbon uptake during the first two winters following fertilization (Figure 4) could be a result of the strong mobility of N within the ecosystem. This was emphasized by the reduction of [N] during the experiment. Also, the highest N availability in the soil is expected to occur during the fall with the re-wetting of the soil after the summer drought (Joffre, 1990;Morris et al., 2019). This increased N availability then caused a faster development of the vegetation biomass (Luo et al., 2020) and higher gross and net carbon uptake rates. At the same time, the high mobility of N, due to the strong precipitation in the fall, likely leads to strong uptake via the re-greening vegetation and strong leaching of N into deeper soil layers where the herbaceous vegetation cannot access it anymore, but trees potentially could. This might lead to a continuous decrease in leaf [N] in the herbaceous layer and the strong fertilization effect being restricted to the first two winter periods. In contrast, P, which is less mobile, stayed in the system in the NPT plot and led to persistently high levels of P in the herbaceous layer at NPT.
The increased N availability can alter N:P stoichiometry in plants and result in P limitation (Du et al., 2020;Li et al., 2016;Peñuelas et al., 2013). Under these conditions, plants try to increase their P uptake (Güsewell, 2004). In a recent review, Oldroyd and Leyser (2020) describe how plants sense P availability and how they react to conditions of P limitation. The main root adaptations are increased lateral growth and elongated root hairs to more efficiently extract P, which has the highest availability in the topsoil. Nair et al. (2019) analyzed root biomass and root length density in the CT, NT, and NPT plots and found increased root biomass but no increase in root length in the NT plot; while in the NPT plot, the root mass did not increase, but the root length density did, which is in agreement with the findings of Oldroyd and Leyser (2020).
Plants increase transpiration to better extract nutrients, especially P, from the soil (Cernusak et al., 2011;Cramer et al., 2008;Huang et al., 2017;Kröbel et al., 2012;Pang et al., 2018;Rose et al., 2018). Our results agree with this finding and showed that both treatments increased their LAI, but transpiration increased most in the NT plot. A consequence of increased transpiration under P limitation is the reduction in WUE. G 1 and WUE T showed that WUE increased more in the NPT plot than in the NT plot, even though GPP increased similarly during the spring (Table 4). This is supported by the higher δ 13 C at NPT, which indicates that the ratio of leaf-internal CO 2 concentration (c i ) to ambient CO 2 concentration (c a ) is smaller (Cernusak et al., 2013;Seibt et al., 2008). The decrease in the c i /c a points toward stronger stomatal closure, less transpiration and thus higher WUE (Medlyn et al., 2011). The increase in WUE at NPT is consistent between methods and across scales that is, δ 13 C at leaf level and G 1 and WUE T at ecosystem scale. This highlights the importance of P availability for WUE across scales.
The N fertilization induced a P limitation in the NT plot and the increasing transpiration might (i) increase P uptake directly via a direct water stream into the plant or by (ii) allowing increased photosynthesis and fueling soil microbes to extract P from organic matter through rhizosphere priming (Kuzyakov, 2002). This is supported by the recent work of Chen et al. (2020), who showed that during the first years following N fertilization, phosphatase activity increases and indicates a progressive attenuation of the N-induced P limitation through the plant and microbial activity. This would also be supported by the increased root biomass at NT . Additionally, the increased LAI green in the NT plot combined with similar [P] as in CT shows that more P was incorporated into plant material. This additional P at NT must be extracted from the soil or recycled from organic material, maybe through one of the suggested processes.
When sampling the abovementioned differences between treatments at the ecosystem scale, we observed a substantial IAV in WUE and its associated water and carbon fluxes (Figure 4) between dry and wet spring periods ( Figure S6). This indicates that ecosystem responses to nutrient availability and stoichiometric imbalance are modulated by water availability. Lower WUE in the NT plot, compared to NPT, resulted in faster senescence toward the end of spring, as shown by Luo et al. (2020). A reduced growing season length and lower WUE (at NT) potentially reduces the carbon uptake compared to an ecosystem losing less water with a higher WUE (NPT). Meteorological drivers, especially the ones related to water availability, strongly impact the productivity of this ecosystem during the spring and winter Sippel et al., 2017) and influence the phenology of herbaceous plant species at the site (Luo et al., 2018). On top of the IAV, there was a strong seasonal variability with the strongest fertilization effects occurring when conditions were promoting the development of the vegetation (i.e., the winter and spring), while the fertilization effects and their IAV strongly reduce during the summer (Figure 7). This, again, emphasizes that the herbaceous layer drives the fertilization effects rather than the trees, which are the only active component of the ecosystem in the summer and they only have a canopy cover of about 20%. The relationship between LAI green and GPP sat and G s especially underpin that for low LAI green conditions (LAI < 0.5 m 2 m −2 ), the treatment differences are also low.
The strong linear relationship of LAI green to A max and LAI green * [N] to A max suggests that the increased maximum net carbon uptake was driven by increased LAI in both treatments. The effect of nitrogen was through the release of the N limitation which promoted the development of more biomass and potentially, also the increase of chlorophyll content in the leaves. As shown in Luo et al. (2020), the increased biomass development was already happening during the re-greening when the fertilized treatments developed and increased their LAI more rapidly than the CT. This also explains why the largest treatment effects in carbon fluxes were found during the winter and spring periods and emphasize that the study site was N limited.
A positive relationship between leaf nitrogen content and photosynthetic capacity is described at leaf level (Evans, 1989;Kattge et al., 2009) but also at ecosystem scale (Kergoat et al., 2008;Musavi et al., 2016). When using the product of leaf nitrogen content and LAI (i.e., LAI green * [N]) and photosynthetic capacity (i.e., GPP sat ) the relationship between LAI green * [N] and GPP sat improves for all treatments (with more explained variance) than that involving only LAI green . The flattening of the curve at higher LAI values could be a consequence of lower leaf levels getting less light. Additionally, leaf properties like the maximum rate of Rubisco carboxylase activity (Vc max ) might result in light saturation. This would be in agreement with the relationship of LAI green and LAI green * [N] with G s , which reached its maximum value of approximately 7 mm s −1 after LAI green reached approximately 1.0 m 2 m −2 . The increase of mean springtime LAI green in the fertilized treatments was approximately 0.7 m 2 m −2 , which corresponds with an increase of 40% compared to CT (1.72 m 2 m −2 ). The reduced evaporation and increased transpiration during spring could be driven by increased LAI because it increases the transpiration-surface per m 2 ground area and additionally intercepts more light which is then not able to heat the soil below the leaves (Pott & Hüppe, 2007;Wei et al., 2017). The surface conductance reaches its maximum already at an LAI of 1.0 m 2 m −2 and at a lower value as compared to GPP sat (1.8 m 2 m −2 ). With increasing LAI, leaves could start to shade other leaves below, reduce their surface temperature and thus VPD within the herbaceous canopy. Similarly, the increasing LAI results in reduced aerodynamic conductance ( Figure S9) which in turn hinders the exchange of moisture out of the canopy and thus potentially reduces VPD. In addition, these processes are modulated by modified WUE as a result of nutrient availability and stoichiometric imbalance of the N:P ratio.

Conclusion
Systematic differences in ecosystem responses were detected between NT and NPT as a consequence of nutrient application. The fertilizer application at NT led to a wide N:P ratio (i.e., a stoichiometric imbalance and potential P limitation in the herbaceous layer), which increased transpiration to extract more P from the soil. In addition, springtime biomass production and carbon fluxes increased compared to the control. Altogether, springtime WUE increased in the NT plot, but not as strongly as in the NPT plot, where not only GPP increased but also the transpiration was reduced compared to NT. WUE T only increased in the NPT plot but decreased in the NT plot.
The functional relationships between LAI green and A max and GPP sat show better fits when using LAI green * [N] and emphasize the importance of [N] for explaining the carbon fluxes, while for G s using LAI green * [N] even reduced the goodness of fit. A lower N:P ratio at CT and NPT resulted in the lowest surface conductance while it was strongly increased at NT with its high N:P ratio.
Overall, the changes in WUE in the fertilized treatments were first driven by the increase in carbon uptake through increased LAI. Second, between-treatment differences were driven by the narrow N:P ratio in the NPT plot, which reduced transpiration compared to NT and increased WUE. Further, the strength of all responses was highly modulated by water availability. Our results suggest the importance of compensating P limitations in case of N:P imbalance to increase WUE in semi-arid ecosystems.

Data Availability Statement
The half-hourly eddy covariance, meteorological, nutrient, isotopic signatures, and green LAI data are made available via http://doi.org/10.5281/zenodo.4453567 . For the availability of the hyperspectral images we refere to  and the data availability statement within. of the MaNiP project, M. Pilar Martín and Vicente Burchard-Levine thank the Spanish Ministry of Economy and Competitiveness for financing the flights with hyperspectral imagery through the FLUXPEC project (CGL2012-34383) and SynerTGE (CGL2015-G9095-R), and Javier Pacheco-Labrador and Mirco Migliavacca acknowledge the German Aerospace Center (DLR) and the German Federal Ministry of Economic Affairs and Energy that provided support within the framework of the EnMAP project (Contract No. 50EE1621). Victor Rolo was funded by a "Talento" fellowship (TA18022) funded by the regional government of Extremadura (Spain). Open access funding enabled and organized by Projekt DEAL.