Next Article in Journal
Evaluation of Satellite-Derived Estimates of Lake Ice Cover Timing on Linnévatnet, Kapp Linné, Svalbard Using In-Situ Data
Next Article in Special Issue
A Decrease in the Daily Maximum Temperature during Global Warming Hiatus Causes a Delay in Spring Phenology in the China–DPRK–Russia Cross-Border Area
Previous Article in Journal
Uncertainty-Guided Depth Fusion from Multi-View Satellite Images to Improve the Accuracy in Large-Scale DSM Generation
Previous Article in Special Issue
Integrating Multi-Source Remote Sensing to Assess Forest Aboveground Biomass in the Khingan Mountains of North-Eastern China Using Machine-Learning Algorithms
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Remote Sensing
Volume 14
Issue 6
10.3390/rs14061310
remotesensing-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

by
Beatriz Martínez
*,
Sergio Sánchez-Ruiz
,
Manuel Campos-Taberner
,
F. Javier García-Haro
and
M. Amparo Gilabert
Environmental Remote Sensing Group, Departament de Física de la Terra i Termodinàmica, Facultat de Física, Universitat de València, 46100 Burjassot, Spain
*
Author to whom correspondence should be addressed.
Remote Sens. 2022, 14(6), 1310; https://doi.org/10.3390/rs14061310
Submission received: 14 January 2022 / Revised: 4 March 2022 / Accepted: 4 March 2022 / Published: 8 March 2022
(This article belongs to the Special Issue Remote Sensing Application for Promoting Ecosystem Services and Land Degradation Management in Mid-Latitude Ecotone (MLE))
Browse Figures
Versions Notes

Abstract

:
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 fluxes 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 specific 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 significant pixels. At pixel level, larger significant 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 fluxes trends were correlated with precipitation. The temporal analysis based on daily TS demonstrated an increase of accuracy in the trend estimates since more significant pixels were obtained as compared to annual resolution studies (72% as to only 17%).

1. 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 then, long time series of EO-based products from different models have been provided at different spatial and temporal scales. This is the case of the already mentioned MODIS products, MOD17, at different temporal and spatial resolutions [4], the daily Soil Moisture Active Passive (SMAP) Level-4 carbon (L4_C) product at 9 km [5] from NASA, or Gross Dry Matter Productivity (GDMP) product at 1 km [6] from the Copernicus Global Land Service. All of them are based on the production efficiency model (PEM) concept proposed by Monteith [7].
Most studies on the response of vegetation dynamics to climate change in Spain have focused until now on the NDVI time series (TS) [8] due to its well-recognized lineal relationship with photosynthesis as well as with GPP and NPP across diverse ecosystems [9]. Although the robustness of vegetation indices (VIs) along with the availability of long TS have led to a vast number of global carbon-climate studies [9,10], VIs only have a moderate capacity to predict ecosystem structural and functional attributes. It is expected that GPP and NPP provide more information regarding the vegetation status thanks to the PAR and the light-use efficiency variables, which have shown a more rapid response than fAPAR to different environmental factors related to energy balance, water availability, and nutrient levels [11].
However, the carbon flux models mentioned above do not offer completely satisfactory results in Spain; some because of how they implement the water stress (MOD17 and L4_C), and others because they do not implement it (GDMP). Particularly, it was demonstrated that a different characterization of water stress, such as the one used in the current study, leads to more accurate results in the region [12] (see a summary of the methodology in Section 2.3) and, therefore, provides a better quality input for the temporal analysis. Moreover, a higher temporal variability may offer the possibility of detecting more accurately the presence of cloudy skies that leads to a carbon decrease related to the irradiation reduction and not to the vegetation type. Hence, inappropriate daily GPP and NPP characterization may result in inaccurate global GPP and NPP values at annual scales.
The development of effective methodologies for the analysis of TS is one of the most important challenging issues for the remote sensing community due to the dynamic nature of terrestrial ecosystems [13,14]. Many studies have tried to identify the best performing method among the different approaches proposed in the literature (including statistical, time–frequency-based, etc.). Based on previous comparative studies, ensemble empirical mode decomposition (EEMD) [15], wavelet transform (WT) [14], break for additive trend and season (BFAST) [13], and detecting breakpoints and estimating segments in trend (DBEST) [16] have proved their effectiveness for vegetation change analysis using normalized difference vegetation index (NDVI) TS [15].
In this study, the multiresolution analysis based on the wavelet transform (MRA-WT) was used because of its results during the last decades as a time–frequency analysis tool for complex nonstationary signals in several fields (e.g., environment, medicine, finance) [17]. Particularly, it was useful in the study of nonstationary TS for vegetation temporal dynamics [14,18,19].
The interannual variability of vegetation dynamics plays a key role in the description of year-to-year variations of plant photosynthesis and production, determining land carbon uptake and ecosystem capacity to be productive [20]. The interest in ecosystem productivity responses to climate change grew substantially in the research community during the last decade [8,21]. The relationships of vegetation interannual variability and conservative climatic factors, particularly air temperature and precipitation, have been studied in other regions [19,22].
The main objective of this study is to analyze (1) the spatial and interannual variations of carbon fluxes in Peninsular Spain across 15 years from 2004 to 2018 using daily GPP and NPP EO-based estimates and (2) their relationships with changes in air temperature and precipitation, since the hydrological response to climate change over Spain due to the increased air temperature and decreased precipitation is already confirmed [8].

2. 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.

2.1. 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.

2.2. 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).

2.3. GPP

The PEM used to estimate the daily GPP (g m–2 d–1) images is
GPP = PAR fAPAR εmax εW εT
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 fAPAR (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 SEVIRI/MSG fAPAR product (LSA-425) [32]. In this study, the fAPAR TS were filtered and reconstructed using an optimized locally weighted regression and smoothing scatterplots (LOESS) method [33].
εW is the water stress coefficient CWS 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):
CWS = 0.5 + 0.5 AET/PET
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 Rg (kJ m–2 d–1), as given by the LSA-203 product, and the daily T (in °C):
PET = Rg (0.025 T + 0.08)/2450
ε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.

2.4. NPP

Forest NPP was obtained based on the combination of the GPP calculated through equation (1), and the growth (RG) and maintenance (RM) 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
NPP = f (GPP − RG) − v RM
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
f = [1 − exp (–v LAI)]/[1 − exp (–LAI)]
and
v = 50 GSV ρ C/σ
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 m3/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].

3. 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 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].
Ten sites representative of the study area (Table 1) were selected to examine their seasonal variations. One site was chosen for each vegetation type indicated in Figure 1c,d. All sites present ecoclimatic features representative of the corresponding vegetation type except site F, which presents higher air temperature.
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.

4. 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).

4.1. 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, 2j (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
g (t) = Am (t) + Σmj =1 Dj (t)
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 vc = 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: A1 (1.5 days), A2 (3 days), A3 (6 days), A4 (12 days), A5 (24 days), A6 (48 days), A7 (96 days), A8 (192 days), A9 (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 (A9 and A4, 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 (V4–9). The plot on the bottom shows the interannual (A9) component of the original TS.

4.2. 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 (QGPP), NPP (QNPP), SPI (QSPI), and T (QT) 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 (H0) of this test (no trend exists) may not be rejected, whereas a significant test value (1 value) assumes H0 is rejected, and a trend is considered. H0 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.
QGPP, QNPP, and QT 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 QSPI 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 (A9 for GPP, NPP, and T, and A4 for SPI).

5. Results

5.1. Spatial Patterns of Vegetation

Figure 3 shows the spatial distribution of average annual GPP and NPP for the study period (2004–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).

5.2. Change Detection Analysis

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 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.
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 needle-leaved 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 site (I) and evergreen broadleaved forest (F) site. A strong correlation between NPP and annual precipitation is observed for these sites, I and F.
Figure 5a,b illustrate the GPP and NPP trend images for the period 2004–2018 (QGPP and QNPP). Red values refer to negative changes, whereas green values are indicative of positive variations. No significant variation during the 15-year period is observed in 28% (5%) of the GPP (NPP) values. Highly significant GPP (NPP) decrease is observed in 16% (25%) of the Spain pixels. Highly significant GPP (NPP) increase is observed in 56% (70%) of the Spain pixels. The associated errors with QGPP (Figure 5c) and QNPP (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.
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.
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), whereas a productivity loss and nonfavorable precipitation conditions are found in croplands (QGPP < 0 and QSPI < 0 values located in areas 7, 8, and 9 in Figure 7b).

6. 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.

6.1. 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.

6.2. 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 QGPP and QNPP 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 Figure 5a and Figure 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 (Figure 5a and Figure 7b) may be the consequence of the abandonment of old irrigated fields due to relocation and socioeconomic changes benefiting from the positive precipitation trend (QSPI > 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 (rSPI-GPP > 0.59) and NPP (rSPI-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 (QGPP(10−2) < −0.02 d−1; QNPP(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 (QGPP(10−2) = −0.021 ± 0.015 d−1; QNPP(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 (QGPP(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 (rGPP-T = 0.12) evidences the precipitation (rGPP-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 (QGPP(10−2) = 0.161 ± 0.006 d−1; QNPP(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.

7. 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.

Author Contributions

Conceptualization, methodology, software, formal analysis, investigation, resources, data curation, writing, visualization, supervision, project administration, and funding acquisition: B.M., S.S.-R., M.C.-T., F.J.G.-H. and M.A.G. All authors have read and agreed to the published version of the manuscript.

Funding

This research was supported thanks to grant CGL2016-75238-R (project ESCENARIOS) funded by Ministerio de Ciencia e Innovación (MCIN) and grant PID2020-118036RB-I00 (project ECCE EO) funded by MCIN/AEI/10.13039/501100011033 and by “ERDF A way of making Europe”. The European Organization for the Exploitation of Meteorological Satellites (EUMETSAT) with grant number LSA SAF CDOP-3 has also supported this study.

Acknowledgments

Part of the data was kindly provided by the Spanish Meteorological Agency (www.aemet.es). The Third Spanish National Forest Inventory data were freely downloaded from the Spanish Ministry of Ecological Transition (https://www.miteco.gob.es/es/biodiversidad/servicios/banco-datosnaturaleza/informacion-disponible/ifn3.aspx, accessed on 13 January 2022). We thank Jeff Burkey for implementing the code to compute Mann–Kendall test with the Theil–Sen slope in MATLAB. We thank Taesam Lee for implementing the code to compute SPI in MATLAB.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

References

  1. Xiao, J.; Chevallier, F.; Gomez, C.; Guanter, L.; Hicke, J.A.; Huete, A.R.; Ichii, K.; Ni, W.; Pang, Y.; Rahman, A.F.; et al. Remote sensing of the terrestrial carbon cycle: A review of advances over 50 years. Remote Sens. Environ. 2019, 233, 111383. [Google Scholar] [CrossRef]
  2. Waring, H.R.; Running, S.W. Forest Ecosystems, 3rd ed.; Academic Press: San Diego, CA, USA, 2007. [Google Scholar]
  3. Heinsch, F.A.; Zhao, M.; Running, S.W.; Kimball, J.S.; Nemani, R.R.; Davis, K.J.; Bolstad, P.V.; Cook, B.D.; Desai, A.R.; Ricciuto, D.M.; et al. Evaluation of remote sensing based terrestrial productivity from MODIS using regional tower eddy flux network observations. IEEE Trans. Geosci. Remote Sens. 2006, 44, 1908–1923. [Google Scholar] [CrossRef] [Green Version]
  4. Running, S.; Zhao, M. User’s Guide Daily GPP and Annual NPP (MOD17A2H/A3H) and Year-end GapFilled (MOD17A2HGF/A3HGF) Products NASA Earth Observing System MODIS Land Algorithm (For Collection 6.1); University of Montana: Missoula, MT, USA, 2021; pp. 1–33. [Google Scholar]
  5. Jones, L.A.; Kimball, J.S.; Reichle, R.H.; Madani, N.; Glassy, J.; Ardizzone, J.V.; Colliander, A.; Cleverly, J.; Desai, A.R.; Eamus, D.; et al. The SMAP Level 4 Carbon Product for Monitoring Ecosystem Land-Atmosphere CO2 Exchange. IEEE Trans. Geosci. Remote Sens. 2017, 55, 6517–6532. [Google Scholar] [CrossRef]
  6. Cglops-1 Product User Manual Dry Matter Productivity (Dmp) Gross Dry Matter Productivity (Gdmp) Collection 1 km Version 2; Copernicus Global Land Operations: Brussels, Belgium, 2019.
  7. Monteith, J.L. Solar Radiation and Productivity in Tropical Ecosystems. J. Appl. Ecol. 1972, 9, 747–766. [Google Scholar] [CrossRef] [Green Version]
  8. Vicente-Serrano, S.M.; Martín-Hernández, N.; Reig, F.; Azorin-Molina, C.; Zabalza, J.; Beguería, S.; Domínguez-Castro, F.; El Kenawy, A.; Peña-Gallardo, M.; Noguera, I.; et al. Vegetation greening in Spain detected from long term data (1981–2015). Int. J. Remote Sens. 2020, 41, 1709–1740. [Google Scholar] [CrossRef]
  9. Tucker, C.J.; Slayback, D.A.; Pinzon, J.E.; Los, S.O.; Myneni, R.B.; Taylor, M.G. Higher northern latitude normalized difference vegetation index and growing season trends from 1982 to 1999. Int. J. Biometeorol. 2001, 45, 184–190. [Google Scholar] [CrossRef]
  10. Saleska, S.R.; Wu, J.; Guan, K.; Araujo, A.C.; Huete, A.; Nobre, A.D.; Restrepo-Coupe, N. Dry-season greening of Amazon forests. Nature 2016, 531, E4–E5. [Google Scholar] [CrossRef]
  11. Garbulsky, M.F.; Peñuelas, J.; Papale, D.; Ardö, J.; Goulden, M.L.; Kiely, G.; Richardson, A.D.; Rotenberg, E.; Veenendaal, E.M.; Filella, I. Patterns and controls of the variability of radiation use efficiency and primary productivity across terrestrial ecosystems. Glob. Ecol. Biogeogr. 2010, 19, 253–267. [Google Scholar] [CrossRef]
  12. Gilabert, M.A.; Moreno, A.; Maselli, F.; Martínez, B.; Chiesi, M.; Sánchez-Ruiz, S.; García-Haro, F.J.; Pérez-Hoyos, A.; Campos-Taberner, M.; Pérez-Priego, O.; et al. Daily GPP estimates in Mediterranean ecosystems by combining remote sensing and meteorological data. ISPRS J. Photogramm. Remote Sens. 2015, 102, 184–197. [Google Scholar] [CrossRef]
  13. Verbesselt, J.; Hyndman, R.; Newnham, G.; Culvenor, D. Detecting trend and seasonal changes in satellite image time series. Remote Sens. Environ. 2010, 114, 106–115. [Google Scholar] [CrossRef]
  14. Martínez, B.; Gilabert, M.A. Vegetation dynamics from NDVI time series analysis using the wavelet transform. Remote Sens. Environ. 2009, 113, 1823–1842. [Google Scholar] [CrossRef]
  15. Kong, Y.L.; Meng, Y.; Li, W.; Yue, A.Z.; Yuan, Y. Satellite image time series decomposition based on EEMD. Remote Sens. 2015, 7, 15583–15604. [Google Scholar] [CrossRef] [Green Version]
  16. Jamali, S.; Jönsson, P.; Eklundh, L.; Ardö, J.; Seaquist, J. Detecting changes in vegetation trends using time series segmentation. Remote Sens. Environ. 2015, 156, 182–195. [Google Scholar] [CrossRef]
  17. Rhif, M.; Abbes, A.B.; Farah, I.R.; Martínez, B.; Sang, Y. Wavelet transform application for/in non-stationary time-series analysis: A review. Appl. Sci. 2019, 9, 1345. [Google Scholar] [CrossRef] [Green Version]
  18. Huang, M.; Piao, S.; Sun, Y.; Ciais, P.; Cheng, L.; Mao, J.; Poulter, B.; Shi, X.; Zeng, Z.; Wang, Y. Change in terrestrial ecosystem water-use efficiency over the last three decades. Glob. Chang. Biol. 2015, 21, 2366–2378. [Google Scholar] [CrossRef] [PubMed]
  19. Martínez, B.; Gilabert, M.A.; García-Haro, F.J.; Faye, A.; Meliá, J. Characterizing land condition variability in Ferlo, Senegal (2001-2009) using multi-temporal 1-km Apparent Green Cover (AGC) SPOT Vegetation data. Glob. Planet. Chang. 2011, 76, 152–165. [Google Scholar] [CrossRef]
  20. Le Quéré, C.; Raupach, M.R.; Canadell, J.G.; Marland, G.; Bopp, L.; Ciais, P.; Conway, T.J.; Doney, S.C.; Feely, R.A.; Foster, P.; et al. Trends in the sources and sinks of carbon dioxide. Nat. Geosci. 2009, 2, 831–836. [Google Scholar] [CrossRef]
  21. Baldocchi, D.; Chu, H.; Reichstein, M. Inter-annual variability of net and gross ecosystem carbon fluxes: A review. Agric. For. Meteorol. 2018, 249, 520–533. [Google Scholar] [CrossRef] [Green Version]
  22. Zhang, L.; Ren, X.; Wang, J.; He, H.; Wang, S.; Wang, M.; Piao, S.; Yan, H.; Ju, W.; Gu, F.; et al. Interannual variability of terrestrial net ecosystem productivity over China: Regional contributions and climate attribution. Environ. Res. Lett. 2019, 14, 014003. [Google Scholar] [CrossRef]
  23. Moreno, A.; Gilabert, M.A.; Martínez, B. Mapping daily global solar irradiation over Spain: A comparative study of selected approaches. Sol. Energy 2011, 85, 2072–2084. [Google Scholar] [CrossRef]
  24. WMO. Standardized Precipitation Index User Guide; WMO: Geneva, Switzerland, 2012. [Google Scholar]
  25. Pérez-Hoyos, A.; García-Haro, F.J.; San-Miguel-Ayanz, J. A methodology to generate a synergetic land-cover map by fusion of different: Land-cover products. Int. J. Appl. Earth Obs. Geoinf. 2012, 19, 72–87. [Google Scholar] [CrossRef]
  26. Farr, T.; Rosen, P.; Caro, E.; Crippen, R.; Duren, R.; Hensley, S.; Kobrick, M.; Paller, M.; Rodriguez, E.; Roth, L.; et al. The shuttle radar topography mission. Rev. Geophys. 2007, 45, 1–33. [Google Scholar] [CrossRef] [Green Version]
  27. IGN Documento Técnico SIOSE 2011; Ministerio de Fomento, Dirección general del Instituto Geográfico Nacional (IGN): Madrid, Spain, 2011.
  28. Sánchez-Ruiz, S.; Chiesi, M.; Fibbi, L.; Carrara, A.; Maselli, F.; Gilabert, M.A.M.A. Optimized Application of Biome-BGC for Modeling the Daily GPP of Natural Vegetation Over Peninsular Spain. J. Geophys. Res. Biogeosci. 2018, 123, 531–546. [Google Scholar] [CrossRef]
  29. Product User Manual Down-welling Surface Shortwave Flux (DSSF); LSA-SAF: Lisbon, Portugal, 2011.
  30. Martínez, B.; Gilabert, M.A.; Sánchez-Ruiz, S.; Campos-Taberner, M.; García-Haro, F.J.; Brümmer, C.; Carrara, A.; Feig, G.; Grünwald, T.; Mammarella, I.; et al. Evaluation of the LSA-SAF gross primary production product derived from SEVIRI/MSG data (MGPP). ISPRS J. Photogramm. Remote Sens. 2020, 159, 220–236. [Google Scholar] [CrossRef]
  31. Roujean, J.; Breon, F. Estimating PAR Absorbed by Vegetation from Bidirectional Reflectance Measurements. Remote Sens. Environ. 1995, 51, 375–384. [Google Scholar] [CrossRef]
  32. Algorithm Theoretical Basis Document for Vegetation Parameters (VEGA); LSA-SAF: Lisbon, Portugal, 2016.
  33. Moreno, Á.; García-Haro, F.; Martínez, B.; Gilabert, M.A. Noise Reduction and Gap Filling of fAPAR Time Series Using an Adapted Local Regression Filter. Remote Sens. 2014, 6, 8238–8260. [Google Scholar] [CrossRef] [Green Version]
  34. Potter, C.S.; Klooster, S.; Myneni, R.; Genovese, V.; Tan, P.N.; Kumar, V. Continental-scale comparisons of terrestrial carbon sinks estimated from satellite data and ecosystem modeling 1982-1998. Glob. Planet. Chang. 2003, 39, 201–213. [Google Scholar] [CrossRef]
  35. Maselli, F.; Papale, D.; Puletti, N.; Chirici, G.; Corona, P. Combining remote sensing and ancillary data to monitor the gross productivity of water-limited forest ecosystems. Remote Sens. Environ. 2009, 113, 657–667. [Google Scholar] [CrossRef] [Green Version]
  36. Sánchez-Ruiz, S.; Moreno, A.; Piles, M.; Maselli, F.; Carrara, A.; Running, S.; Gilabert, M.A. Quantifying water stress effect on daily light use efficiency in Mediterranean ecosystems using satellite data. Int. J. Digit. Earth 2017, 10, 623–638. [Google Scholar] [CrossRef]
  37. Jensen, M.E.; Haise, H.R. Estimating evapotranspiration from solar radiation. J. Irrig. Drain. Div. 1965, 89, 15–41. [Google Scholar] [CrossRef]
  38. Running, S.; Hunt, E. Generalization of a forest ecosystem process model for other biomes, BIOME-BGC, and an application for global-scale models. In Scaling Physiological Processes: Leaf to Globe: A volume in Physiological Ecology; University of Montana: Missoula, MT, USA, 1993; pp. 141–158. [Google Scholar] [CrossRef]
  39. White, M.A.; Thornton, P.E.; Running, S.W.; Nemani, R.R. Parameterization and Sensitivity Analysis of the BIOME–BGC Terrestrial Ecosystem Model: Net Primary Production Controls. Earth Interact. 2000, 4, 1–85. [Google Scholar] [CrossRef]
  40. Thornton, P.E.; Hasenauers, H.; White, M.A. Simultaneous estimation of daily solar radiation and humidity from observed temperature and precipitation: An application over complex terrain in Austria. Agric. For. Meteorol. 2000, 104, 255–271. [Google Scholar] [CrossRef]
  41. Chiesi, M.; Maselli, F.; Moriondo, M.; Fibbi, L.; Bindi, M.; Running, S.W. Application of BIOME-BGC to simulate Mediterranean forest processes. Ecol. Modell. 2007, 206, 179–190. [Google Scholar] [CrossRef]
  42. Maselli, F.; Chiesi, M.; Moriondo, M.; Fibbi, L.; Bindi, M.; Running, S.W. Modelling the forest carbon budget of a Mediterranean region through the integration of ground and satellite data. Ecol. Modell. 2009, 220, 330–342. [Google Scholar] [CrossRef]
  43. Sánchez-Ruiz, S.; Moreno-Martínez, Á.; Izquierdo-Verdiguier, E.; Chiesi, M.; Maselli, F.; Gilabert, M.A. Growing stock volume from multi-temporal landsat imagery through google earth engine. Int. J. Appl. Earth Obs. Geoinf. 2019, 83, 101913. [Google Scholar] [CrossRef]
  44. DGCN III Inventario Forestal Nacional 1997–2007; Dirección General de Conservación de la Naturaleza, Ministerio de Medio Ambiente: Madrid, Spain, 2006.
  45. Federici, S.; Vitullo, M.; Tulipano, S.; De Laurets, R.; Seufert, G. An approach to estimate carbon stocks change in forest carbon pools under the UNFCCC: The Italian case. Iforest—Biogeosci. For. 2008, 1, 86–95. [Google Scholar] [CrossRef] [Green Version]
  46. Sánchez-Ruiz, S.; Maselli, F.; Chiesi, M.; Fibbi, L.; Martínez, B.; Campos-Taberner, M.; García-Haro, F.J.; Gilabert, M.A. Remote sensing and bio-geochemical modeling of forest carbon storage in Spain. Remote Sens. 2020, 12, 1356. [Google Scholar] [CrossRef]
  47. Percival, D.B.; Walden, A.T. Wavelet Methods For Time Series Analysis; Cambridge University Press: Cambridge, UK, 2000. [Google Scholar]
  48. Ben Abbes, A.; Bounouh, O.; Farah, I.R.; de Jong, R.; Martínez, B. Comparative study of three satellite image time-series decomposition methods for vegetation change detection. Eur. J. Remote Sens. 2018, 51, 607–615. [Google Scholar] [CrossRef] [Green Version]
  49. Kendall, M.G. Rank Correlation Methods, 4th ed.; Charles Griffiin: London, UK, 1985. [Google Scholar]
  50. Mann, H.B. Non-parametric tests against trend. Econometrica 1945, 13, 245–259. [Google Scholar] [CrossRef]
  51. Sen, P.K. Estimates of the Regression Coefficient Based on Kendall’s Tau. J. Am. Stat. Assoc. 1968, 63, 1379–1389. [Google Scholar] [CrossRef]
  52. de Jong, R.; de Bruin, S.; de Wit, A.; Schaepman, M.E.; Dent, D.L. Analysis of monotonic greening and browning trends from global NDVI time-series. Remote Sens. Environ. 2011, 115, 692–702. [Google Scholar] [CrossRef] [Green Version]
  53. Paniagua, L.L.; García-Martín, A.; Moral, F.J.; Rebollo, F.J. Aridity in the Iberian Peninsula (1960–2017): Distribution, tendencies, and changes. Theor. Appl. Climatol. 2019, 138, 811–830. [Google Scholar] [CrossRef]
  54. Song, C.; Pietsch, S.A.; Kim, M.; Cha, S.; Park, E.; Shvidenko, A.; Schepaschenko, D.; Kraxner, F.; Lee, W.K. Assessing Forest Ecosystems across the Vertical Edge of the Mid-Latitude Ecotone Using the BioGeoChemistry Management Model (BGC-MAN). Forests 2019, 10, 523. [Google Scholar] [CrossRef] [Green Version]
  55. Gilabert, M.A.; Sánchez-Ruiz, S.; Moreno, A. Annual gross primary production from vegetation indices: A theoretically sound approach. Remote Sens. 2017, 9, 193. [Google Scholar] [CrossRef] [Green Version]
  56. Cao, M.; Prince, S.D.; Small, J.; Goetz, S.J. Remotely Sensed Interannual Variations and Trends in Terrestrial Net Primary Productivity 1981-2000. Ecosystems 2004, 7, 233–242. [Google Scholar] [CrossRef]
  57. Nemani, R.R.; Keeling, C.D.; Hashimoto, H.; Jolly, W.M.; Piper, S.C.; Tucker, C.J.; Myneni, R.B.; Running, S.W. Climate-driven increases in global terrestrial net primary production from 1982 to 1999. Science 2003, 300, 1560–1563. [Google Scholar] [CrossRef] [Green Version]
  58. Sánchez-Ruiz, S.; Martínez, B.; Campos-Taberner, M.; García-Haro, F.J.; Gilabert, M.A. Análisis de tendencia en la GPP anual sobre la España peninsular. In Proceedings of the XVIII Congreso de la Asociación Española de Teledetección, Valladolid, Spain, 24–27 September 2019. [Google Scholar]
  59. Lasanta, T.; Arnáez, J.; Pascual, N.; Ruiz-Flaño, P.; Errea, M.P.; Lana-Renault, N. Space–time process and drivers of land abandonment in Europe. Catena 2017, 149, 810–823. [Google Scholar] [CrossRef]
  60. Expósito, A.; Berbel, J. Agricultural irrigation water use in a closed basin and the impacts on water productivity: The case of the Guadalquivir River Basin (Southern Spain). Water Switz. 2017, 9, 136. [Google Scholar] [CrossRef] [Green Version]
  61. Cleverly, J.; Boulain, N.; Villalobos-Vega, R.; Grant, N.; Faux, R.; Wood, C.; Cook, P.G.; Yu, Q.; Leigh, A.; Eamus, D. Dynamics of component carbon fluxes in a semi-arid Acacia woodland, central Australia. J. Geophys. Res. Biogeosci. 2013, 118, 1168–1185. [Google Scholar] [CrossRef]
  62. Oliver, G.; Florín, M. The wetlands of La Mancha, central Spain: Opportunities and problems concerning restoration. In Bases Ecológicas Para la Restauración de Humedales en la Cuenca Mediterránea; Consejería de Medio Ambiente, Junta de Andalucía: Sevilla, Spain, 1995; ISBN 84-87294-78-2. [Google Scholar]
  63. Zhao, M.; Running, S.W. Drought-Induced Reduction in Global Terrestrial Net Primary Production from 2000 Through 2009. Science 2010, 329, 940–943. [Google Scholar] [CrossRef] [Green Version]
  64. Kim, H.N.; Jin, H.Y.; Kwak, M.J.; Khaine, I.; You, H.N.; Lee, T.Y.; Ahn, T.H.; Woo, S.Y. Why does Quercus suber species decline in Mediterranean areas? J. Asia-Pac. Biodivers. 2017, 10, 337–341. [Google Scholar] [CrossRef]
  65. Peña-Gallardo, M.; Vicente-Serrano, S.M.; Camarero, J.J.; Gazol, A.; Sánchez-Salguero, R.; Domínguez-Castro, F.; El Kenawy, A.; Beguería-Portugés, S.; Gutiérrez, E.; de Luis, M.; et al. Drought sensitiveness on forest growth in peninsular Spain and the Balearic Islands. Forests 2018, 9, 524. [Google Scholar] [CrossRef] [Green Version]
  66. Ruiz-Labourdette, D.; Nogués-Bravo, D.; Ollero, H.S.; Schmitz, M.F.; Pineda, F.D. Forest composition in Mediterranean mountains is projected to shift along the entire elevational gradient under climate change. J. Biogeogr. 2012, 39, 162–176. [Google Scholar] [CrossRef]
  67. del Barrio, G.; Puigdefabregas, J.; Sanjuan, M.E.; Stellmes, M.; Ruiz, A. Assessment and monitoring of land condition in the Iberian Peninsula, 1989–2000. Remote Sens. Environ. 2010, 114, 1817–1832. [Google Scholar] [CrossRef]
  68. Arrogante-Funes, P.; Novillo, C.J.; Romero-Calcerrada, R. Monitoring NDVI inter-annual behavior in mountain areas of mainland Spain (2001–2016). Sustainability 2018, 10, 4363. [Google Scholar] [CrossRef] [Green Version]
  69. Frantz, J.M.; Bugbee, B. Acclimation of plant population to shade: Photosynthesis, respiration, and carbon use efficiency. J. Am. Soc. Hortic. Sci. 2005, 130, 918–927. [Google Scholar] [CrossRef] [Green Version]
  70. Cañizares, M.; Moreno, A.; Sánchez-Ruiz, S.; Gilabert, M.A. Variabilidad de la eficiencia en el uso del carbono a partir de datos MODIS. Rev. Teledetección. 2017, 48, 1–12. [Google Scholar] [CrossRef] [Green Version]
Remotesensing 14 01310 g001 550
Figure 1. Spatial distribution of the average annual precipitation (a) and the average annual T (b) for the 2004–2018 period. Black crosses refer to the location of the sites used for the local analysis (A–J). (c) Vegetation type map [25] used in the estimation of the GPP: GRA (grasslands), SHR (shrublands), CRO (rainfed croplands), ICRO (irrigated crops), MCRO (mosaic cropland), BF (broadleaved forest), NF (needle-leaved forest), MF (mixed forest), SPV (sparse vegetation). (d) Specific forest type map used to derive the NPP: EBF (evergreen broadleaved forest), LDBF (low-altitude 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.
Figure 1. Spatial distribution of the average annual precipitation (a) and the average annual T (b) for the 2004–2018 period. Black crosses refer to the location of the sites used for the local analysis (A–J). (c) Vegetation type map [25] used in the estimation of the GPP: GRA (grasslands), SHR (shrublands), CRO (rainfed croplands), ICRO (irrigated crops), MCRO (mosaic cropland), BF (broadleaved forest), NF (needle-leaved forest), MF (mixed forest), SPV (sparse vegetation). (d) Specific forest type map used to derive the NPP: EBF (evergreen broadleaved forest), LDBF (low-altitude 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.
Remotesensing 14 01310 g001
Remotesensing 14 01310 g002 550
Figure 2. Example of MRA-WT decomposition for the GPP (g m−2 d−1) TS of a deciduous broadleaved pixel. The top panel refers to the original GPP TS. The center panel refers to the variability of the original signal V4–9. The bottom panel shows the interannual A9 component.
Figure 2. Example of MRA-WT decomposition for the GPP (g m−2 d−1) TS of a deciduous broadleaved pixel. The top panel refers to the original GPP TS. The center panel refers to the variability of the original signal V4–9. The bottom panel shows the interannual A9 component.
Remotesensing 14 01310 g002
Remotesensing 14 01310 g003 550
Figure 3. Average annual GPP (a) and NPP (b) for the period 2004–2018.
Figure 3. Average annual GPP (a) and NPP (b) for the period 2004–2018.
Remotesensing 14 01310 g003
Remotesensing 14 01310 g004 550
Figure 4. Annual GPP (a) for all sites and NPP (b) only for forest sites described in Table 1.
Figure 4. Annual GPP (a) for all sites and NPP (b) only for forest sites described in Table 1.
Remotesensing 14 01310 g004
Remotesensing 14 01310 g005a 550Remotesensing 14 01310 g005b 550
Figure 5. GPP (a) and NPP (b) normalized trends (QGPP and QNPP) at a significance level of p-value < 0.1 across the period 2004–2018. The error associated with QGPP (c) and QNPP (d) is also included. 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 5. GPP (a) and NPP (b) normalized trends (QGPP and QNPP) at a significance level of p-value < 0.1 across the period 2004–2018. The error associated with QGPP (c) and QNPP (d) is also included. 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.
Remotesensing 14 01310 g005aRemotesensing 14 01310 g005b
Remotesensing 14 01310 g006 550
Figure 6. SPI slope (QSPI; p-value < 0.3) (a) and temperature slope (QT; p-value < 0.1 (b) derived from daily spatialized meteorological station values.
Figure 6. SPI slope (QSPI; p-value < 0.3) (a) and temperature slope (QT; p-value < 0.1 (b) derived from daily spatialized meteorological station values.
Remotesensing 14 01310 g006
Remotesensing 14 01310 g007 550
Figure 7. Areas with positive and negative GPP trends based on different precipitation and temperature behaviors for natural vegetation (a) and croplands (b).
Figure 7. Areas with positive and negative GPP trends based on different precipitation and temperature behaviors for natural vegetation (a) and croplands (b).
Remotesensing 14 01310 g007
Table
Table 1. Vegetation type, main ecoclimatic features, and coordinates of the ten selected sites with their average annual precipitation (PRE), average annual T, and elevation above sea level (h). SPV (sparse vegetation), GRA (grasslands), SHR (shrublands), CRO (rainfed croplands), ICRO (irrigated crops), EBF (evergreen broadleaved forest), LDBF (low-altitude deciduous broadleaved forest), HDBF (high-altitude deciduous broadleaved forest), LENF (low-altitude evergreen needle-leaved forest), and HENF (high-altitude evergreen needle-leaved forest).
Table 1. Vegetation type, main ecoclimatic features, and coordinates of the ten selected sites with their average annual precipitation (PRE), average annual T, and elevation above sea level (h). SPV (sparse vegetation), GRA (grasslands), SHR (shrublands), CRO (rainfed croplands), ICRO (irrigated crops), EBF (evergreen broadleaved forest), LDBF (low-altitude deciduous broadleaved forest), HDBF (high-altitude deciduous broadleaved forest), LENF (low-altitude evergreen needle-leaved forest), and HENF (high-altitude evergreen needle-leaved forest).
SiteVegetation TypeLatitude (°)Longitude (°)PRE (mm y−1)T (°C)h (m)
ASPV37.07−2.3624117.2463
BGRA38.77−5.5039616.9365
CSHR38.95−0.8746615.9763
DCRO39.34−1.8340414.2756
EICRO41.710.9045014.0224
FEBF36.40−5.5484117.9278
GLDBF42.95−1.61141912.8619
HHDBF43.12−4.78115011.41044
ILENF38.12−2.7663415.7696
JHENF40.74−2.0949711.81099
Table
Table 2. Average annual GPP and NPP. Only pixels with coherent classification between Figure 1c and Figure 1d were used in the calculations of the means and the standard deviations (between brackets). SPV (sparse vegetation), GRA (grasslands), SHR (shrublands), CRO (rainfed croplands), ICRO (irrigated crops), EBF (evergreen broadleaved forest), LDBF (low-altitude deciduous broadleaved forest), HDBF (high-altitude deciduous broadleaved forest), LENF (low-altitude evergreen needle-leaved forest), and HENF (high-altitude evergreen needle-leaved forest).
Table 2. Average annual GPP and NPP. Only pixels with coherent classification between Figure 1c and Figure 1d were used in the calculations of the means and the standard deviations (between brackets). SPV (sparse vegetation), GRA (grasslands), SHR (shrublands), CRO (rainfed croplands), ICRO (irrigated crops), EBF (evergreen broadleaved forest), LDBF (low-altitude deciduous broadleaved forest), HDBF (high-altitude deciduous broadleaved forest), LENF (low-altitude evergreen needle-leaved forest), and HENF (high-altitude evergreen needle-leaved forest).
Vegetation TypeGPP (kg m–2 y–1)NPP (kg m–2 y–1)
SPV0.6 (0.3)
GRA0.8 (0.3)
SHR0.7 (0.2)
CRO0.6 (0.3)
ICRO1.4 (0.4)
EBF1.6 (0.5)0.4 (0.3)
LDBF2.0 (0.3)0.8 (0.2)
HDBF1.8 (0.4)0.5 (0.2)
LENF1.2 (0.3)0.20 (0.16)
HENF1.1 (0.3)0.20 (0.14)
Table
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).
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).
SiteGPP Q G P P 10 − 2   rSPI-GPPrGPP-TNPP Q N P P 10 − 2   rNPP-SPIrNPP-T
A0.10−0.93 ± 0.090.590.12−−−−
B0.530.26 ± 0.050.72−0.52−−−−
C0.65−0.056 ± 0.0140.74−0.56−−−−
D0.33−0.17 ± 0.030.63−0.64−−−−
E1.600.10 ± 0.020.180.27−−−−
F2.22−0.202 ± 0.0110.56−0.430.53−0.280 ± 0.0150.61−0.53
G2.110.15 ± 0.020.550.650.910.173 ± 0.0170.520.64
H2.190.161 ± 0.0060.160.420.770.111 ± 0.0050.170.29
I1.08−0.021 ± 0.0150.64−0.530.17−0.23 ± 0.040.75−0.69
J0.93−0.029±0.0110.37−0.210.13−0.039 ± 0.0130.27−0.30
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Martínez, B.; Sánchez-Ruiz, S.; Campos-Taberner, M.; García-Haro, F.J.; Gilabert, M.A. Exploring Ecosystem Functioning in Spain with Gross and Net Primary Production Time Series. Remote Sens. 2022, 14, 1310. https://doi.org/10.3390/rs14061310

AMA Style

Martínez B, Sánchez-Ruiz S, Campos-Taberner M, García-Haro FJ, Gilabert MA. Exploring Ecosystem Functioning in Spain with Gross and Net Primary Production Time Series. Remote Sensing. 2022; 14(6):1310. https://doi.org/10.3390/rs14061310

Chicago/Turabian Style

Martínez, Beatriz, Sergio Sánchez-Ruiz, Manuel Campos-Taberner, F. Javier García-Haro, and M. Amparo Gilabert. 2022. "Exploring Ecosystem Functioning in Spain with Gross and Net Primary Production Time Series" Remote Sensing 14, no. 6: 1310. https://doi.org/10.3390/rs14061310

APA Style

Martínez, B., Sánchez-Ruiz, S., Campos-Taberner, M., García-Haro, F. J., & Gilabert, M. A. (2022). Exploring Ecosystem Functioning in Spain with Gross and Net Primary Production Time Series. Remote Sensing, 14(6), 1310. https://doi.org/10.3390/rs14061310

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop