Altered seasonal sensitivity of net ecosystem exchange to controls driven by nutrient balances in a semi-arid savanna

Semi-arid ecosystems dominate the variability and trend of the terrestrial carbon sink. They are sensitive to environmental changes following anthropogenic influences, such as an altered ratio of nitrogen (N) to phosphorus (P) due to increasing N deposition. Semi-arid savannas with different vegetation compositions have complex carbon flux dynamics, and their responses to environmental change are not yet well understood. We analyzed a long-term (2016–2022/2023) dataset of flux, biometeorological, and vegetation data (satellite and ground measurements) from a manipulated semi-arid savanna to examine how altered nutrient levels and stoichiometric balance affect the seasonal sensitivity of net ecosystem exchange (NEE) to its drivers. We used singular spectrum analysis to extract the seasonal signals of all variables and assessed the key drivers of NEE over the study period as a whole and in different seasons, using Pearson correlation and information theory. We found that the addition of both N and N + P to the ecosystem increased seasonal NEE variability, driven by greenness of the herbaceous layer. Analyzing 7 years of combined data, we found that the water limitation in summer and the energy limitation in winter outcompeted the fertilization effect. We learned that the effects of nutrient addition on NEE–control relationships became clearer through analyzing different phenological seasons. In the summer, N + P addition led to a potential change in species composition and productivity, resulting in a stronger interaction between the herbaceous layer and NEE. During the transitional seasons (i.e., drydown and autumn), which mark the senescence and regreening of the herbaceous layer, we found NEE to be less sensitive towards meteorological drivers like relative humidity, radiation, and air temperature with N addition. In the future, the increasing NEE variability might become even more pronounced with increasing N deposition and a changing climate.

1 Introduction

Terrestrial ecosystems play a major role in the global carbon cycle, with the ability to store significant amounts of carbon (Friedlingstein et al., 2022). While forests and wetlands contribute most to the terrestrial carbon sink, semi-arid ecosystems dominate its trend and interannual variability (Ahlström et al., 2015; Poulter et al., 2014; Zhang et al., 2016). Semi-arid ecosystems typically take up carbon from the atmosphere during the wet season and are either dormant or emit carbon during the dry season (Metz et al., 2023). Net ecosystem exchange (NEE) describes this balance between carbon uptake through photosynthesis, typically expressed as gross primary productivity (GPP), and carbon release through ecosystem respiration (R_eco). NEE in semi-arid regions varies strongly from year to year, depending on the climatic conditions and water availability (Haverd et al., 2017; Piao et al., 2020).

Despite their important role in the global carbon cycle, semi-arid ecosystems and their dynamics are still not well understood. Long-term in situ measurements from these regions are scarce. Particularly, eddy covariance (EC) measurements, which provide high-frequency and continuous ecosystem trace gas and water flux data (Baldocchi, 2020), are underrepresented in these regions (Jung et al., 2020). Consequently, semi-arid ecosystems remain poorly represented in terrestrial biosphere models (Fawcett et al., 2022; MacBean et al., 2021) due to their complex structure and high spatio-temporal variability, which is difficult to generalize.

Recently, efforts have been made to examine drivers of NEE in semi-arid savannas in order to better understand their role in the global carbon cycle (Baldocchi and Arias Ortiz, 2024; Kannenberg et al., 2024; Ma et al., 2007, 2016; Zhang et al., 2010). Water-related variables like precipitation and soil moisture availability are amongst the main NEE drivers (Archibald et al., 2009; Baldocchi and Arias Ortiz, 2024; Del Grosso et al., 2018; Huang et al., 2016b; Morgan et al., 2016), as they promote plant photosynthesis (Parton et al., 2012) and enhance heterotrophic respiration rates (Ma et al., 2016). Furthermore, photosynthetically active radiation (PAR), vapor pressure deficit (VPD), and air temperatures can strongly impact NEE (Archibald et al., 2009; Baldocchi and Arias Ortiz, 2024; Del Grosso et al., 2018). Other biotic factors, like soil microbial communities and organic matter, also play an important role in the ecosystem carbon cycle and contribute to R_eco (Austin and Vivanco, 2006; Bastida et al., 2016; Hu et al., 2014). These drivers can differ for different vegetation types. In semi-arid regions, savannas are a typical ecosystem type. They comprise coexisting vegetation layers (e.g., tree and grass), which interact in complex ways (Higgins et al., 2000; House et al., 2003). The layers differ in their rooting depths (Moreno et al., 2005; Rolo and Moreno, 2012), water use strategies (Cubera and Moreno, 2007; Miller et al., 2010; Steiner et al., 2024), and phenological and life cycle strategies (Whitecross et al., 2017). The herbaceous vegetation in two-layer ecosystems is often underestimated in its importance for the ecosystem water and carbon fluxes (Dubbert et al., 2014). On the Iberian Peninsula, dehesas (or montados in Portuguese) – savanna-like agroecosystems shaped by long-term human management – are widespread (Den Herder et al., 2017). Dehesas are open oak woodlands with an herbaceous layer that consists mainly of annual grasses and sometimes crops. The tree layer is evergreen (Moreno, 2008), whereas the herbaceous layer typically follows an annual cycle of growth, senescence, and regreening (Ma et al., 2007; Perez-Priego et al., 2015). As savannas are typically characterized by changing resource limitations throughout the year (Luo et al., 2020; Ries and Shugart, 2008), limited by water in the dry season and by nutrients and energy in the wet season (Moreno, 2008; Morris et al., 2019; Nair et al., 2019; Whitley et al., 2011), the dominant environmental drivers change with the seasons. The complex interactions between tree and grass layers, the dominant environmental drivers of NEE, which change with the seasons, along with changing limitations, result in highly complex ecosystem carbon flux dynamics, which remain to be fully understood.

Semi-arid ecosystems face numerous human-induced environmental changes, including stoichiometric imbalances between nitrogen (N) and phosphorus (P). These imbalances arise from increasing N inputs into ecosystems due to fertilizer use and the combustion of fossil fuels (Steffen et al., 2015) without a corresponding increase in P inputs (Penuelas et al., 2013). Few studies so far have dealt with the impact of altered nutrient levels on NEE and its drivers in semi-arid regions. The availability and stoichiometric balance of N and P influence ecosystem functioning and plant traits (Reichstein et al., 2014), water use efficiency (El-Madany et al., 2021; Huang et al., 2016a), canopy structure (Migliavacca et al., 2017), composition of species (Sardans et al., 2012), and the seasonality of vegetation activity (Luo et al., 2020). However, different plant types react differently to changes in nutrient availability, due to variations in generation times and buffering capacities (Pardo et al., 2011). Therefore, the understanding of the response of complex tree–grass ecosystems to changes in N and P availability and their stoichiometric balance is still poor.

In this study, we took advantage of the unique long-term dataset collected in a semi-arid dehesa, Majadas de Tiétar, in southwestern Spain. A large-scale nutrient addition experiment has been running here since 2015, providing an exceptional opportunity to study the long-term influence of altered N : P ratios on ecosystem functioning (El-Madany et al., 2021). Three EC flux towers have been set up, with the footprint of one tower receiving N fertilization, another one receiving N + P fertilization, and the third serving as control. Previous studies found that both treatments increased the annual carbon uptake of the ecosystem and that N + P addition increased the water use efficiency of the ecosystem more than N-only addition, which could be attributed to higher transpiration rates and a changed root strategy in the N-only-fertilized plot (El-Madany et al., 2021; Nair et al., 2019). Nutrient addition led to not only a higher seasonal amplitude of maximum GPP and a faster increase in GPP during the regreening period, but also a faster senescence during the drydown period, indicating changes in plant structure and physiology (Luo et al., 2020).

Here, we analyzed a 7-year (2016–2022/23) time series of daily values of environmental and biogenic variables from Majadas de Tiétar, combining flux data, meteorological measurements, digital repeat photography, and satellite data to address the following question: how do altered nutrient levels and stoichiometric balance affect

• annual NEE and its variability,

• the relationship between NEE and its key controls,

• the relationship between NEE and its key controls in different seasons, and

• the sensitivity of NEE to its controls over time?

The relationships between NEE and its controls vary across different timescales (Mahecha et al., 2007). To disentangle these timescales from the time series and eliminate noise from the high-frequency measurements, we can use decomposition methods (Linscheid et al., 2020). On short timescales, the NEE sensitivity follows the diurnal cycle of the sun, showing a great dependency on radiation. Ecosystem-level responses, in contrast, often develop on scales of months, seasons, or years (Ma et al., 2016). Therefore, we extracted the seasonal signals of all variables from the time series with singular spectrum analysis, a data-driven time series decomposition method. To characterize the relationship between NEE and its environmental and biogenic controls on the seasonal scale, we used a conventional correlation metric: the Person's correlation coefficient. As these relationships can be quite complex and nonlinear, we additionally used information-theory-based methods. These non-parametric metrics, such as mutual information, can identify non-monotonic couplings as well as leading and lagging effects and are therefore particularly suited for driver identification in observational flux data from complex systems (Chamberlain et al., 2020; Knox et al., 2018; Sturtevant et al., 2016). As NEE controls vary in their importance throughout the year due to a high seasonality of the ecosystem, we split the dataset into phenological seasons defined by vegetation responses to different limitations.

2 Material and methods

2.1 Site description

The Majadas de Tiétar research site is located in western Spain (39°56^′25^′′ N 5°46^′29^′′ W) (Fig. 1a). The local ecosystem consists of an herbaceous stratum and scattered evergreen oak trees (98 % Quercus ilex). The tree density is around 20–25 trees per hectare, with a mean diameter at breast height of 46 cm (El-Madany et al., 2018). The fractional canopy cover of trees is 23 % and canopy height is on average 8.7 m (Bogdanovich et al., 2021). The tree leaf area index (LAI) is around 0.35 m² m⁻². Although the grass layer has a peak LAI in spring, it is quite spatially variable with values between 0.5 and 2.5 m² m⁻² due to the seasonal temporal variability (described in detail below) of grass growth (Migliavacca et al., 2017). The site is managed and continuously used for grazing livestock at a low density of 0.3 cows per hectare (El-Madany et al., 2018). In the driest months (July–September), the farmers move the cattle to nearby mountain grasslands (personal communication).

The climate at the site is semi-arid with an annual precipitation of around 650 mm with strong interannual variability. Almost 85 % of the annual precipitation falls in the wet season between October and April, whereas the rest of the year is dry with occasional rains (Fig. 1c) (El-Madany et al., 2021). According to Nair et al. (2024), we defined five different seasons. Spring is the main growing season and usually starts around March and ends in late May. Then the drydown period starts and the grasses start to become senescent due to depletion of soil moisture, increase in temperature, radiation, and vapor pressure deficit. The summer (typically between end of June until end of September) is characterized by long-lasting dryness and a dormant/dead herbaceous layer. With the onset of precipitation (usually in October), the autumn starts and the herbaceous layer regreens (Nair et al., 2024). The winter months (December–February) are energy-limited. The onset and offset of the different seasons vary from year to year, depending on water availability and winter temperature (Luo et al., 2020). The mean annual temperature is 16.7 °C with an average minimum temperature of −4.7 °C and maximum temperature of 41.1 °C (2004–2019) (El-Madany et al., 2021). Dominant wind directions are west-southwest and east-northeast (El-Madany et al., 2018).

Three EC towers at ecosystem level were operated simultaneously at the site (Fig. 1b). The ecosystem is heterogeneous with a high variability in plant species in the herbaceous layer (at a scale of centimeters) and tree cover (at a scale of meters). It becomes homogeneous on the scale of a few hundred meters. The daytime flux footprints of the three towers are homogeneous at the relevant scale and do not overlap with each other under prevailing meteorological conditions (El-Madany et al., 2018). The control tower (Fluxnet ID: ES-LMa) has been operated since 2003; it is hereafter referred to as CT. The north tower (Fluxnet ID: ES-LM1) was set up at a distance of 450 m from CT in the northwestern direction and the south tower (Fluxnet ID: ES-LM2) was located 630 m in the southern direction from CT (El-Madany et al., 2018). Since 2015, a large-scale fertilization experiment has been conducted at the site, where N fertilizer is added in the footprint of the north tower (hereafter referred to as NT) and N and P fertilizer are applied in the footprint of the south tower (hereafter referred to as NPT) (El-Madany et al., 2021). The N and P fertilizations were applied around similar times at the sites each year, with some exceptions due to weather or logistics restrictions (e.g., pandemic). N was added at 100, 20, 50, 24, and 12 kg N ha⁻¹ at both sites by end of winter 2015, 2016, 2017, 2021, and 2023, respectively, and P was added at 50, 10, 25, 6, 6, and 6 kg P ha⁻¹ at NPT in autumn 2014, 2015, 2016, 2019, 2020, and 2022, respectively. This timing of the application of N and P was selected to have maximal possibility for use by vegetation in the next growing season after each addition. Next to each flux tower is a radiometric tower setup measuring radiation components above the tree and grass layers, alternating every 15 min, providing half-hourly measurements for each layer.

Figure 1 (a) Site location on the Iberian Peninsula. (b) Location of the three eddy covariance towers. Nitrogen-added tower (NT) is in blue, control tower (CT) is in purple, and nitrogen + phosphorous-added tower (NPT) is in light blue (base map from Google Earth, 2025). The tower locations were chosen to ensure that, under prevailing wind conditions, their measurement footprints do not overlap. Footprint climatologies can be found in Fig. 1 in El-Madany et al. (2018). (c) Average monthly precipitation total and average temperature measured at 15 m over 2016–2023.

2.2 Eddy covariance and biometeorological data

Each of the three EC towers continuously measures sensible heat flux (H), latent heat flux (LE), and CO₂ flux. Each system is equipped with an R3-50 sonic anemometer (Gill Instruments Limited, Lymington, UK) to measure three-dimensional wind components and sonic temperature and an LI-7200 infrared gas analyzer (LI-COR Bioscience, Lincoln, Nebraska, USA) to measure CO₂ and H₂O mixing ratios. The measurement heights are the same at 15 m above ground (El-Madany et al., 2021). The flux and meteorological data were collected as described by El-Madany et al. (2018). The raw high-frequency data were processed with EddyPro v.7.0.9 (Fratini and Mauder, 2014). The post-processing was done in R using the REddyProc package (Wutzler et al., 2018). The storage corrections of the CO₂ flux were made with profile measurements from seven points on the EC towers. A friction velocity (ustar) threshold was applied following Papale et al. (2006) and data with ustar values below the defined threshold were removed. Missing and low-quality data were gap-filled (Mauder and Foken, 2011; Reichstein et al., 2005) for calculating the annual budgets. Additional atmospheric variables that we used are air temperature (Ta) and relative humidity (Rh) measured at two heights (2 and 15 m), vapor pressure deficit (VPD), air pressure (air_press), and friction velocity (ustar). Furthermore, we incorporated radiometric components such as longwave downward radiation (LWDR), shortwave downward radiation (SWDR), and photosynthetically active radiation (PAR). Soil measurements were comprised of soil temperature measured from two sensors in open pasture (Tsoil_op) and below oak tree canopy (Tsoil_bc) and soil heat flux from two sensors in open pasture (SHF_op) and below canopy (SHF_bc). We used the average measurements of two sensors when both were available and otherwise the measured values of one sensor. Additionally, to calculate soil water content, we used the different measurements integrated over the top 20 cm of the soil and weighted them by canopy a cover of 20 % to obtain ecosystem soil water content values (SWCn).

Some small data gaps existed in meteorological variables and were filled with the average of the remaining two towers and interpolation. However, the PAR sensor at CT had a long-term malfunction, so we used the PAR time series from NT as a substitute, since the incoming radiation was not expected to differ substantially between the towers (El-Madany et al., 2018). In addition, we calculated evaporative fraction (EF) as the ratio between LE and available energy (EF = LE $/$ (LE + H)) (Gentine et al., 2007; Tong et al., 2022). EF is a normalized measure of the surface energy partitioning and can serve as a diagnostic of vegetation water status (Nutini et al., 2014). We calculated it at half-hourly time steps from only positive LE and H values to avoid introducing extreme outliers into the analysis. EF is strongly linked to meteorological variables like soil moisture, VPD, and net radiation (Gentine et al., 2007; Tong et al., 2022), as well as to vegetation cover and LAI (Gentine et al., 2007). A full overview of analyzed variables is shown in Table 1.

2.3 Vegetation indices

We used three different vegetation indices to represent vegetation greenness: green chromatic coordinates and albedo (both indices derived from in situ data) and the normalized difference vegetation index (NDVI) (derived from satellite data). Green chromatic coordinates (GCCs) are an effective measure for describing greenness variation in semi-arid ecosystems (Luo et al., 2018, 2020). We used daily mean GCC values extracted from the images collected every 30 min by digital cameras (Stardot NetCam 5MP) which were installed at the top of each ecosystem EC tower facing north. The cameras were set up according to the protocol of the PhenoCam network (https://phenocam.sr.unh.edu/webcam/tools/, last access: 3 July 2024) and collected red, blue, and green (RGB) images (Luo et al., 2018). GCC was computed by dividing green digital numbers (G_DN) by the sum of red (R_DN), blue (B_DN), and green digital numbers (Richardson et al., 2009):

$$\begin{array}{}\text{(1)}& \mathrm{GCC}={\displaystyle \frac{G_{\mathrm{DN}}}{R_{\mathrm{DN}}+B_{\mathrm{DN}}+G_{\mathrm{DN}}}}.\end{array}$$

At each site, we selected two regions of interest in which we calculated GCC, one capturing the grass layer (gcc_gr) and one capturing the trees (gcc_tr). The data derived from RGB images can be found on the website of the PhenoCam network (IDs: ES-LM1, ES_LMa, and ES_LM2 for the NT, CT, and NPT, respectively). At each site we chose the masks GR_1000 for the grass layer and EB_1000 for the tree layer.

We further calculated albedo as the ratio of outgoing shortwave radiation to incoming shortwave radiation, measured at the radiometric tower setup at each site. We distinguished ecosystem albedo (Alb_eco), tree albedo (Alb_tr), and grass albedo (Alb_gr) to account for reflectance and to use as another proxy for vegetation greenness and water status of the plants. We used daily averages from only daytime hours (11:00–15:00 UTC)to guarantee a high solar zenith angle for reliable measurements. Furthermore, cloudy days were filtered out and only time steps where the ratio of downward radiation to extra-terrestrial radiation at the top of the atmosphere was 0.7 or more were kept (Wood et al., 2015).

Finally, we used NDVI from the FluxnetEO dataset as a proxy describing the amount and health of vegetation cover (Tucker, 1979). The dataset complements ground measurements by providing satellite-based vegetation indices, surface reflectance, and land surface temperatures for a 2 km radius around a flux site (Walther et al., 2022). We use NDVI from MODIS (Moderate Resolution Imaging Spectroradiometer) with a daily temporal resolution (Walther et al., 2022). NDVI is calculated from the normalized difference between the reflectance of near-infrared (NIR) and red-light bands (Tucker, 1979).

Table 1 Flux, meteorological, and soil variables and vegetation indices used in this study. Soil water content was calculated based on the shadow fraction estimated from the solar zenith angle and a canopy cover of 20 %.

2.4 Data analysis

2.4.1 Aggregation to daily data

To make sure that the driver identification is not confounded by gap-filling techniques based on meteorological measurements, we only use non-gap-filled, measured flux data. To ensure that there are only high-quality measured values, we selected data with quality flag = 0 (flagging policy according to Mauder and Foken, 2004). Consequently, the data coverage of the measured half-hourly time series is quite low (around 30 %) and especially heterogeneous during the nighttime. Therefore, we calculated from the biometeorological and flux data daily mean values by aggregating only daytime measurements to avoid the bias. Daytime includes only values measured after sunrise and before sunset, identified using the suncalc package in R (Thieurmel, 2017). This does not apply to vegetation indices as they were calculated as described above. GPP and R_eco were not assessed in this study as partitioning methods depend on other environmental factors that would also confound the analysis of NEE controls.

If not stated differently, the following analyses cover the 7-year period from 2016–2022 because during this time all variables were available. For the assessment of NEE variability and budgets, we utilized data spanning 8 years (2016–2023) because this extended dataset was available and incorporating additional years enhances the robustness of observed trends.

2.4.2 Time series decomposition with singular spectrum analysis

Decomposition methods assume that observed time series are composed of additively superimposed sub-signals, each shaped by different scales of variability (Mahecha et al., 2010). Consequently, the time series represents the sum of a trend, oscillatory components at various scales, and noise (Liu et al., 2022).

Here we used singular spectrum analysis (SSA) for the decomposition. SSA is entirely data-driven and non-parametric and is therefore free of the bias of function selection (Golyandina et al., 2001; Liu et al., 2022; Mahecha et al., 2007). This makes it advantageous compared to other decomposition methods like Fourier and wavelet analysis (Baldocchi et al., 2021). SSA is more flexible in grouping components of similar frequencies than wavelet decomposition (Liu et al., 2022) and able to detect aperiodic or non-harmonic sub-signals from short and noisy signals (Golyandina and Zhigljavsky, 2013; Mahecha et al., 2007). Since it is fully phase–amplitude modulated and relatively robust to instationarities of the signal mean and variance, it is suitable for nonstationary signals (Allen and Smith, 1996; Golyandina and Zhigljavsky, 2013; Yiou et al., 2000). Even fragmented time series can be handled with it, as SSA can be used for filling gaps according to the first reconstructed component, which is a low-frequency signal (Kondrashov and Ghil, 2006). This makes it particularly useful for flux data (Mahecha et al., 2007).

SSA consists of four steps: embedding, decomposition, grouping, and reconstruction. In the first step, a one-dimensional time series y(t) is embedded into a two-dimensional lagged matrix X, by shifting a moving window of a certain window length (L) along the time series. In the second step, a singular value decomposition of X is performed, and it is decomposed into its orthogonal components by determining eigenvalues and eigenvectors corresponding to principal components. The eigenvalues of the covariance matrix $\mathbf{X}\cdot {\mathbf{X}}^{\prime }$ are then ordered in decreasing magnitude. In the next step, the components are grouped, as some sub-signals consist of a set of components with complementary oscillatory frequency. In the last step, by inverting the ranked principal components, the reconstructed components of the original time series are computed. These reconstructed components show how much of the variability of the original time series is associated with the different timescales. A more detailed description of the method can be found in Golyandina et al. (2018).

Here we used the Rssa package in R (Golyandina and Korobeynikov, 2014) for our analysis. To support our hypothesis that daily-scale NEE variations are predominantly influenced by radiation, with a neglectable effect of nutrient addition, we conducted a preliminary analysis extracting the daily signals of NEE and all potential driving variables from half-hourly measurements. Detailed procedures and results of this analysis are provided in the Supplement (Text S1, Figs. S1, S2).

For our analysis, we extracted the seasonal signals of the daily time series of all variables shown in Table 1. First, as required by SSA, we gap-filled the time series with the Rssa package's internal function, igapfill, which fills gaps using the low-frequency component of the time series itself (i.e., not based on meteorological measurements). For gap filling, as a window length (L) of $n/\mathrm{2.5}$ is recommended (Mahecha et al., 2007), we selected a gap-filling window length of L = 1000 for 2557 data points from 7 years of daily data. By conducting a sensitivity analysis, we found that adding a 3-month margin at the beginning and end of a time series can help to reduce edge effects during the gap-filling process (details see Text S2).

To extract the seasonal signal, we reconstructed the components of the frequency bin from 15 to 366 (days). We selected L=732 (2 years) based on the criteria that L should be less than $n/\mathrm{2}$ and ideally an integer multiple of the period length to be extracted in order to ensure a clear signal (Biriukova et al., 2021; Golyandina and Zhigljavsky, 2013). Frequency contributions of less than 0.2 were defined as noise (Liu et al., 2022). For the grouping, we used the automated method provided by the Rssa package, which identifies groups using a hierarchical clustering algorithm based on the w-correlation matrix. The w-correlation matrix shows the weighted correlations between reconstructed components (von Buttlar, 2014; Golyandina and Korobeynikov, 2014).

For analyzing the changes in seasonal NEE variability and budgets, we used data from 2016–2023. Accordingly, L was set at 1169 ($L=n/\mathrm{2.5}$, with n=3105). To account for seasonal variability, we calculated the standard deviation of the reconstructed NEE signal for each year to quantify the variation amplitude.

2.4.3 Pearson correlation coefficient

To identify the key drivers of NEE, we first computed in R the Pearson correlation coefficients (r) between NEE and all the investigated variables (Table 1) from the reconstructed seasonal signal using the daily datasets. It is calculated as follows:

$$\begin{array}{}\text{(2)}& r={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n(x_i-\underset{\mathrm{¯}}{x})(y_i-\underset{\mathrm{¯}}{y})}{\sqrt{{\sum }_{i=\mathrm{1}}^n(x_i-\underset{\mathrm{¯}}{x}{)}^{\mathrm{2}}{\sum }_{i=\mathrm{1}}^n(y_i-\underset{\mathrm{¯}}{y}{)}^{\mathrm{2}}}}},\end{array}$$

with n as the time series length, x_i and y_i as the single time step values within the time series, and $\underset{\mathrm{¯}}{x}$ and $\underset{\mathrm{¯}}{y}$ as the sample means. We calculated values for each tower and then ranked r according to their absolute value to identify the main drivers of NEE.

2.4.4 Information theory

To consider collinear relationships and potential lagging effects between NEE and its controls, we extended our analysis using information theory. Metrics of mutual information (MI) are powerful tools for understanding nonlinear and feedback-driven relationships in complex ecosystems (Chamberlain et al., 2020; Knox et al., 2018). MI is a non-parametric method that can disentangle interactions on different scales (Chamberlain et al., 2018; Knox et al., 2018; Sturtevant et al., 2016), by describing the average tendency for joint states of two variables X and Y to co-occur (Fraser and Swinney, 1986). This means it quantifies the amount of information that two variables X and Y hold in common, or the reduction of uncertainty of one variable, given the knowledge of the other (Chamberlain et al., 2020; Knox et al., 2021). It is a normalized measure of the statistical dependence of Y on X and no prior knowledge about their relationship is needed (Liu et al., 2022). Larger values indicate either higher dependence or a stronger interaction between the variables. With Shannon entropy (H_x) we can quantify the uncertainty in a system:

$$\begin{array}{}\text{(3)}& H_x=-{\sum }_{x_t}p\left(x_t\right){\log }_{\mathrm{2}}p\left(x_t\right),\end{array}$$

with p(x) as the marginal probability distribution of X and X_t as the different states of X in the time series t. Here we discretized the states of continuous variables into 10 fixed-interval histogram bins, as Sturtevant et al. (2016) and Ruddell and Kumar (2009) showed that 10 histogram bins ensure sufficient resolution for a robust estimate. MI was then calculated with both the marginal and joint probability distributions of X and Y, p(x,y):

$$\begin{array}{}\text{(4)}& \mathrm{MI}={\sum }_{x_t,y_t}p\left(x_t,y_t\right){\log }_{\mathrm{2}}{\displaystyle \frac{p\left(x_t,y_t\right)}{p\left(x_t\right)p\left(y_t\right)}}.\end{array}$$

To make the MI between NEE and different potential drivers comparable, we used a normalized form of MI:

$$\begin{array}{}\text{(5)}& {\mathrm{MI}}_{\mathrm{sync}}={\displaystyle \frac{\mathrm{MI}}{H_y}}.\end{array}$$

We refer to this relative MI as synchronous MI (MI_sync), as it depicts the interaction between X and Y at the concurrent time step. An additional strength of MI lies in its capability to account for the temporal direction (τ) of the interaction between X and Y (Liu et al., 2022):

$$\begin{array}{}\text{(6)}& \begin{array}{rl}{\mathrm{MI}}_{\max }& =max\left({\mathrm{MI}}_{{\mathrm{sync}}_{\left(\mathit{\tau }\right)}}\right)\\ & =max\left({\displaystyle \frac{{\sum }_{x_{t-\mathit{\tau }}}{\sum }_{y_t}p\left(x_{t-\mathit{\tau }},y_t\right){\log }_{\mathrm{2}}\frac{p\left(x_{t-\mathit{\tau }},y_t\right)}{p\left(x_{t-\mathit{\tau }}\right),p\left(y_t\right)}}{-{\sum }_{y_t}p\left(y_t\right){\log }_{\mathrm{2}}p\left(y_t\right)}}\right)\end{array}.\end{array}$$

When τ is positive or negative, the interaction between X and Y is characterized as asynchronous, with τ showing the lead or lag in Y relative to X, respectively. We chose 60 d as the maximum value for τ to check if the potential driving variable (Y) is leading NEE (X) or vice versa (Liu et al., 2022). We then picked the highest MI value (MI_max) in this window and the respective day of its occurrence. If MI_sync > MI_max, the interaction is synchronous. If MI_sync < MI_max, the interaction is asynchronous. If τ<0, Y lags X. If τ>0, Y leads X and can therefore be characterized as a driver or control of X. Significance thresholds were calculated from the 95th percentile (p<0.05) of 1000 Monte Carlo random walks of the independent variable (Chamberlain et al., 2020; Ruddell and Kumar, 2009). We calculated MI measures and confidence thresholds in R, based on functions by Chamberlain et al. (2020).

We determined MI_sync and MI_max for the 7-year time series (2016–2022) from the reconstructed seasonal signal. Gap-filled time steps by SSA were removed before the calculation of both r and MI measures. We kept them only for NDVI, as the gap filling is based on the original time series and does not depend on other variables (Walther et al., 2022). Gap filling therefore does not confound the analysis of potential drivers.

2.4.5 Phenological seasons

As the NEE controls vary in their importance in different seasons (Baldocchi and Arias Ortiz, 2024), we calculated MI_sync for each season to better capture how the nutrient addition and stoichiometric balance change the importance of different drivers over the study period. As this ecosystem's strong seasonality is reflected in vegetation activity, we assigned seasons using PhenoCam imagery. We defined phenological seasons following Nair et al. (2024). Phenological transition dates were extracted using GCC at all three sites according to changes between stationary and rising or declining greenness (Luo et al., 2018). Then, phenological transition dates averaged across the three sites for each year were calculated. According to these dates, each day of the 7-year time series was assigned to one season, describing different phases of net vegetation activity (i.e., spring, drydown, summer, autumn, and winter). Figure 2 illustrates a typical annual cycle of the seasons at Majadas de Tiétar.

We calculated MI_sync values for each pair (NEE and potential driving variable) of interest in each season across all 7 years. In addition, we estimated yearly MI_sync for each single season (35 data points) to evaluate how sensitivity of NEE to drivers developed over time. To isolate the fertilization effect on the importance of different drivers for NEE, we calculated the differences in the MI_sync values of each season in each year between the fertilized plots and the control plot, i.e., NT–CT and NPT–CT, referred to as MI_diff. We plotted the MI_diff values for each season along the 7-year period and calculated linear regressions to confirm whether there are significant trends in the importance of drivers. The significance level was set at p<0.05. Variables with MI_max<0.2 were discarded.

Figure 2 A typical annual cycle of midday green chromatic coordinates (GCCs) derived from the grass layer at the control plot in 2018, showing the five phenological seasons – winter, spring, drydown, summer, and autumn. Spring is the main growing season (first peak in May), the grasses become senescent during drydown and dormant in summer, regreening starts in autumn (second peak around November) with the onset of rains, and winter is radiation and temperature limited.

3 Results

3.1 Seasonal NEE variability

At CT, the tower not experiencing any manipulation, the annual ecosystem NEE derived from EC measurements was positive for the 2016–2023 period, with an average annual NEE budget of 90.8 ± 48.0 gC m⁻² yr⁻¹. This indicates that the ecosystem acted as a CO₂ source. With fertilization treatment, the measured ecosystem NEE shifted towards CO₂ neutrality, with annual averages of 34.1 ± 66.7 and 23.1 ± 69.5 gC m⁻² yr⁻¹ at NPT and NT, respectively. Annual NEE budgets fluctuated between positive and negative values at the fertilized plots, while CT consistently showed positive NEE every year. In 2017, 2022, and 2023, we observed high positive NEE values (i.e., stronger CO₂ source) at all three plots. Conversely, in 2016, 2018, and 2021, fertilized areas exhibited higher CO₂ uptake, acting as stronger CO₂ sinks (Fig. 3). This illustrates the high interannual variability of the CO₂ fluxes in this ecosystem and the substantial impact of fertilization. The nutrient addition led to higher seasonal variability of NEE, as shown by the greater yearly standard deviation of the seasonally reconstructed signals. The variability at NT and NPT further exhibited an increasing pattern over time (Fig. 3). In 2017, NEE had comparatively low seasonal variability at all three sites, which might be attributed to the extraordinary dryness in that year.

Figure 3 Annual net ecosystem exchange (NEE) budgets (gC m⁻² yr⁻¹) and yearly standard deviations calculated using the seasonally reconstructed signals at the three sites for the 2016–2023 period. CT – control site, NT – nitrogen-fertilized site, and NPT – nitrogen and phosphorus-fertilized site.

3.2 Key NEE controls

We identified key controls of NEE at the three plots, comparing the results of two different statistical methods: Pearson correlation coefficient r considers only linear relationships between variables; mutual information (MI) accounts for collinear relationships. MI_sync and r values show synchronous relationships. MI_max values can account for leading and lagging interactions by identifying the day of the highest interaction between the potential driver and NEE within a 60 d window.

For all plots, grass-layer GCC and NDVI (at ecosystem level) were the most important predictors of NEE (Fig. 4). Both r and MI identified these proxies representing vegetation greenness as the most significant drivers. They were followed by EF (i.e., the fraction of heat transport that is done by LE), which is influenced by meteorological variables (such as soil moisture, net radiation, and VPD) as well as vegetation properties, like LAI.

At CT, T_soil_bc and Ta15 further exhibited strong interactions with NEE using both r and MI_sync (Fig. 4a, b). Variables describing water availability (VPD, SWCn, and Rh) were ranked in the middle ranges by MI_sync. The MI analysis provided deeper insights into the interactions between the environment and NEE by considering leading and lagging effects, as shown by MI_max (Fig. 4c). NDVI showed the highest interaction with NEE at a time lag of 16 d, and gcc_gr had a lag of 7 d. When considering leading and lagging effects, EF became relatively less important. Soil temperatures were identified amongst the five most important controls. SWCn was also important with a 20 d lag. Air temperature and VPD showed the highest interactions with lags of around a month. Radiation-related variables like PAR and SWDR exhibited long lag times in their highest interaction with NEE (60 d and 53 d, respectively). All MI values can be found in the Supplement (Table S1).

At NT, soil temperatures, VPD, SWCn, and air temperatures were among the most significant controls identified by both synchronous methods, following the primary drivers vegetation greenness and EF. NDVI showed the highest interaction with NEE with a lag of 12 d, followed by gcc_gr with a lag of 6 d. Soil temperatures exhibited the highest interactions with lags of around a month, while air temperatures showed the highest interaction at a lag of 26 d. Moisture-related variables all showed similar time lags (16–20 d). EF had the highest interaction with NEE at a lag of 2 weeks. Shortwave-radiation-related variables showed a strongly lagged effect (i.e., PAR 59 d, SWDR 57 d) (Fig. 4f).

At NPT, both r and MI_sync detected soil temperatures, air temperatures, and VPD as the most important NEE controls behind gcc_gr and NDVI (Fig. 4g, h). NDVI and gcc_gr led NEE interactions with the strongest interactions at lags of 2 weeks and 10 d, respectively, followed by soil temperatures and air temperatures with the highest interactions at lags of around a month (Fig. 4i). EF showed the highest interaction at a lag of 12 d. Other moisture-related variables like VPD, SWCn, and Rh were also detected to be in the middle ranks by MI_max, with time lags of 20–26 d. PAR and SWDR showed the highest interactions with NEE at time lags of around 50 d (Fig. 4i).

MI and r agreed in the detection of the most important drivers, thereby proving that information theory is applicable to our case. Therefore, in the remainder of this paper, we focus on values obtained using MI, as MI is able to detect collinear relationships as well as leading and lagging effects. We discuss variables with MI_max>0.2 in the following sections to concentrate on the information provided by variables with greater explanatory value.

Figure 4 Pearson correlation coefficient (r) (a, d, g), synchronous mutual information (MI_sync) (b, e, h), and maximum mutual information within a 60 d window (MI_max) (c, f, i) between net ecosystem exchange (NEE) and potential drivers over the 7-year period (2016–2022) at the control plot CT (a–c), the nitrogen-fertilized plot NT (d–f), and the nitrogen- and phosphorus-fertilized plot NPT (g–i). The color scale in the MI_max plots indicates the day when MI_max occurs; positive values indicate that the variable leads NEE, negative values vice versa.

3.3 Effect of fertilization on NEE sensitivity to its controls

The relationships between NEE, biogenic, and environmental variables were asynchronous, as indicated by MI_sync < MI_max for all variables. Therefore, we focused on MI_max to describe the differences in NEE sensitivity to various controls across towers.

N fertilization appeared to shorten the reaction time of NEE to changes in NDVI. GCC at the grass level showed higher explanatory value for NEE at NPT and NT (MI_max=0.37) compared to CT (MI_max=0.33). EF showed only slight differences in interaction strengths among the sites (Fig. 5, Table S1). Relative humidity at two heights showed the lowest interaction with NEE at CT (MI_max=0.24), while the fertilized sites had slightly higher values (0.26–0.27). The reaction time of NEE to relative humidity appeared to decrease with fertilization. VPD appeared to have the higher explanatory value for NEE at NT and NPT and slightly less at CT (MI_max=0.27). The interaction between NEE and air temperatures was slightly higher at the fertilized plots compared to the control.

Regarding radiation variables, PAR seemed to have a higher interaction with NEE at NT than at NPT and CT. Similarly, SWDR showed slightly higher interaction with NEE at NT, while at NPT and CT it was equally strong.

In terms of soil variables, soil temperatures exhibited the strongest interaction with NEE. While soil temperatures below canopy (T_soil_bc) were almost the same across sites (MI_max = 0.33), the importance of soil temperatures under open air conditions was lower at CT than at the fertilized plots. SWCn showed the highest explanatory value for NEE at NT (Fig. 5, Table S1). An overview plot with all variables including the ones with MI_max<0.2 is provided in the Supplement (Fig. S3).

Nutrient addition did not show a substantial effect on the sensitivity of ecosystem NEE to different drivers over the 7-year scale when considering the whole time series together. In the next step, we examine the different seasons in greater detail.

Figure 5 (a) Synchronous (MI_sync, grey) and maximum (MI_max, colors) mutual information at the control site (CT), the nitrogen-fertilized site (NT), and the nitrogen- and phosphorus-fertilized site (NPT) at the seasonal scale. The color scale indicates the day when MI_max occurs, with positive values indicating that the variable leads net ecosystem exchange (NEE) and vice versa. (b) MI_sync (dotted lines) and MI_max (solid lines) values at the three sites. Variables with MI_max < 0.2 are not shown here.

3.4 Identifying driver importance in different phenological seasons

We split the 7-year dataset into five different phenological seasons based on the grass-layer GCC and calculated MI_sync between NEE and each of the drivers. This analysis showed that the most important drivers differed between seasons and treatments (Table 2).

In winter, the water vapor transfers of available energy, represented by EF, show a strong interaction with NEE at NPT and CT. Furthermore, NDVI and tree-layer albedo, as well as radiation parameters PAR and SWDR, were important in explaining NEE variations.

In spring (i.e., the main growing season), NDVI and GCC at grass and tree levels showed the strongest interactions with NEE, indicating that NEE was dominated by photosynthetic activity (GPP) during this season. Furthermore, soil temperatures showed strong interactions with NEE at CT and NPT, but not at NT.

During the drydown phase, NEE was dominated by NDVI across treatments, with gcc_gr also showing strong interactions with NEE at CT and NPT. At NT, VPD exhibited a strong link with NEE, which was not as dominant at the other plots.

In summer, soil temperatures showed high interactions with NEE, possibly relating to soil respiration. Additionally, SWDR and PAR were important in explaining NEE variations during this season. At CT and NT, gcc_tr became important, which was logical as the grass layer becomes senescent in the summer and is dormant in terms of ecosystem carbon flux. At NPT, gcc_gr showed a higher interaction with ecosystem NEE than gcc_tr.

In autumn, regreening starts with the onset of rains, and NDVI and gcc_gr showed strong interactions with NEE, as GPP starts to dominate NEE again, driven by photosynthetic activity. Additionally, soil temperatures and air temperatures show a strong link with NEE in autumn. The strongest link between soil temperatures and NEE is found at CT.

Table 2 The five most important drivers in each phenological season at each tower derived using synchronous mutual information. CT: control site, NT: nitrogen-fertilized site, NPT: nitrogen + phosphorus-fertilized site.

3.5 Changes in NEE sensitivity over time

Using yearly MI_sync for each single season, we observed that with N addition, NEE became less sensitive to certain variables during autumn (i.e., the regreening phase), the drydown phase, and winter over time (Fig. 6). Specifically, in autumn the sensitivity of ecosystem NEE to changes in air temperature (Ta15), shortwave radiation (SWDR and PAR), and NDVI decreased significantly over the 7-year period. In the drydown phase, the sensitivity of ecosystem NEE to changes in relative humidity (Rh02 and Rh15) and soil heat flux (SHF_op) also decreased significantly. In winter, however, we observed a significant increase in the sensitivity of NEE to variations in PAR, SWDR, and grass-layer GCC.

With the addition of N + P, significant changes in NEE sensitivity over time were observed in all seasons except the drydown phase (Fig. 6). In autumn, the fertilization with N and P led to a significant decrease in NEE sensitivity to air and soil temperatures (Ta02 and T_soil_op), PAR, and VPD. In spring, which is the main growing season, NEE sensitivity to variations in PAR increased significantly. In summer, NEE became significantly more sensitive to changes in grass-layer GCC. In winter, NEE shows a significant increase in sensitivity to changes in both gcc_gr and SWCn.

4 Discussion

4.1 Nutrient addition increases seasonal NEE variability, driven by grass layer

The annual NEE average from EC measurements is about 91 ± 45 gC m⁻² yr⁻¹, suggesting the unfertilized site acts as a carbon source with a high interannual variability. Similar mean and variability were found in other semi-arid savannas, such as in Kruger National Park in South Africa (75 ± 105 gC m⁻² yr⁻¹) (Archibald et al., 2009). However, semi-arid savannas can also act as carbon sinks. In California, a similar oak savanna (i.e., Tonzi Ranch) was observed to be a carbon sink (values from −144 to −35 gC m⁻² yr⁻¹), while the neighboring grassland (i.e., Vaira Ranch) was found to be a carbon source (−88 to 189 gC m⁻² yr⁻¹) (Ma et al., 2007). In Dakar, Senegal, a Sahelian savanna ecosystem acted as a carbon sink with an average annual NEE budget of −180 ± 29 gC m⁻² yr⁻¹ (Wieckowski et al., 2024). Another natural West African savanna in the south of Burkina Faso has been found to be a strong sink of CO₂ (−864 to −1299 g CO₂ m⁻² yr⁻¹), while two degraded sites nearby were CO₂ sources (118 to 605 g CO₂ m⁻² yr⁻¹) (Berger et al., 2019). Our results indicate that both nutrient addition cases enhance seasonal NEE variability compared to the control. Additionally, seasonal variability increases over time at the fertilized plots. Looking at the difference between annual NEE maximum and annual NEE minimum, we notice substantially increasing trends at both sites, with a significant trend at the NPT plot (Fig. 7). These trends are slightly more driven by increasing maximum values. We argue that this nutrient effect is dominated by grass layer, which substantially controls the NEE dynamics in this system. Our analysis supports the fact that grass-layer GCC and NDVI are most important in explaining NEE variations across treatments (Figs. 4, 5). Both indices represent grass-layer greenness, as the larger fraction of the surface consists of annual grasses (Bogdanovich et al., 2021) and remotely sensed NDVI is dominated by the herbaceous layer.

At the study sites, the added nutrients mostly stay in the herbaceous layer (El-Madany et al., 2021), and it is therefore more affected by the nutrient manipulation than the tree layer. Nutrient addition has been found to increase root biomass and root length density (Nair et al., 2019) and to enable greater absorption for and allocation to leaves. In the leaves, N enhances the photosynthetic capacity (Fleischer et al., 2013), which supports the faster increase in maximum GPP and biomass in the fertilized plots, as confirmed by Luo et al. (2020). NT and NPT show higher productivity and therefore higher biomass compared to CT (Luo et al., 2020). As the grass layer is senescent in summer, the higher biomass results in a greater amount of dead biomass, which will then be respirated by soil microbes (Manzoni et al., 2020; Moyano et al., 2013) as soon as there is sufficient water available (Huxman et al., 2004). This indicates that there is a higher carbon turnover at the fertilized plots, leading to an increased range of NEE within a year (Fig. 7). It agrees with findings from Ma et al. (2016), which, at a Californian oak grass savanna, found that the amount of grass litter determines the size of the fast carbon pool in consecutive seasons.

Figure 6 Linear regressions of the seasonal synchronous mutual information difference (MI_diff) between NT and CT (top four panels) and between NPT and CT (bottom four panels) in different phenological seasons. CT refers to the control site, NT to the nitrogen-fertilized site, and NPT to the nitrogen + phosphorus-fertilized site. Find variable definitions in Table 1. Only the relationships with significant trends are shown. Significance level is set at p<0.05. Variables with overall MI < 0.2 at all towers are not shown.

Evergreen tree species have relatively constant foliage amount throughout the year and are able to use their deeper roots to access lower water resources in the soil (Baldocchi et al., 2004; Rolo and Moreno, 2012). The herbaceous layer is strongly dependent on rainfall variations as it accesses water from the topsoil with a dense near-surface root system (Ward et al., 2013). This layer is therefore much more sensitive to intra- and inter-annual climate variations (Luo et al., 2020). This is probably the reason why the seasonal NEE variability was very low at all sites in 2017. We attribute this to extraordinary dryness in that year, as dryness can lead to severe decreases in both GPP and R_eco in this type of ecosystem (Ma et al., 2007).

Figure 7 Annual range of net ecosystem exchange (NEE) (i.e., maximum NEE minus minimum NEE) in gC m⁻² yr⁻¹ calculated using the seasonally reconstructed signals at the control site (CT), the nitrogen-fertilized site (NT), and the nitrogen + phosphorus-fertilized site (NPT). The range at NPT (p value = 0.049) significantly increased over 8 years. The ranges at NT (p value = 0.116) and CT (p value = 0.270) did not significantly increase.

4.2 Key controls of seasonal NEE

Our results indicate that proxies for vegetation greenness (NDVI and GCC at the grass layer derived from satellite and PhenoCam data, respectively) are the primary factors influencing the seasonal NEE signal in this ecosystem across treatments (Fig. 4). However, depending on the season, other variables such as air temperatures, VPD, moisture-related variables, and soil temperatures can also be important.

Many studies identify NDVI, a proxy for vegetation greenness and photosynthesis, as a primary predictor of NEE (Del Grosso et al., 2018; Hermance et al., 2015; Morgan et al., 2016). NDVI, generally derived from satellite data, represents ecosystem greenness, and its connection with ecosystem CO₂ fluxes has been intensively studied (Barnes et al., 2016; Hermance et al., 2015; Morgan et al., 2016; Running and Nemani, 1988). However, quantifying the importance of coexisting vegetation layers is more complex and less understood. Digital repeat cameras and the vegetation greenness indices derived from them provide a powerful tool for analyzing the greenness of different plant types (Migliavacca et al., 2011; Petach et al., 2014; Richardson et al., 2009; Yan et al., 2019) and their influence on ecosystem fluxes (Luo et al., 2018; Moore et al., 2016, 2017; Wingate et al., 2015). Our analysis confirms that grass-layer dynamics are dominant in controlling seasonal ecosystem NEE at this site. Furthermore, we suggest that alternatively to NDVI, the enhanced vegetation index (EVI) can be considered as a representation of vegetation greenness, as it is found to covary closely with carbon flux in semi-arid ecosystems (Maluleke et al., 2024).

In situ measurements of vegetation greenness, however, are not available at all EC sites. We found that EF (i.e., evaporative fraction), representing the fraction of available energy transported by LE, is the third most important driver across treatments and methods on the synchronous scale. EF is strongly influenced by net radiation and water-related variables like soil moisture and VPD (Gentine et al., 2007; Tong et al., 2022). Since it depends on the portion of LE that is transpired by plants, it is also impacted by LAI (Gentine et al., 2007). EF therefore serves as a bridge between meteorological and vegetation controls. We suggest that at semi-arid sites where GCC measurements are not available, EF, calculated from measured LE and H, can serve as an important predictor of NEE.

In water-limited semi-arid ecosystems, NEE variations are typically dominated by soil-water-related variables such as SWCn and precipitation (Archibald et al., 2009; Baldocchi and Arias Ortiz, 2024; Huang et al., 2016b; Morgan et al., 2016). These variables usually exert a greater influence than radiation and temperature (Del Grosso et al., 2018; Kannenberg et al., 2024). Water availability promotes plant photosynthesis (Parton et al., 2012), but rain pulses can also enhance heterotrophic respiration rates (Morgan et al., 2016). While we do not use precipitation data for the MI analysis, as the measurements tend to be zero on many days and cannot be used in MI (Gong et al., 2014), SWCn can capture topsoil moisture and indicate precipitation pulses. Additionally, EF can serve as a proxy for these pulses. In our analysis we identify EF as one of the most important NEE drivers, while other moisture-related variables (e.g., SWCn, VPD, Rh) are generally ranked lower in importance compared to air and soil temperatures (Fig. 4). This might point to the relationship of EF with LAI being the dominant one in this context, as the vegetation indices showed greater importance than other water-related variables.

Air temperature can directly affect the speed of the enzyme responsible for carbon fixation and the rate of photosynthetic electron transport (Leuning, 2002; Xu and Baldocchi, 2003). Additionally, air temperature impacts the availability of photosynthetic enzymes, membrane fluidity, and the expression of associated proteins (Yamori et al., 2014). However, our results show that soil temperatures, both under oak trees and in open areas, play a significant role in explaining seasonal ecosystem NEE variations (Fig. 4), exceeding the importance of air temperatures. Soil respiration, one of the components in R_eco, is highly sensitive to soil temperature (Conant et al., 2000), and elevated soil temperatures are associated with increased soil respiration in semi-arid ecosystems (Richardson et al., 2012). Soil temperatures influence heterotrophic respiration, which constitutes a substantial part of ecosystem NEE at our site (Casals et al., 2011). The high importance of soil temperatures suggests that R_eco dominates ecosystem processes. This is especially relevant as the trend of increasing soil temperatures is stronger than that of air temperatures in the Mediterranean, particularly in grasslands with low soil moisture availability (Wang et al., 2024).

Radiation parameters, in particular PAR, do not appear to play a crucial role at the seasonal scale. While other studies have identified it as a major control of NEE in semi-arid ecosystems (Baldocchi and Arias Ortiz, 2024), we argue that PAR predominantly influences the daily NEE signal (Fig. S2), but its importance diminishes on seasonal timescales.

Overall, we observe only marginal differences between the treatments when considering the 7-year period (2016–2022) together. The added nutrients, particularly N, are primarily absorbed by the herbaceous layer (El-Madany et al., 2021) that senesces annually. Consequently, some of the added nutrients may be lost from the system, diminishing the long-term effect of the fertilization. By calculating MI_sync and MI_max for 1 year post-fertilization (March 2016–February 2017), we observe greater differences between the three plots in MI values and lag times (Fig. S4). Additionally, the ecosystem is strongly water-limited in the summer and energy-limited in the winter (Luo et al., 2018; Nair et al., 2019). These limitations can be more pronounced than nutrient limitations in the respective seasons, overshadowing the effects of added nutrients when analyzing the entire dataset together. Therefore, we divided the dataset into five phenological seasons to gain deeper insights into how added nutrients and altered stoichiometric balance affect seasonal NEE.

4.3 Fertilization effects in different phenological seasons

Looking into phenological seasons gives a deeper insight into how environmental variables influence seasonal NEE and how N : P levels affect this relationship. We find that nutrient addition has an effect on NEE–control relationships when other limitations are not too strong.

In winter, the ecosystem is energy-limited (Luo et al., 2018); therefore radiation components (i.e., PAR and SWDR) are important predictors for NEE. Tree albedo shows strong interactions with NEE at CT and NT, and NDVI shows strong interactions with NEE at both fertilized plots. Plant growth is enhanced by added nutrients (Luo et al., 2020) and flourishes due to abundant water availability (Lee et al., 2010) in this season. Also, N + P addition can lead to increased species diversity due to diminished nutrient limitation, facilitating the coexistence of multiple species (Köbel et al., 2024). Additionally, EF shares high mutual information with NEE variations. This is likely because respiration does not change significantly during this period and VPD is relatively low, leading to a strong coupling between NEE and LE. The stomatal conductance of tree transpiration is not too strong in winter, as soil water is abundant (Klein et al., 2013).

In the primary growing season, spring, NEE is typically dominated by GPP. The key drivers during this season across sites are the NDVI and the GCC of both the herbaceous and tree layers (Table 2). Water is usually abundant promoting plant photosynthesis during moderate temperatures in this time (Baldocchi and Arias Ortiz, 2024). These conditions are further supported by increased day length and higher radiation levels (Luo et al., 2018). The rise in incoming radiation, extended daylight hours, and elevated temperatures, coupled with the increased atmospheric evaporative demand (i.e., higher VPD), lead to a strong correlation between precipitation and both GCC and GPP, as observed in various Mediterranean ecosystems (Diodato and Bellocchi, 2008; Luo et al., 2018; Ma et al., 2007).

In the water-limited seasons, the nutrient effect is minimal as the grass layer is dormant and nutrients are not available due to lack of water. During the drydown period, soil moisture (i.e., SWCn) decreases drastically due to increasing air temperatures and scarce rainfall (Battista et al., 2018; Luo et al., 2018). This induces annual grasses to become senescent, leading to a loss of chlorophyll content (Luo et al., 2018). The rate of this senescence can determine whether NEE becomes positive or negative during this time. NDVI and grass-layer GCC, the most important predictors of NEE in this season across sites, can provide insights into the drydown rate. At NT grass-layer GCC is less important, which we attribute to a more rapid drydown, causing the grass layer to enter dormancy earlier than at other sites (Luo et al., 2020). This is because N addition promotes faster water usage (Luo et al., 2020), accelerating the decrease in SWCn and thereby hampering photosynthesis. It leads to a higher transpiration at NT compared to the other sites, potentially due to rhizosphere priming to increase P mobilization through microbes, as adding only N to the system leads to a P deficiency (El-Madany et al., 2021). In addition, N fertilization can alter species diversity and composition, likely selecting for species that senesce early (Wang and Tang, 2019). The higher interaction of soil temperatures with NEE in this season compared to the wetter seasons shows that R_eco starts dominating NEE, as R_eco is strongly connected to soil temperatures (Metz et al., 2023). VPD is a stronger control of NEE at NT compared to the other two plots. Transpiration is highest at NT, as plants transpire more to obtain limited P from the soil (El-Madany et al., 2021; Pang et al., 2018; Rose et al., 2018). It is therefore more sensitive to changes in VPD.

In summer, the driest period in the ecosystem, R_eco dominates NEE and thus we find a strong interaction between NEE, soil temperature, and soil heat flux (i.e., SHF_op and SHF_bc). PAR is important for predicting seasonal NEE, showing the strongest interaction at CT. The importance of PAR is lower at NT and lowest at NPT. N + P addition increases the light-use efficiency most because P has a positive effect on photochemical quenching in leaves and on active fluorescence measurements (Martini et al., 2019; Singh and Reddy, 2014), leading to less dependence of NEE on radiation parameters at that site. At CT and NT, tree-layer GCC is important as the grass layer becomes senescent in the summer and is dormant in terms of ecosystem carbon flux. Since the greenness of the oak trees is relatively constant throughout the year, GPP is mainly determined by the tree layer in the summer months (Luo et al., 2018). However, gcc_gr shows a higher interaction with NEE than gcc_tr at NPT. Even though the majority of the grass layer is mostly dead in this season, there are some perennial species (e.g., Cynodon dactylon) remaining green for longer in summer that can regreen after rain events (personal communication with local collaborators). Therefore, N + P addition very likely leads to a consequential change in species composition (Köbel et al., 2024) with an increase in perennial species or an increase in their productivity. So far, it has been found that N + P addition can lead to an increasing number of forbs (Köbel et al., 2024), which tend to senesce later than other herbaceous species at the site (Luo et al., 2020). Nevertheless, the occurrence of summer-green species following nutrient addition will have to be investigated further.

During the regreening of the herbaceous layer starting in autumn, NDVI shows the strongest interaction with NEE at the fertilized plots – but not at the control plot. This aligns with previous studies showing that the green-up happens faster and the maximum GPP is higher at the fertilized plots, resulting from larger resource utilization at NT and/or improved resource-use efficiency at NPT (Luo et al., 2020). With the increase in soil moisture in early autumn, a greater quantity of organic and inorganic nutrients becomes available to plants (Agehara and Warncke, 2005; Luo et al., 2020). N availability in the soil is expected to be highest at this time (Morris et al., 2019), leading to higher net carbon uptake rates (El-Madany et al., 2021). Leaves quickly expand and pigments rapidly increase during this green-up period (Croft et al., 2015). At CT, the green-up happens later than at the fertilized plots, and NEE is dominated for a longer time by R_eco rather than photosynthetic activity (Luo et al., 2020). Our results indicate that soil temperatures below oak trees are more important than those in open areas during this season (Table 2). The carbon pools under oak trees are the largest, providing substantial material for heterotrophic decomposition (Casals et al., 2009). During autumn, after a prolonged dry season in which a significant amount of litter and organic material has already been decomposed by microbes, litter remains available for further heterotrophic decomposition mainly below the trees. This ongoing decomposition under oak trees contributes to R_eco, especially as the onset of rains enhances microbial activity due to increased water availability (Borken and Matzner, 2009). Additionally, the topsoil remains wet for longer after rain pulses under oak trees compared to open areas, as soil moisture during this season is primarily influenced by evaporation from the relatively bare soil. Therefore, differences in soil respiration between open and shaded pastures can also be attributed to variations in soil moisture.

The analysis of driver importance in different phenological seasons provides significant insights into ecosystem processes. However, some variables must be interpreted with caution. The soil properties at this site are highly heterogeneous, which affects the representativity of variables like soil temperature, soil water content, and soil heat flux in the EC flux footprint (Luo et al., 2018; Paulus et al., 2022). This is particularly relevant given the substantial differences between below-canopy and open-air soil conditions. To address this, we have separated the measurements into areas under the oak tree canopy and sunlit areas (T_soil_bc and SHF_bc; T_soil_op and SHF_op). Despite this effort, the local soil heterogeneity remains complex, influenced by varying proportions of sand, clay, and soil organic carbon (Casals et al., 2011; Weiner et al., 2018). Therefore, it is important to consider that these measures may not fully capture the sensitivity differences in the ecosystem.

4.4 Future implications

In winter, the ecosystem has abundant water availability, and energy becomes the primary limiting factor after nutrients were added. With the N-only addition, we observe that NEE becomes significantly more sensitive to changes in the radiation components, PAR and SWDR (Fig. 6). However, the addition of N + P results in a significant increase in sensitivity to changes in soil water content rather than to radiation components. N + P addition enhances water-use efficiency in the ecosystem (El-Madany et al., 2021; Martini et al., 2019), and consequently, water can be used more efficiently for photosynthesis – with similarly low radiation levels – and increased water availability could lead to higher GPP. N fertilization primarily affects the herbaceous layer (El-Madany et al., 2021), and our results agree with this, showing a significantly increased sensitivity of NEE to grass-layer greenness in winter at NT and an even steeper increase at NPT (Fig. 6). At the N + P plot there are more nutrients available at a higher N : P stoichiometric balance.

In spring, the sensitivity to tree-layer greenness decreases with N fertilization. An experimental study by Biro et al. (2024) supports these findings, demonstrating that N addition results in decreased tree growth due to competition with grass, which also intensively forages for P. The study suggests that grasses likely prevail in belowground competition, primarily due to their substantial root biomass allocation and investment in nutrient-acquiring enzymes, such as phosphatase. These adaptations enable grasses to efficiently sequester both N and P from the soil, thereby outcompeting trees for these essential nutrients (Biro et al., 2024; Rolo and Moreno, 2012). With the addition of N + P, we observe that the NEE sensitivity to PAR increases significantly in spring (Fig. 6). Water and nutrients are abundant in this season at NPT, making the availability of energy more crucial.

In water-limited seasons, ecosystem processes behave quite differently, and we observe less effect of nutrient addition. With N + P addition, there is no significant trend in NEE sensitivity to its drivers, except for a significantly increased sensitivity to grass-layer greenness in summer. This agrees with our previous findings that in summer gcc_gr is amongst the most important drivers at NPT (Table 2). This reflects changes in the species composition with N + P fertilization, particularly enhancing the growth and diversity of forbs and perennial species. We argue that long-term N + P addition could even lead to an increased productivity, resulting in an increasing importance of grass-layer greenness for ecosystem NEE.

The significant increase of the yearly NEE range at NPT over time (Fig. 7) is very likely caused by the increased NEE sensitivity to drivers in spring and summer, as the minimum NEE (usually occurring in spring) becomes more negative and maximum NEE (usually occurring in summer) becomes more positive. Consequently, the increased NEE sensitivity to changes in PAR in spring and increased sensitivity to gcc_gr in summer might enhance the size of this annual range.

In the drydown phase, we observe that with N addition, the sensitivity of ecosystem NEE to changes in relative humidity (i.e., Rh02 and Rh15) and SHF_op decreases significantly. This indicates that the ecosystem might become more resistant to variations in these variables in the future.

In autumn, both fertilized sites become less sensitive to changes in atmospheric variables such as the radiation components PAR and SWDR, air temperatures, and VPD, compared to the control plot. This indicates that water availability is predominantly important for NEE with added nutrients, while sensitivity to other variables decreases. It is possible that either the vegetation or the microbes become less restricted by these variables.

We conclude that with more N input from human activities entering terrestrial ecosystems (Peñuelas et al., 2013), savannas may become less sensitive to environmental factors like humidity, radiation, and temperature during the transitional seasons (i.e., drydown and regreening). These seasons determine the start and end of an active grass layer and therefore dominate the annual carbon balance of the ecosystem. In addition, we expect the NEE variability to increase even more in the future with more N deposition and a changing climate. To note, the results of this study cannot explain how long-term nutrient addition affects the ecosystem's resistance to extreme events. As the timing of the application of N and P was chosen to increase the possibility that the sites could be used by vegetation in the next growing season, the observed changes in NEE and driver importance at NT and NPT might have been smaller if fertilization had been applied at different timing.

5 Conclusions

We analyzed a long-term (2016–2022/23) dataset of flux, biometeorological, satellite, and PhenoCam data from the semi-arid experimental site, Majadas de Tiétar, to evaluate the importance of different drivers for NEE across different nutrient levels and balances. To detect the most important drivers, we used only daytime daily values of observed data in order to extract the seasonal signals of all variables using singular spectrum analysis.

With both Pearson correlation and mutual information analysis, we show that the grass-layer greenness drives seasonal variations in NEE across all treatments and that the addition of both N and N + P increases the seasonal NEE variability. We find that soil temperatures are more important in explaining NEE variations than previously expected. When looking into the 7-year data, the water–energy limitation cycles overshadow the nutrient addition effect. Dividing the dataset into phenological seasons reveals how environmental variables and nutrient manipulation influenced NEE on a seasonal scale. Altered nutrient levels affect NEE–control relationships when water and energy limitations are not too strong, particularly during the primary growing season in spring, where NDVI and grass-layer GCC are key drivers. In autumn, NDVI shows the strongest interaction with NEE at fertilized plots, indicating faster green-up and higher GPP due to enhanced nutrient availability. During drier seasons, nutrient effects are less pronounced as the grass layer becomes dormant.

N and N + P additions significantly alter the sensitivity of NEE to environmental controls over time. In winter, N addition increases NEE sensitivity to radiation, while N + P addition increases its sensitivity to changes in soil water content. In spring, N + P addition increases sensitivity to PAR. The herbaceous layer primarily benefits from nutrient additions, leading to increased sensitivity of NEE to grass-layer greenness and decreased sensitivity to tree-layer greenness. During water-limited seasons, nutrient effects were minimal, except for increased importance of grass-layer GCC in summer at NPT, indicating an increase in abundance and/or productivity with N + P treatment due to changed species composition and higher biodiversity. We conclude that with increasing anthropogenic N deposition, the carbon dynamics of savannas might become even more variable in the future and might also become more resistant to variations in some atmospheric variables in the transitional seasons, which are important periods in determining the annual carbon balance of the ecosystem. However, their responses to extreme events in the future remain to be explored.