Exploring Ecosystem Functioning in Spain with Gross and Net Primary Production Time Series

: The main objective of this study is to analyze the spatial and temporal variability of gross and net primary production (GPP and NPP) in Peninsular Spain across 15 years (2004–2018) and determine the relationship of those carbon ﬂuxes with precipitation and air temperature. A time series study of daily GPP, NPP, mean air temperature, and monthly standardized precipitation index (SPI) at 1 km spatial resolution is conducted to analyze the ecosystem status and adaptation to changing environmental conditions. Spatial variability is analyzed for vegetation and speciﬁc forest types. Temporal dynamics are examined from a multiresolution analysis based on the wavelet transform (MRA-WT). The Mann–Kendall nonparametric test and the Theil–Sen slope are applied to quantify the magnitude and direction of trends (increasing or decreasing) within the time series. The use of MRA-WT to extract the annual component from daily series increased the number of statistically signiﬁcant pixels. At pixel level, larger signiﬁcant GPP and NPP negative changes ( p -value < 0.1) are observed, especially in southeastern Spain, eastern Mediterranean coastland, and central Spain. At annual temporal scale, forests and irrigated crops are estimated to have twice the GPP of rainfed crops, shrublands, grasslands, and sparse vegetation. Within forest types, deciduous broadleaved trees exhibited the greatest annual NPP, followed by evergreen broadleaved and evergreen needle-leaved tree species. Carbon ﬂuxes trends were correlated with precipitation. The temporal analysis based on daily TS demonstrated an increase of accuracy in the trend estimates since more signiﬁcant pixels were obtained as compared to annual resolution studies (72% as to only 17%).


Introduction
A full description of the carbon cycle requires detailed information on spatiotemporal patterns of surface-atmosphere carbon fluxes. One of the main carbon fluxes characterizing terrestrial ecosystems and biodiversity is the gross primary production (GPP), i.e., the amount of carbon absorbed by vegetation to perform photosynthesis, since it establishes the main carbon and energy inputs to ecosystems for producing food, wood, and fiber [1]. However, approximately half of the GPP is respired by plants to provide the energy that supports their growth and maintenance. Net primary production (NPP) is the net carbon gain by vegetation and equals the difference between GPP and plant autotrophic respiration [2]. NPP is recognized as the most relevant variable to characterize the state of ecosystems since it represents the fundamental source of energy for all organisms in an ecosystem and is a driver of the most essential of ecosystem services [2].
The role that Earth observation (EO) data play in quantifying carbon fluxes evolved substantially in the last two decades with the beginning of the operational Moderate Resolution Imaging Spectroradiometer (MODIS) GPP product from the National Aeronautics and Space Administration (NASA) Earth Observing System (EOS) program [3]. Since

Material
Four datasets of 1 km spatial resolution and their interrelationship were analyzed: daily GPP, NPP, and mean air temperature (T), along with monthly standardized precipitation index (SPI), all during the 2004-2018 period. This period was chosen due to the availability of some of the data: solar irradiance was available only from 2004, while air temperature and precipitation only until 2018.
SPI and T were obtained from the spatialization of ground measurements. GPP and NPP were simulated respectively by means of a Monteith-like PEM [7] and the Biome-BGC model, both previously optimized for the study area. In the following subsections, the data needed to obtain these datasets, as well as the data processing, are briefly described.

SPI and T
Daily precipitation and minimum and maximum air temperature measurements from ground stations were provided by the Spanish Meteorological Agency (http://www. aemet.es, accessed on 13 January 2022) for the 2004-2018 period. Daily images at 1 km spatial resolution for the same period were obtained from the measurements by ordinary kriging [23]. T was calculated as the mean value between minimum and maximum air temperature. SPI images were computed from daily precipitation images. The SPI quantifies precipitation anomalies by transforming observed precipitation into a gamma distribution for a specific time period [24]. In this study, the 12-month SPI over the period 2004-2018 was computed in order to enhance the identification and understanding of vegetation changes.

Vegetation and Forest Type Maps
The 14-class and 1 km spatial resolution land cover map generated by [25] was used to obtain a vegetation type map, which was used to simulate GPP (Section 2.3). Elevation extracted from the Shuttle Radar Topography Mission 3 arc second digital elevation model [26] and land cover data from the Spanish Land Use Information System (SIOSE) [27] were combined to obtain a forest type map [28], which was used to simulate forest respiration (Section 2.4).

GPP
The PEM used to estimate the daily GPP (g m -2 d -1 ) images is where PAR (MJ m -2 d -1 ) accounts for the photosynthetically active radiation in the 0.4-0.7 µm spectral range. It was computed as 46% of the daily solar irradiance, which was given by the MSG DIDSSF product (LSA-203) [29] reprocessed and available from 2004 up to the present at https://landsaf.ipma.pt/ (accessed on 13 January 2022). This product was very satisfactorily validated in the study area and is used as the input of the MSG/SEVIRI GPP product LSA-411 [30]. The daily f APAR (fraction of absorbed photosynthetically active radiation) at 1 km spatial resolution was obtained by applying the algorithm proposed in [31] to MODIS products MCD43A1 and MCD43A2. This algorithm is actually used to derive the SE-VIRI/MSG f APAR product (LSA-425) [32]. In this study, the f APAR TS were filtered and reconstructed using an optimized locally weighted regression and smoothing scatterplots (LOESS) method [33].
ε W is the water stress coefficient C WS proposed by [34] and applied satisfactorily over Mediterranean regions [12,35,36]. It is calculated as a linear relationship with the water balance between actual and potential evapotranspiration (AET and PET, respectively): If AET is not available and precipitation is lower than PET then it can be assumed that AET is equal to precipitation. In this study, the accumulated precipitation and PET are used because it demonstrated a better performance in the study area [36]. The accumulation period depends on the vegetation type [28]. PET (mm·m -2 ·d -1 ) is calculated, according to [37] parameterization, from the daily global irradiation R g (kJ m -2 d -1 ), as given by the LSA-203 product, and the daily T (in • C): ε max (g MJ -1 ) is the ecosystem light use efficiency under optimal conditions and depends on the vegetation type. The vegetation type map described in Section 2.2 was used to assign the ε max value to each pixel according to the literature [11]. Finally, ε T is the Tmin_scalar used in MOD17 product [4], a linear ramp function of minimum air temperature that reduces ε max accounting for low temperatures, depending on the vegetation type.
This methodology was previously tested in the study area against eddy covariance observations obtaining coefficients of determination between 0.5 and 0.96 in [12], where more details on the estimation of GPP from remotely sensed data can be found.

NPP
Forest NPP was obtained based on the combination of the GPP calculated through equation (1), and the growth (R G ) and maintenance (R M ) respirations simulated by version 4.2 of the ecosystem process model Biome-BGC [38,39]. Daily precipitation, maximum and minimum air temperature images, and the DIDSSF product described above were used. Daylength, daylight average partial pressure of water vapor, and daily average air temperature were simulated by the microclimate simulation model MT-CLIM [40] version 4.3. The needed ecophysiological parameters were obtained in a previous study over an area of similar eco-climatic features [41]. Elevation and forest type described in Section 2.2 were also used. The rooting depth calibration, soil texture maps description, and more details on how to use Biome-BGC over the study area can be found in [28].
In this study, forest NPP is calculated as where f is the actual to potential tree cover ratio and v is the actual to potential growing stock volume (GSV) ratio. GPP calculated from Equation (1) needs to be corrected to account only for the fraction of forest present in the pixel. This expression was proposed in [42], based on the ecosystems distance to climax concept. These ratios are calculated as and v = 50 GSV ρ C/σ (6) where LAI is the annual maximum leaf area index simulated by Biome-BGC, ρ is the basic wood density, C is a biomass expansion factor, σ is the annual maximum dead stem carbon simulated by Biome-BGC, and the factor 50 accounts for the unit conversion to m 3 /ha and the transformation from carbon mass to dry mass (2 kg/kg). The growing stock volume GSV was obtained in a previous study combining remotely sensed and forest inventory data [43]. NPP was only calculated in forest pixels because GSV is only available for trees in the Spanish Forest Inventory described in [44]. For the main Mediterranean forests, ρ and C can be found in [45]. This methodology was previously assessed in the study area using forest inventory observations. Coefficients of determination between 0.5 and 0.7 in [46] were obtained. A detailed description on modeling of forest carbon storage can be found in [46].

Study Area
Peninsular Spain is located in the southwest of Europe, approximately between (44 • , −10 • ) and (36 • , 3 • ) lat/lon. The ecosystem's diversity, mainly caused by its orography (its elevation ranges from sea level to 3479 m), is remarkable. Its geographical location leads to a range of climates from Atlantic in the northwest to semiarid in the southeast. Its annual precipitation follows the same gradient, ranging from around 2000 mm in the northwest to around 200 mm in the southeast. Figure 1a,b show the spatial distribution of the average annual precipitation and the average annual T for the 2004-2018 period.

Study Area
Peninsular Spain is located in the southwest of Europe, approximately between (44°, −10°) and (36°, 3°) lat/lon. The ecosystem's diversity, mainly caused by its orography (its elevation ranges from sea level to 3479 m), is remarkable. Its geographical location leads to a range of climates from Atlantic in the northwest to semiarid in the southeast. Its annual precipitation follows the same gradient, ranging from around 2000 mm in the northwest to around 200 mm in the southeast. Figure 1a,b show the spatial distribution of the average annual precipitation and the average annual T for the 2004-2018 period. According to the land cover map by [25] (Figure 1c), 52% of the pixels of the study area are classified as croplands, 21% as forests, and 25% as shrublands, grasslands, or sparse vegetation. Figure 1d shows the forest type map used in the simulation of forest respirations [28]. Forests in Peninsular Spain present great spatial complexity and GSV irregularity. According to the Third Spanish Forest Inventory, 56% of the forest areas (d) Specific forest type map used to derive the NPP: EBF (evergreen broadleaved forest), LDBF (lowaltitude deciduous broadleaved forest), HDBF (high-altitude deciduous broadleaved forest), LENF (low-altitude evergreen needle-leaved forest), and HENF (high-altitude evergreen needle-leaved forest). These data are described in Section 2.
According to the land cover map by [25] (Figure 1c), 52% of the pixels of the study area are classified as croplands, 21% as forests, and 25% as shrublands, grasslands, or sparse vegetation. Figure 1d shows the forest type map used in the simulation of forest respirations [28]. Forests in Peninsular Spain present great spatial complexity and GSV irregularity. According to the Third Spanish Forest Inventory, 56% of the forest areas covered by tree species are broadleaved forests with Quercus ilex as the most common species, while the other 34% are needle-leaved forests with majority of Pinus halepensis [44].

Site
Vegetation Type Site A is a sparse vegetation region located in the southeast of Spain, at the north of Desierto de Tabernas, in the province of Almería, characterized by pronounced dry conditions causing desertification. Site B is a grassland located at the east of the province of Badajoz, in the interior of Spain. Site C is a shrubland located at the southwest of the province of València, in the east of Spain. Site D is a rainfed cropland located in the province of Cuenca. Site E is an irrigated cropland located in the province of Lleida, in the Ebro basin. Forest sites (F-J) are located in or near to national parks. Site F is an evergreen broadleaved forest, particularly, the most extensive cork forest in Spain and one of the largest in the world dominated by Quercus suber. Site G is a deciduous broadleaved forest located in the second largest and best preserved beech forest in Europe, the national park Selva de Irati, in the province of Navarra (north of Spain). Site H is a high-altitude deciduous broadleaved forest located at the southeast of Picos de Europa. Site I is a low-altitude evergreen needle-leaved forest located in the national park of Sierras de Cazorla, Segura y las Villas, the Spanish largest continuous area of pine forest with representatives of nearly all pine species found in the Iberian Peninsula, the most abundant being Pinus nigra. Site J is a high-altitude evergreen needle-leaved forest located in the north of the national park Serranía de Cuenca, in the province of Guadalajara. All the described sites were selected according to the vegetation type provided by the land cover map produced by [25], which uses land cover data from 1999 to 2006.

Methods
Two procedures are used for the temporal analysis at pixel level. First, the temporal dynamics of GPP, NPP, SPI, and T are assessed using the MRA-WT procedure (Section 4.1). Second, the trend (with magnitude and direction) of the annual component is obtained through the combination of the Mann-Kendall test and the Theil-Sen slope estimator (Section 4.2).

Interannual Component from MRA-WT
The MRA-WT allows for a fast implementation of the discrete wavelet transform by means of a decomposition of the TS into timescales based on powers of two, 2 j (j = 1, . . . , m), where j refers to the different levels of decomposition and m to the highest level considered [47]. As a result of the MRA, the original signal can be reconstructed as where A is the approximation component and D the detail component. The MRA can be understood as a sequence of the A component at the lowest level and complementary information given by the details D time series, which provide information of high-frequency contributions at each temporal resolution. The temporal resolution of each level depends on the center frequency of the selected wavelet and the temporal resolution of the TS. In our case, the Meyer wavelet (with central frequency v c = 0.67213 Hz) is chosen due to its potential for vegetation dynamic analysis shown in previous studies [14,48]. For example, in the case of daily temporal resolution (GPP, NPP, and T in the current study), the approximation component of each level accounts for all the vegetation variations at the following specific time resolutions: A 1 (1.5 days), A 2 (3 days), A 3 (6 days), A 4 (12 days), A 5 (24 days), A 6 (48 days), A 7 (96 days), A 8 (192 days), A 9 (385 days). The detail components provide information about the portion of the signal that can be attributed to variations between the resolutions [j − 1, j]. The sum of detail components D between different levels provides the variability (V) of the signal attributed to the temporal frequencies between these particular levels. Details of this methodology can be found in [14].
The MRA-WT is applied until level nine for daily TS (GPP, NPP, and T) and until level four for monthly SPI series since it gives us the approximation component (A 9 and A 4 , respectively) to study the interannual changes. Figure 2 shows an example of the MRA-WT decomposition applied to the GPP series for a deciduous broadleaved forest pixel. The top panel shows the original TS, whereas the center plot refers to the total variation removed from the interannual part corresponding to days between 6 and 385 days (V 4-9 ). The plot on the bottom shows the interannual (A 9 ) component of the original TS.

Trend Analysis
Once the interannual components are derived, the Mann-Kendall nonparametric t [49,50] and the Theil-Sen slope estimator [51] are applied to statistically detect trends the GPP (QGPP), NPP (QNPP), SPI (QSPI), and T (QT) interannual TS. The Mann-Kendall t has been extensively used to identify the direction and statistical significance of mon tonic trend within each pixel of the image TS [52]. A nonsignificant Mann-Kendall t value (0 value) means the null hypothesis (H0) of this test (no trend exists) may not

Trend Analysis
Once the interannual components are derived, the Mann-Kendall nonparametric test [49,50] and the Theil-Sen slope estimator [51] are applied to statistically detect trends in the GPP (Q GPP ), NPP (Q NPP ), SPI (Q SPI ), and T (Q T ) interannual TS. The Mann-Kendall test has been extensively used to identify the direction and statistical significance of monotonic trend within each pixel of the image TS [52]. A nonsignificant Mann-Kendall test value (0 value) means the null hypothesis (H 0 ) of this test (no trend exists) may not be rejected, whereas a significant test value (1 value) assumes H 0 is rejected, and a trend is considered. H 0 was rejected at p-value under 0.1 significance levels for the daily time series, such as carbon fluxes and temperature TS (GPP, NPP, and T). For SPI, this condition was relaxed to p-values under 0.3 in order to have a larger number of pixels to be compared with the GPP and NPP trend values.
Q GPP , Q NPP , and Q T are, respectively, calculated as the Theil-Sen slope obtained for GPP, NPP, and T divided by the mean value of GPP, NPP, and T, and multiplied by 100. This normalization is not needed with Q SPI since SPI is already normalized by definition [24].
The MRA-WT procedure (Section 4.1) along with the trend analysis (Section 4.2) are performed at pixel level in order to quantify the long-term changes in Peninsular Spain. A local analysis is also carried out for the ten selected sites (see Table 1 in Section 3). The interrelation between carbon fluxes and meteorological variables (air temperature and precipitation) trends are investigated at this scale by calculating correlations between their annual components (A 9 for GPP, NPP, and T, and A 4 for SPI). Figure 3 shows the spatial distribution of average annual GPP and NPP for the study period (2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018). Both GPP and NPP show a strong latitudinal gradient increasing from low to high latitudes. Exceptionally, the most southern region, where site F is located, exhibits average annual GPP and NPP similar to those in the northern regions. Table 2 shows average annual GPP and NPP for the vegetation and forest types indicated in Figure 1c,d. The results reveal that the greatest GPP is found for broadleaved forests (1.6, 1.8, and 2.0 kg m -2 y -1 for EBF, HDBF, and LDBF, respectively), while the smallest GPP is found for sparse vegetation and rainfed crops, both presenting GPP = 0.6 kg m -2 y -1 . The greatest forest NPP is also found for broadleaved forests (0.4, 0.5, and 0.8 kg m -2 y -1 for EBF, HDBF, and LDBF, respectively). Needle-leaved forests reach NPP values up to half of the broadleaved forests NPP (0.2 kg m -2 y -1 ).  ure 1c,d. The results reveal that the greatest GPP is found for broadleaved forests (1.6, 1.8, and 2.0 kg m -2 y -1 for EBF, HDBF, and LDBF, respectively), while the smallest GPP is found for sparse vegetation and rainfed crops, both presenting GPP = 0.6 kg m -2 y -1 . The greatest forest NPP is also found for broadleaved forests (0.4, 0.5, and 0.8 kg m -2 y -1 for EBF, HDBF, and LDBF, respectively). Needle-leaved forests reach NPP values up to half of the broadleaved forests NPP (0.2 kg m -2 y -1 ).   Figure 4 shows the annual GPP and NPP for the selected sites (Table 1) in order to illustrate general behaviors between land covers. A significant gradient is observed in the GPP values, where the highest GPP values are found in broadleaved forests, the intermediate ones in needle-leaved forests, and the lowest ones in shrublands and nonirrigated crops. Particularly, sites B, D, F, and J (much less pronounced) reproduce well the years characterized by low vegetation activity due to drought episodes (e.g., 2005, 2012, and 2015). In fact, site F provides the highest GPP and NPP decrease throughout the period considered (see Q GPP and Q NPP in Table 3). This effect is also identified in the annual NPP values for the evergreen site F, whereas a flat trend dominates in the rest of the forest sites.

Change Detection Analysis
GPP values, where the highest GPP values are found in broadleaved forests, the intermediate ones in needle-leaved forests, and the lowest ones in shrublands and nonirrigated crops. Particularly, sites B, D, F, and J (much less pronounced) reproduce well the years characterized by low vegetation activity due to drought episodes (e.g., 2005, 2012, and 2015). In fact, site F provides the highest GPP and NPP decrease throughout the period considered (see QGPP and QNPP in Table 3). This effect is also identified in the annual NPP values for the evergreen site F, whereas a flat trend dominates in the rest of the forest sites.  Table 1.
The normalized trends of GPP (QGPP) and NPP (QNPP) and annual averages are presented in Table 3 together with their correlations between SPI and T annual components. The grassland (B), shrubland (C), rainfed crop (D), and low-altitude evergreen needleleaved forest (I) sites are the most sensitive to annual precipitation. A strong negative correlation with temperature is also observed for sites B, C, D, and I, as opposed to the positive correlation found for deciduous broadleaved forest sites (G and H). The poorest correlations are obtained for the irrigated crop (site E) and the sparse vegetation (A). Site A is practically a desert and has the lowest mean annual GPP and the lowest interannual GPP variation. In case of NPP, a similar behavior is found for all forest sites with a notable decrease of biomass (QNPP < -0.28 x 10 -2 d -1 ) shown by the low-altitude evergreen forest  Table 1.
The normalized trends of GPP (Q GPP ) and NPP (Q NPP ) and annual averages are presented in Table 3 together with their correlations between SPI and T annual components. The grassland (B), shrubland (C), rainfed crop (D), and low-altitude evergreen needleleaved forest (I) sites are the most sensitive to annual precipitation. A strong negative correlation with temperature is also observed for sites B, C, D, and I, as opposed to the positive correlation found for deciduous broadleaved forest sites (G and H). The poorest correlations are obtained for the irrigated crop (site E) and the sparse vegetation (A). Site A is practically a desert and has the lowest mean annual GPP and the lowest interannual GPP variation. In case of NPP, a similar behavior is found for all forest sites with a notable decrease of biomass (Q NPP < -0.28 x 10 -2 d -1 ) shown by the low-altitude evergreen forest site (I) and evergreen broadleaved forest (F) site. A strong correlation between NPP and annual precipitation is observed for these sites, I and F. Table 3. GPP and NPP normalized trends (in d -1 ) (Q GPP and Q NPP ) and average annual values for (in kg m -2 y -1 ) the sites in Table 1. Correlations of annual GPP and NPP with interannual SPI and mean annual T are included. The correlation is derived considering anomaly time series for GPP and NPP. R > |0.6| in bold (p-value < 0.05).   (Figure 5d) are also computed. The images reveal that larger errors are observed for greater changes, and the opposite is also true. Therefore, there has not been a spatial pattern found in the relative image, and relative errors below 10% (absolute error < 2 × 10 −4 ) are observed for the major part of the study area. this figure, increasing productivity is observed in areas 1 and 2. Conversely, the negative trends in the GPP come with an associated rainfall diminution along time (Figure 6a). As an example, in eastern Spain, croplands at the south of the Teruel province (area 3 in Figure  5a) and north of the Castilla y León region (area 5 in Figure 5a) can be observed. This behavior is also observed in site F (national park of Los Alcornocales in Cádiz province, area 6 in Figure 5a), with a decrease in GPP and NPP, as supported by results shown in Table 3.  White color refers to pixels with a slope that is very small or pixels with a nonsignificant Mann-Kendall test value (0 value). For NPP, white color also refers to non-forest pixels. To better understand the nature of observed changes, Figure 7 identifies areas with positive and negative GPP trends based on different precipitation and temperature behavior. Thus, spot areas of natural vegetation (Figure 7a) with negative production trends and potentially affected by long-and short-term degradation processes (dark blue, light blue, and red colors) are located in Almería, Murcia, and València provinces. Areas with a negative GPP trend and rainfall patterns shifting towards conditions of increasing aridity (e.g., Almería region) [53] are those associated with QGPP < 0, QT < 0, and QSPI < 0 (dark blue color). Regions with negative GPP trends values in spite of increasing rainfall, that should White color refers to pixels with a slope that is very small or pixels with a nonsignificant Mann-Kendall test value (0 value). For NPP, white color also refers to non-forest pixels. Figure 5a demonstrates a positive GPP trend in northeastern and western Spain. In this figure, increasing productivity is observed in areas 1 and 2. Conversely, the negative trends in the GPP come with an associated rainfall diminution along time ( Figure 6a). As an example, in eastern Spain, croplands at the south of the Teruel province (area 3 in Figure 5a) and north of the Castilla y León region (area 5 in Figure 5a) can be observed. This behavior is also observed in site F (national park of Los Alcornocales in Cádiz province, area 6 in Figure 5a), with a decrease in GPP and NPP, as supported by results shown in Table 3.

Site
To better understand the nature of observed changes, Figure 7 identifies areas with positive and negative GPP trends based on different precipitation and temperature behavior. Thus, spot areas of natural vegetation (Figure 7a) with negative production trends and potentially affected by long-and short-term degradation processes (dark blue, light blue, and red colors) are located in Almería, Murcia, and València provinces. Areas with a negative GPP trend and rainfall patterns shifting towards conditions of increasing aridity (e.g., Almería region) [53] are those associated with Q GPP < 0, Q T < 0, and Q SPI < 0 (dark blue color). Regions with negative GPP trends values in spite of increasing rainfall, that should favor the production rise unless the canopy is in a very degraded condition, are selected based on Q GPP < 0, Q T < 0, and Q SPI > 0 values (light blue color). Positive production trends along with positive precipitation variations mainly correspond to most mountain areas located in northeastern and northwestern Spain (dark green). Figure 7b exhibits the cropland areas changes. In spite of the precipitation trend, large crop areas with positive productivity (Q GPP > 0) are mainly concentrated on irrigated and nonirrigated lands in the Guadalquivir and Ebro basins (light and dark green colors in northwestern and southern Spain). A productivity increase is observed in the Ebro basin (Q GPP > 0 and Q SPI > 0), whereas a productivity loss and nonfavorable precipitation conditions are found in croplands (Q GPP < 0 and Q SPI < 0 values located in areas 7, 8, and 9 in Figure 7b). White color refers to pixels with a slope that is very small or pixels with a nonsignificant Mann-Kendall test value (0 value). For NPP, white color also refers to non-forest pixels. To better understand the nature of observed changes, Figure 7 identifies areas with positive and negative GPP trends based on different precipitation and temperature behavior. Thus, spot areas of natural vegetation (Figure 7a) with negative production trends and potentially affected by long-and short-term degradation processes (dark blue, light blue, and red colors) are located in Almería, Murcia, and València provinces. Areas with a negative GPP trend and rainfall patterns shifting towards conditions of increasing aridity (e.g., Almería region) [53] are those associated with QGPP < 0, QT < 0, and QSPI < 0 (dark blue color). Regions with negative GPP trends values in spite of increasing rainfall, that should favor the production rise unless the canopy is in a very degraded condition, are selected based on QGPP < 0, QT < 0, and QSPI > 0 values (light blue color). Positive production trends along with positive precipitation variations mainly correspond to most mountain areas located in northeastern and northwestern Spain (dark green). Figure 7b exhibits the cropland areas changes. In spite of the precipitation trend, large crop areas with positive productivity (QGPP > 0) are mainly concentrated on irrigated and nonirrigated lands in the Guadalquivir and Ebro basins (light and dark green colors in northwestern and southern Spain). A productivity increase is observed in the Ebro basin (QGPP > 0 and QSPI > 0),

Discussion
This paper has explored the dynamics of GPP and NPP variables in Peninsular Spain, which is located in the mid-latitude ecotone (MLE) [54], the zone between 30°-60° latitude that comprises a transition belt between forests (temperate forest in Spain) and drylands. GPP and NPP are variables linked to carbon fluxes and biomass production, which can be used to analyze the ecosystem state. Thus, the analysis of GPP and NPP time series on a

Discussion
This paper has explored the dynamics of GPP and NPP variables in Peninsular Spain, which is located in the mid-latitude ecotone (MLE) [54], the zone between 30 • -60 • latitude that comprises a transition belt between forests (temperate forest in Spain) and drylands. GPP and NPP are variables linked to carbon fluxes and biomass production, which can be used to analyze the ecosystem state. Thus, the analysis of GPP and NPP time series on a long-term basis may be useful to detect degradation processes and loss of ecosystem services as a consequence of environmental changes and social pressures that characterize the MLE.

Spatial Patterns of Vegetation
The spatial distribution and magnitude of annual GPP in this study were consistent with the estimates of other studies [30,55]. Forests annual NPP followed a similar pattern to the one exhibited by GPP, but its magnitude is reduced to more than half its value because the methodology used to simulate NPP in this study only accounts for the trees in the pixel, while GPP accounts for all vegetation included in the pixel (see Section 2). This reduction is especially noticed for needle-leaved forests, probably explained as an underestimation by the methodology which was recently reported [46]. Previous studies stated that spatial variations of GPP and NPP are mainly controlled by climate conditions, vegetation types, and their spatial distribution [56,57]. Thus, broadleaved forests constitute the main potential natural vegetation of carbon storage over much of temperate Europe, the strategy of which is to deal with the lack of sunlight and cold temperatures in winter. This causes broadleaved trees to usually exhibit the maximum peak in summer due to the highest PAR, and undergo consequent loss of leaves in winter. The estimates of this study were also able to capture the intermediate annual GPP and NPP values given by evergreen forests (site F) in spite of remaining green all year. The smallest GPP was observed for sparse vegetation due to the small fraction of the pixel actually covered by vegetation.

Change Detection Analysis
The temporal analysis based on daily TS offered an opportunity to infer a higher accuracy in the trend estimates since more significant pixels are obtained as compared to annual resolution studies. This can be checked with GPP: the combination of daily interannual time series (derived from the MRA-WT) with the Mann-Kendall and the Theil-Sen slope methods provided 72% of the considered pixels with statistically significant trend in this study, as opposed to only 17% found in a previous study over the same area using only annual series [58], even if three more years were used (18 instead of 15).
Trend results are consistent with similar studies performed over different periods [8]. The associated errors with Q GPP and Q NPP reveal that the precision does not depend on the vegetation type cover, since a specific spatial pattern was not found. Furthermore, larger significant negative changes and associated errors are observed, especially in southeastern Spain, eastern Mediterranean coastland, and central Spain, as was already found in [14] for 1989-2002, where a remarkable loss of productivity in southeastern and eastern Spain was observed. Particularly, the sparse natural vegetation in Almería province (almost a desert) is the most affected land cover with a significant loss of production during this period. Positive GPP trends are generally coincident with positive SPI trends (see Figures 5a and 6a), but in some areas they may be partly assured by the abandonment of old cultivated fields in mountain areas that have been transformed to natural vegetation [59]. For example, positive patches in the Ebro basin (Figures 5a and 7b) may be the consequence of the abandonment of old irrigated fields due to relocation and socioeconomic changes benefiting from the positive precipitation trend (Q SPI > 0 in Figure 6a) for the same period [59]. Most of these areas are classified as arable lands and grasslands (http://sigpac.aragon.es/visor/, accessed on 13 January 2022) with a high response to favorable meteorological conditions. Moreover, a high overlap with the conservative areas included in the Natura 2000 ecological network (https://www.miteco.gob.es/, accessed on 13 January 2022) is found. It comprises up to 40% of agricultural lands, predominantly marginal areas with low productivity. Conversely, the productivity increase in the Guadalquivir basin (light green color in Figure 7b) coexists with unfavorable precipitation conditions (Figure 6a). This could be attributed to the remarkable expansion in irrigated areas and modernization of irrigation systems [60]. In other areas, vegetation recovery from forest fires that occurred in the beginning of the period may explain this behavior. This is the case of a 34,000 ha area affected by the 2004 forest fire in Riotinto (Huelva province) and 22,000 ha burnt in 2005 in La Riba (Guadalajara province), as well as croplands and needle-leaved forests in Murcia and Comunitat Valenciana.
Precipitation is one of the most important meteorological factors that drive GPP in natural vegetation. Negative trends in GPP demonstrate its close relationship with precipitation ( Figure 5a). A similar pattern can be observed in the SPI trend (Figure 6a). These observations are particularly clear in arid and semiarid areas where cloudy skies are not frequent, and precipitation is a controlling factor of vegetation phenology and productivity [61]. In this case, the lack of water provision ( Figure 6a) has a direct impact on carbon absorption (mainly controlled by PAR) and carbon assimilation after autotrophic respiration. Negative trends in agricultural areas (areas 7, 8, and 9 in Figure 7b) mainly belong to regions where average rainfall decreased and its seasonal distribution was altered. For instance, La Mancha wetlands (area 7) has been widely recognized to suffer from the depletion of the groundwater due to an excessive water extraction from the aquifers [62] and an intense use of agricultural exploitation practices. At local scale, the strong correlation between GPP (r SPI-GPP > 0.59) and NPP (r SPI-GPP > 0.53) with SPI found for almost all sites (Table 3) reinforces the statement that precipitation is a key factor of vegetation dynamics in temperate regions [57]. A remarkable reduction in photosynthetic activity and carbon assimilation (Q GPP (10 −2 ) < −0.02 d −1 ; Q NPP (10 −2 ) < −0.23 d −1 ) was observed on interannual scale at sites C (shrublands), D (rainfed crops), F (evergreen broadleaved forest), and I (low-altitude evergreen needle-leaved forest). In these cases, a negative correlation with temperature may possibly be a determinant factor. Increasing temperatures can cause high evaporation and larger vapor pressure deficit, which in turn creates water deficit and, subsequently, a negative impact on optimum plant growth which seems to more directly affect carbon fluxes [63]. Particularly, site F (Mediterranean evergreen cork oak) was reported to experience a considerable increase of mean annual air temperature and rainfall extremes during recent decades, projecting a higher frequency of droughts and intense rain events. Productivity losses reported in this study are supported by a generalized decline of the cork oak in the Mediterranean basin found in different studies [64]. Site F shows a severe reduction in annual GPP in 2005, 2012, 2013, and 2017. As reported by the Spanish Meteorological Agency (www.aemet.es, accessed on 13 January 2022), 2005 has been recognized as the driest year, followed by 2017 and 2012 with 40%, 27%, and 15% less annual precipitation than the mean value, respectively. In this regard, [63] reported that drought is the driving factor of global NPP diminution. Moreover, different causes, such as the increase in the outbreaks of tree disease due to more favorable conditions for pathogen development, extreme overgrazing by livestock leading to a reduced regeneration rate, and an excessive development of forest resources [64] may also explain these findings. For site I (conifer forest species, mainly Pinus nigra ssp.), a reduction in productivity was observed (Q GPP (10 −2 ) = −0.021 ± 0.015 d −1 ; Q NPP (10 −2 ) = −0.23 ±0.04 d −1 ), which was more important in the production of new biomass. The large positive correlation with SPI and the large negative correlation with air temperature may indicate that this ecosystem is particularly vulnerable to climate variations. Vulnerability of zonal ecosystems, particularly forests, over the MLE is very high given its significant role in carbon and hydrological cycling studies [65]. Evidence of major changes for Mediterranean conifer species located at low elevations was also found in [66].
Site A (sparse vegetation) presented the largest negative GPP trend (Q GPP (10 −2 ) = −0.93 ± 0.09 d −1 ) and a poor correlation with air temperature. The trend is supported by the land aridity observed in recent decades due to water scarcity, temperature, and radiation factors [8]. The small correlation between GPP and air temperature (r GPP-T = 0.12) evidences the precipitation (r GPP-SPI = 0.59) as a limiting growth factor compared to air temperature. In sparse vegetation areas, vegetation conduces to the expansion of root system for the increasing absorption of minerals and water in an attempt to increase resilience. Figure 6a revealed that a very advanced degradation process may be explained by the severely decreased precipitation in the studied period. As reported by [67], less precipitation in aridity areas is compensated in the long term by a proportional loss of density in vegetation cover.
Increase of productivity (positive GPP (NPP) trend) was observed in sites B (grasslands), E (irrigated crops), G (low-altitude DBF), and H (high-altitude DBF). While site B revealed a positive correlation between GPP and SPI and negative correlation between GPP and air temperature, site G provided a positive correlation between GPP (NPP) and both SPI and air temperature. Site B belongs to a Mediterranean savanna-like ecosystem with scattered trees (open holm oak woodland) over annual grassland, where the favorable condition of precipitation observed in Figure 6a may account for the increasing of GPP (reaching GPP values up to 5 g m -2 d -1 ). Site G belongs to a broadleaved forest where favorable precipitation and air temperature conditions lead to the predominance of positive GPP (NPP) trends rather than negative ones ( Figure 6 and Table 3). The smallest correlations between GPP (NPP) and SPI were found in sites E and H. Despite these low correlations, site E showed positive GPP trend mainly controlled by irrigation systems. Although this site exhibited productivity values up to 15 g m -2 d -1 , it may be slightly affected by water restrictions due to water resources scarcity as it was observed in lower peaks over drought events (e.g., 2015). In this case, the response of the crop to this situation would depend on the drought magnitude and resilience of the affected area. Site H showed a higher correlation between GPP (NPP) and air temperature than between GPP (NPP) and SPI. It is probable that warming is amplified with elevation resulting in rapid changes in temperature, humidity, and water in mountains areas [68]. However, although these areas may be affected by drought periods, the development of the vegetation generally continues in spring and summer because enough water is still available. Thus, significant positive trends may be found (Q GPP (10 −2 ) = 0.161 ± 0.006 d −1 ; Q NPP (10 −2 ) = 0.111±0.005 d −1 ).
To better understand the effect of these long-term changes, other magnitudes may be also examined regarding other significant information on carbon fluxes. This is the case of the carbon use efficiency (CUE). It is an interesting variable that describes how efficiently plants incorporate the carbon fixed during photosynthesis into biomass gain [69]. It can be calculated as the ratio between NPP and GPP. In a previous study [70], the CUE was analyzed for the Peninsular Spain study area throughout eight years, 2005-2012. Results showed that CUE exhibited a positive correlation with precipitation and a negative correlation with temperature in most ecosystems. In this study, GPP refers to the whole ecosystem, whereas NPP refers only to the trees and is only obtained for forest canopies; then, an analysis of CUE cannot be addressed to compare with previous results. However, considering the annual values of our NPP and GPP (Table 3) as a proxy of CUE, it is observed that, effectively, the correlation between CUE and precipitation is positive (for all the sites except for J), and the correlation between CUE and temperature is negative (for all sites except G), thus confirming a general tendency of CUE to decrease when the ecosystem conditions change towards aridity.

Conclusions
The methodology proposed to infer long-term vegetation changes (MRA-WT) from 15-year daily GPP and NPP 1 km spatial resolution series over Peninsular Spain demonstrated its effectiveness in a multilevel decomposition time series at different temporal scales.
The main key findings of our study are as follows: i. Daily GPP and NPP time series demonstrated to be a high quality input for the temporal analysis, specifically for water stress characterization in sparse vegetation areas (obtained the same precision, 10%, as in dense vegetation areas).
ii. A higher temporal resolution offered an opportunity to improve the accuracy in the trend estimates since more significant pixels (72%) were obtained as compared to annual resolution studies (17%). iii. A negative clear agreement between GPP and precipitation was observed, particularly in southeastern Spain, eastern Mediterranean coastland, and central Spain. An increase in temperature was shown to favor the carbon assimilation of deciduous broadleaved forest in the north of Spain. iv. Evidence of forest vulnerability was confirmed, particularly in Mediterranean conifer species located at low elevations.
As further research, the quantification of limiting factors, such as water stress in the vegetation response and the production of biomass, is particularly relevant. In this sense, the interannual CUE characterization of this region would be of great interest, as would the water used efficiency (WUE), which is defined as the ratio of carbon gain to water consumption (i.e., evapotranspiration) and quantifies the rate of carbon uptaken per unit of water loss.