Next Article in Journal
Radar Interferometry: 20 Years of Development in Time Series Techniques and Future Perspectives
Previous Article in Journal
Democratic Republic of the Congo Tropical Forest Canopy Height and Aboveground Biomass Estimation with Landsat-8 Operational Land Imager (OLI) and Airborne LiDAR Data: The Effect of Seasonal Landsat Image Selection
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Remote Sensing and Bio-Geochemical Modeling of Forest Carbon Storage in Spain

by
Sergio Sánchez-Ruiz
1,*,
Fabio Maselli
2,
Marta Chiesi
2,
Luca Fibbi
2,
Beatriz Martínez
1,
Manuel Campos-Taberner
1,
Francisco Javier García-Haro
1 and
María Amparo Gilabert
1
1
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
2
Istituto per la BioEconomia, Consiglio Nazionale delle Ricerche, 50019 Sesto Fiorentino, Italy
*
Author to whom correspondence should be addressed.
Remote Sens. 2020, 12(9), 1356; https://doi.org/10.3390/rs12091356
Submission received: 20 March 2020 / Revised: 22 April 2020 / Accepted: 23 April 2020 / Published: 25 April 2020

Abstract

:
This study simulates annual net primary production (NPP) of forests over peninsular Spain during the years 2005–2012. The modeling strategy consists of a linked production efficiency model based on the Monteith approach and the bio-geochemical model Biome-BGC. Recently produced databases and data layers over the study area including meteorological daily series, ecophysiological parameters, and maps containing information about forest type, rooting depth, and growing stock volume (GSV), were employed. The models, which simulate forest processes assuming equilibrium conditions, were previously optimized for the study area. The production efficiency model was used to estimate daily gross primary production (GPP), while Biome-BGC was used to simulate growth (RG) and maintenance (RM) respirations. To account for actual forest conditions, GPP, RG, and RM were corrected using the ratio of the remotely-sensed derived actual to potential GSV as an indicator of the actual state of forests. The obtained results were evaluated against current annual increment observations from the Third Spanish Forest Inventory. Coefficients of determination ranged from 0.46 to 0.74 depending on the forest type. A simplified dataset was produced by applying regular increments in air temperature and reductions in precipitation to the original 2005–2012 daily series with the goal of covering the range of variation of the climate projections corresponding to the different climate change scenarios reported in the literature. The modified meteorological series were used to simulate new GPP, RG, and RM through Biome-BGC and corrected using GSV. Precipitation was confirmed as the main limiting factor in the study area. In the regions where precipitation was already a limiting factor during 2005–2012, both the increment in air temperature and the reduction in precipitation contributed to a reduction of NPP. In the regions where precipitation was not a limiting factor during 2005–2012, the increment in air temperature led to an increment of NPP. This study is therefore relevant to characterize the growth of Spanish forests both in current and expected climate conditions.

Graphical Abstract

1. Introduction

The signature of the Kyoto Protocol has renewed interest in the study of forests. Forests act as potential carbon sinks for the mitigation of climate change impacts, but in turn, they are also affected by climate. Specifically, Mediterranean forests are very sensitive to climate change effects. According to [1], an increment of annual mean air temperature of 3–4 °C and a reduction of annual precipitation of 20% would lead to a decrease in photosynthesis and a decline of biomass accumulation during the current century. In this context, the net primary production (NPP) of terrestrial ecosystems is a variable that can reflect these changes, and its monitoring becomes key to apply the needed adaptation and mitigation actions. NPP is a function of gross primary production (GPP), the flux of carbon fixed by plants through photosynthesis, minus plant autotrophic respiration (RA), the flux of carbon released back to atmosphere due to internal plant metabolism. Therefore, NPP quantifies the flux of carbon stored by plants in their structure.
Three main methods used for the estimation of vegetation carbon fluxes and storage are (1) field measurements, (2) remote sensing, and (3) bio-geochemical modeling, which are all briefly described below.
(1) Field measurements are the most accurate and are usually considered as the reference. The eddy covariance (EC) technique allows the estimation of net ecosystem carbon exchange (defined as NPP minus heterotrophic respiration) and total respiration from spectral and micrometeorological continuous measurements [2]. However, it is a site-specific method that only accounts for several hundreds of squared meters around the measurement site. For its best functioning, the required local conditions include flat upwind terrain and lack of horizontal air mass movements. Thus, it is not particularly useful for the study of spatial patterns. Forest inventories (both at regional and national scales) provide periodic estimates of forest attributes such as volume and biomass [3], from which carbon storage can be estimated and categorized by forest type, species, and for different administrative units (e.g., province, region, and state). This methodology can be applied at large scales, but it is expensive to implement in terms of invested resources and time [4]. For example, in the particular case of the study area (peninsular Spain, see Section 3), only one national inventory is performed every 10 years [5].
(2) Remote sensing data derived from measurements of radiation reflected or emitted by the Earth’s surface, which are available over extensive areas, allow carbon fluxes and storage to be estimated [6,7]. Several efforts to estimate GPP have recently been made in the study area [8,9,10] using the Monteith approach [11], but the estimation of NPP is more problematic due to the difficulty in the estimation of RA [12,13]. Optical data for the estimation of carbon storage are available at a suitable frequency and globally, but they suffer from saturation of the signal and the conditions of the atmosphere [4]. Radar data are not affected by atmospheric conditions but also suffer from signal saturation, and the complex processing needed for their use is problematic in topographically rugged areas [14].
(3) Bio-geochemical models such as Biome-BGC [15,16] can simulate vegetation processes (including carbon fluxes and storage) and have been recently applied in Mediterranean regions with success [17,18,19]. These models, however, require a large amount of input data, which can be classified into ecophysiological parameters and drivers (site physical data and meteorological time series). Moreover, Biome-BGC assumes that the considered ecosystem is in equilibrium with its surroundings, and its outputs must be corrected to account for actual ecosystem conditions. Such correction can be carried out by the use of growing stock volume (GSV) [20], which is the volume over bark of all living trees with a diameter at breast height (DBH) ≥ 75 mm [21].
Previous studies have produced the data needed to use bio-geochemical models in the study area despite the difficulty of obtaining them at a spatially distributed level. The methodology for the production of the principal daily meteorological series (minimum and maximum air temperature, precipitation, and incoming solar radiation) by means of ordinary kriging and artificial neural networks (ANN) was developed in [22,23]. Uncertainties of 0.93 °C, 1.94 mm, and 3.16 MJ m-2 d-1 were respectively obtained for air temperature, precipitation, and incoming solar radiation. Other required meteorological daily series (day length, daylight air temperature, humidity) can be obtained through the use of microclimate simulation models, such as MT-CLIM [24]. Ecophysiological parameters for the main vegetation types in Mediterranean areas were set by [17] taking into account the typical water stress effect. Soil texture maps (with standard deviation generally below 4%, although it can rise to ~100% in mountainous terrain) [25], land cover maps (with geometrical error ≤ 3 m) [26], and digital elevation models [27] are also available for the study area. A rooting depth map was obtained as a result of a calibration of Biome-BGC in the study area [19]. The rooting depth is the depth at which plants are able to grow roots [15] and is a crucial input in Biome-BGC for the simulation of the site water balance. A GSV map with RMSE = 64 m3 ha-1 was produced for the study area from a combination of remote sensing and forest inventory data [28].
Following the previous research efforts described above, and based on the capacity of NPP to track the changes in vegetation due to the climate change, the objectives of the present study were: (1) to apply a modeling strategy to assess forest annual NPP at the peninsular Spain scale using remote sensing and ancillary data and (2) to simulate forest annual NPP in future climate conditions (increases in air temperature and reductions in precipitation). The modeling strategy consisted of the combination of remote sensing and bio-geochemical modeling (Monteith approach and Biome-BGC) driven by information layers currently available over peninsular Spain and mostly derived from remote sensing data. It was applied to the five main forest ecosystems in the study area and validated against forest inventory observations. The simulation of future NPP was based on Biome-BGC simulations using meteorological time series modified on the basis of future climate scenarios.

2. Modeling Strategy

To assess forest NPP, a modeling strategy (schematically summarized in Figure 1) that combines two models was applied (for the description of the data involved in the calculation of the different variables, see Section 3.3). The first model was a Monteith-like production efficiency model optimized for the study area [8], which was used to calculate annual forest GPP at 1-km spatial resolution from daily estimates. The model was driven by photosynthetically active radiation (PAR), the fraction of PAR absorbed by the vegetation canopy (fAPAR), precipitation, and air temperature:
GPP = ε max i = 1 N PAR i f APAR , i ε W , i ε T , i
where the subindex i refers to the day of year, while N is the number of days in the year. The maximum light use efficiency εmax was set as 1.2 g MJ−1. PAR was calculated as 46% of incident global solar radiation [29]. fAPAR was computed according to the Roujean and Bréon algorithm [30]. Red (ρRED) and near-infrared (ρNIR) reflectances were first calculated from bidirectional reflectance distribution function (BRDF) parameters k0, k1, and k2 for an optimal angular geometry in the solar principal plane. fAPAR was calculated as a linear function of the renormalized difference vegetation index
RDVI = (ρNIRρRED)/(ρNIR + ρRED)1/2
εW is the water stress coefficient (CWS) as used in [31] that accounts for short-term water stress
CWS = (PRE/PET + 1)/2
where precipitation (PRE) and potential evapotranspiration (PET, calculated from air temperature and incoming global solar radiation through the Jensen-Haise equation [32]) are accumulated for 60 days. εT is the TMIN_scalar used in the MOD17 algorithm [33] to account for thermal stress caused by low temperatures: a linear ramp function of daily minimum air temperature with thresholds depending on the vegetation type.
The second model was Biome-BGC [15,16], a bio-geochemical model that, based on daily meteorological data, site information (i.e., soil, vegetation, and site conditions) and parameters describing the ecophysiological features of the vegetation, is able to simulate fluxes and storage of water, carbon, and nitrogen at different spatial scales in terrestrial ecosystems. In particular, Biome-BGC identifies a quasi-equilibrium condition (the so-called climax) with the local ecoclimatic situation, thus quantifying the initial amounts of all carbon and nitrogen pools. Version 4.2 is available at the Numerical Terradynamic Simulation Group (NTSG) website (http://www.ntsg.umt.edu/project/biome-bgc.php) and was used together with a set of ecophysiological parameters calibrated for the considered forest types [17]. More information on the use and setup of Biome-BGC in the study area can be found in [19].
NPP is the difference between GPP and RA, with RA equal to the sum of growth (RG) and maintenance (RM) respirations. Since Biome-BGC assumes that the considered ecosystem is in equilibrium with its surroundings, its outputs must be corrected to account for actual ecosystem conditions. Note that Equation (1) estimates the GPP of all vegetation growing in each pixel. Therefore, it also must be corrected to account only for forest GPP. The correction applied in the present study was based on the distance to climax of the ecosystems concept [34], from which it is deduced that ratios of actual to potential tree cover or GSV are representative of the ecosystem development status or, in other words, its distance to climax [20]. The correction to the simulated carbon fluxes was carried out following the methodology explained in [20]:
NPP = f (GPP − RG) − v RM
where GPP is estimated by Equation (1), while f and v are, respectively, the ratios of actual to potential tree cover and GSV and are calculated as follows:
f = (1 − ev LAI)/(1 − e−LAI)
and
v = 50 GSV ρ C/σ
where LAI (m3 m−3) and σ (g m−2) are annual maximum leaf area index and dead stem carbon simulated by Biome-BGC, ρ is the basic wood density, C is a biomass expansion factor, and the 50 factor accounts for the transformation from carbon mass to dry mass (2 kg kg−1) and for unit conversion to m3 ha−1.

3. Study Area and Data

3.1. Peninsular Spain

Peninsular Spain, located between −10° and 3° longitude and between 36° and 40° latitude, is a heterogeneous and large region of Southwestern Europe. The climate ranges from Atlantic to semiarid following a NW−SE gradient. Annual precipitation ranges from more than 2000 mm to less than 200 mm following the same gradient. The elevation ranges from sea level to 3479 m. These features favor the existence of a great range of different ecosystems [35,36].
Human activity has also had an impact in the study area [37], resulting in a great heterogeneity in Spanish forest types and characteristics. Of the total area, 55% is covered by wooded lands, while 37% is covered by tree species [38]. Spanish forests can be grouped into five main types, i.e., evergreen broadleaved forest (EBF), low-altitude deciduous broadleaved forest (LDBF), high-altitude deciduous broadleaved forest (HDBF), low-altitude evergreen needleleaved forest (LENF), and high-altitude evergreen needleleaved forest (HENF) [19]. Quercus ilex is the most widespread broadleaved tree species (total broadleaves occupy 46% of the tree-covered area), while Pinus halepensis is the most widespread needleleaved one (total needleleaves occupy 34% of the tree-covered area).

3.2. Reference Forest Observations

Spanish forests’ field data were derived from the third Spanish Forest Inventory (NFI3), performed during 1997–2007 [38]. A total of 92,340 plots distributed in a regular 1 km × 1 km grid were surveyed in the study area. Each plot consists of four concentric circular subplots with 5 m, 10 m, 15 m, and 25 m radii. Trees with 75 mm ≤ DBH < 125 mm, 125 mm ≤ DBH < 225 mm, 225 mm ≤ DBH < 425 mm, and DBH ≥ 425 mm were measured respectively in each subplot at a breast height of 130 cm. For each of the five forest types considered in the present study, a subset of NFI3 plots was selected considering the spatial and species homogeneity and the availability of site information. Specifically, only NFI3 plots with Quercus ilex, Quercus robur, Fagus sylvatica, Pinus pinea, and Pinus nigra were selected for EBF, LDBF, HDBF, LENF, and HENF, respectively. In addition, it was imposed that the 9 pixels in a 3 × 3-pixel window of a forest type map (see Section 3.3) centered in each plot presented the same forest type (except for LDBF—at least 3 of 9—and LENF—at least 8 of 9—due to plot availability). For each plot, GSV and current annual increment (CAI) observations of the four subplots were added to obtain values representative of the whole plot. Table 1 summarizes the data used per forest type.

3.3. Data Used for the NPP Modeling

Most datasets used for the current NPP modeling were produced during previous research exercises. The SIOSE (Land Cover and Use Information System of Spain) land cover database was provided by ©Instituto Geográfico Nacional through http://www.siose.es/web/guest/descargar. It consists of a cartographic scale 1:25000 polygon database generated by photointerpreting and digitalizing reference data, which included SPOT5 and Landsat-5 TM images together with local ground ancillary data [26,39]. The land cover data contains the percentage of surface occupied by each of the present classes in each polygon. With the help of the global 3 arc second digital elevation model from the Shuttle Radar Topography Mission [27] (downloaded from http://edcftp.cr.usgs.gov/pub/data), the original legend of the land cover data was regrouped into the five forest types mentioned in Table 1 using an 800-m threshold to separate low- and high-altitude classes [19]. The resulting forest type map is shown in Figure 2a.
Daily ground measurements of precipitation, and minimum and maximum air temperature were provided by the Spanish Meteorological Agency (www.aemet.es) from ~400 meteorological stations distributed throughout the study area. Ordinary kriging was used to interpolate these data and obtain daily 1-km spatial resolution maps during the 2005–2012 period. These maps were used as inputs in MT-CLIM [24] to simulate daylength, daylight average partial pressure of water vapor, and daily average air temperature. Version 4.3 of the code was used and is available at the NTSG website (http://www.ntsg.umt.edu/project/mt-clim.php).
SEVIRI products LSA-201 (downward surface shortwave flux estimated every 30 min, MDSSF) and LSA-203 (daily downward surface shortwave flux, DIDSSF) [40] for the period 2007–2012 were acquired from the LSA-SAF server (http://landsaf.meteo.pt). The images were reprojected to a regular longitude/latitude 1-km spatial resolution grid. MDSSF was used to calculate DIDSSF when it was unavailable. Incoming global solar radiation was obtained by applying ANN to the abovementioned air temperature and precipitation maps [22] and was used to fill DIDSSF gaps and to calculate DIDSSF for the years 2005 and 2006. To apply the gap-filling procedure, a relationship was used between the two datasets, which were previously validated [22,23].
The main ecoclimatic features of the five Spanish forest types derived from the datasets exposed above are summarized in Table 2.
MODIS products MCD43A1 and MCD43A2 [41] for the period 2005–2012 were retrieved from the online Reverb, courtesy of the NASA EOSDIS Land Processes Distributed Active Archive Center, USGS/Earth Resources Observation and Science Center, Sioux Falls, South Dakota (reverb.echo.nasa.gov). MCD43A1 is a 500-m spatial resolution 8-day composite containing the BRDF parameters k0, k1, and k2. fAPAR was computed through the methodology reported in [30] from the BRDF parameters contained in MCD43A1 images (see Section 2), which were previously reprojected to a regular latitude/longitude 1-km spatial resolution grid. A noise-reduction and gap-filling method dependent on the quality of the measurements [42] was applied, and the 8-day fAPAR was linearly interpolated to obtain daily images.
Soil texture (clay, sand, and silt content) maps [25] were downloaded from the European Soil Database (http://esdac.jrc.ec.europa.eu/resource-type/european-soil-database-soil-properties) and reprojected to regular longitude/latitude 1-km spatial resolution grids comprising the study area.
A 1-km spatial resolution rooting depth map was obtained from a calibration of Biome-BGC in the study area. Daily GPP was first calculated using a production efficiency model optimized for the study area [8] in four sites where EC towers were placed. Several daily GPP series were also simulated using Biome-BGC for the same sites and varying the rooting depth. Both series were compared, and the optimum rooting depth was chosen as the one that resulted in the minimum root mean square error between them. After ascertaining that this methodology improved the accuracy of GPP simulated using Biome-BGC by comparing to EC-derived GPP, the methodology was extended to the whole study area and the rooting depth map was obtained. Detailed methods on its production can be found in [19]. Note that rooting depth used by Biome-BGC does not indicate depth to bedrock. Instead, it refers to an effective soil depth defined as the depth at which plants are able to grow roots [15].
A 30-m spatial resolution GSV map was obtained from the combination of around 50000 NFI3 plot-level GSV estimations and a multitemporal Landsat dataset (~8000 Landsat-5 TM and Landsat-7 ETM+ scenes covering the study area and the NFI3 period, 1997–2007). A total of 805 predictors were calculated from Landsat spectral reflectances including time metrics, texture metrics, and vegetation indices. The most important predictors were identified by means of a guided regularized random forests algorithm (e.g., [43,44]) and used to generate the GSV map. Google Earth Engine [45] was used for data management. Detailed methods on the production of this map can be found in [28]. The 30-m spatial resolution GSV map was aggregated to 1-km spatial resolution (Figure 2c) for coherence with the other datasets used in the present study.

4. Data Processing

4.1. Simulation of Forest NPP in Current Condition

The NPP modeling strategy described in Section 2 was applied at a spatial resolution of 1 km and a daily temporal scale for the study period (2005–2012) using the inputs described in Section 3. In Equation (4), two GPPs were used: the one estimated by the production efficiency model (GPPPEM, Figure 2b) and the one simulated by Biome-BGC (GPPBGC) [19], which also simulated RG and RM. In Equation (6), the GSV map and the data reported in Table 3 [46] were used.
The NPP estimated by the two methods was assessed for the pixels containing the NFI3 plots indicated in Table 1. For each of these pixels, a reference annual NPP was calculated from NFI3 CAI as
NPPNFI3 = 50 CAI ρ C/D
D being the stem carbon allocation ratio (Table 3), and the 50 factor accounting for both the unit conversion to m3 ha−1 and the transformation from carbon mass to dry mass (2 kg kg−1). Annual NPP was calculated for the same locations following the methodology reported in Section 2, but GSV from NFI3 was used instead of the GSV map to test the modeling strategy in the best possible conditions. The calculated NPP was compared to NPPNFI3 by means of the following statistics calculated by forest type: coefficient of determination (R2), mean bias error (MBE), and root mean square error (RMSE).

4.2. Simulation of Forest NPP in Future Condition

The simulation of GPP, RG, and RM through the use of Biome-BGC (as indicated in Section 4.1) was driven by different daily precipitation and minimum and maximum air temperature datasets, which led to different daylight average partial pressure of water vapor, and daily average air temperature through the use of MT-CLIM.
A large dataset of future climate projections over the study area is available through https://www.adaptecca.es. It contains, among others, projections of daily maximum and minimum air temperature and precipitation for different climate change scenarios (RCP 4.5 and RCP 8.5) at both spatial and station levels. However, they are obtained from different models resulting in big differences among the produced datasets. Moreover, the spatial resolution of the spatialized projections (11 km) is too coarse compared to the one used in the present study (1 km), while projections at station-level are restricted to the locations where meteorological stations from the Spanish Meteorological Agency are placed. Furthermore, climatic projections generally present great uncertainties [47], especially for precipitation. Therefore, to address the range of variation of the climate projections corresponding to the different climate change scenarios (e.g., no mitigation, soft mitigation, and hard mitigation) that can be found in the literature, a simplified dataset of possible future changes in meteorological variables (i.e., daily precipitation and maximum and minimum air temperature) was generated for use in the present study.
The original 2005–2012 series were taken as reference and the same variation was applied to each value in each series. In the case of air temperature, absolute increments (ΔT) of 1 °C, 2 °C, 3 °C, and 4 °C were applied to both minimum and maximum series. In the case of precipitation, relative variations were applied in the form of reductions (δPRE) of 10%, 20%, 30%, and 40% of its reference value. These variations were chosen to take into account climatic projections over the study area from both global and regional climatic models for different CO2 emission scenarios [47]. The four modified air temperature series, the four modified precipitation series, and the reference series were combined to obtain 24 new 8-year sets of simulated daily GPP, RG, and RM through Biome-BGC. Finally, the corresponding NPP series were calculated for each set using Equation (4) fed with the same GSV as in the current condition, i.e., from the GSV map.

5. Results

5.1. NPP in Current Condition

The assessment of estimated NPP per forest type is illustrated by the scatterplots in Figure 3. The largest NPPNFI3 explained variance amongst forest types was obtained for EBF with almost 70% (R2 = 0.69); ~60% (R2 = 0.60, R2 = 0.62, and R2 = 0.59 respectively) was obtained for LDBF, HDBF, and LENF; and slightly less than 50% (R2 = 0.46) for HENF. Small scattering was observed for EBF except for a great NPPNFI3 plot, which was heavily underestimated. A general overestimation, more notable for great NPP values, was also observed. LDBF and HDBF exhibited similar behavior: estimates generally followed the 1:1 line until around 300 g m−2 a−1, where scattering appeared (more for LDBF), and they saturated at around 700 g m−2 a−1. LENF presented a constant underestimation, but with small scattering. HENF exhibited more scattering than LENF, but estimates followed the 1:1 line until around 200 g m−2 a−1, where saturation appeared.
The accuracy statistics resulting from the NPP assessment are reported in Table 4. The greatest difference between the two NPPEST to explain the NPPNFI3 variance was found for LENF with 10 percentage points (pp) (from R2 = 0.59 with GPPPEM to R2 = 0.49 with GPPBGC). Around half this difference was found for EBF, LDBF, and HENF (respectively from R2 = 0.69, R2 = 0.60, and R2 = 0.46 with GPPPEM to R2 = 0.74, R2 = 0.56, and R2 = 0.51 with GPPBGC), while only 2 pp were found for HDBF (from R2 = 0.62 with GPPPEM to R2 = 0.60 with GPPBGC). More of the NPPNFI3 variance was explained for broadleaved forests than for needleleaved forests independently of the used GPP. According to the calculated MBE, NPP was overestimated in broadleaved forests but underestimated in needleleaved forests. In addition, the MBE obtained for HDBF and LENF was an order of magnitude greater than for the rest of the forest types. RMSE was high for all forest types compared to their corresponding mean NPPNFI3 values (114 g m−2 a−1, 418 g m−2 a−1, 421 g m−2 a−1, 287 g m−2 a−1, 209 g m−2 a−1 respectively for EBF, LDBF, HDBF, LENF, and HENF).
Figure 4a,b show, respectively, the average annual forest NPP calculated using Equation (4) with GPPPEM and with GPPBGC. The spatial pattern was the same in both cases: NPP was larger in the NW and decreases towards the SE, except in the most southern and the most eastern (this one located in the north too) regions, where NPP was also great. NPP calculated with GPPBGC was generally smaller than NPP calculated with GPPPEM. Table 5 summarizes the basic statistics of the mentioned GPPs and their corresponding NPPs per forest type. GPPPEM maxima are greater than GPPBGC maxima for broadleaved forest, but the opposite relationship occurred for needleleaved forests. GPPPEM means were greater than GPPBGC ones for all forest types. NPP maxima and means exhibited a similar behavior. However, these differences do not seem very relevant, especially the differences in the mean values.

5.2. NPP in Future Condition

Simulated annual forest NPP in possible future scenarios of climate change presents a spatial pattern similar to that exhibited during 2005–2012 (Figure 4a,b). This is shown in Figure 4c, where the average annual forest NPP calculated using Equation (4) with GPPBGC, ΔT = +2 °C, and δPRE = –20% for the corresponding 8-year period is represented. However, it can be appreciated how the regions that already presented great (or small) NPP in Figure 4a,b present even greater (or smaller) NPP in Figure 4c. The basic statistics of Figure 4c are summarized in Table 6. A comparison of Table 6 with NPP with GPPBGC in Table 5 illustrates that minima are all practically the same; maxima are greater in Table 6 for all forest types; means are greater in Table 6 except for HENF, for which it is smaller; and standard deviations are greater in Table 6 for EBF, LDBF, and LENF, while they are nearly equal for HDBF and HENF.
Figure 5 shows the variation of annual forest NPP with air temperature increments (ΔT) and precipitation reductions (δPRE) for a single pixel of each forest type. The main eco-climatic features of the selected pixels together with their geographic coordinates are summarized in Table 7. The pixels were selected for the sake of homogeneity and representativeness: at least a 5 × 5-pixel window (7 × 7 for HENF) with the same forest type surrounding each pixel was imposed, and they are also representative of the main eco-climatic features of each forest type (Table 1). In the case of EBF, the surface is nearly flat and only slight decreases in NPP due to increases in air temperature were found. NPP in LDBF was strongly reduced when both an increment of air temperature of at least 2 °C and a reduction of precipitation of at least 30% were combined, but a softer reduction in NPP was also observed when precipitation was reduced by 40% independently of the variation in air temperature. NPP in HDBF did not present any significant change due to increments in air temperature nor to reductions in precipitation. NPP in LENF was slightly but consistently decreased as precipitation was reduced. HENF presented a dependence of NPP with precipitation too, but the reduction in NPP due to reductions in precipitation was stronger than in the case of LENF.

6. Discussion

6.1. NPP in Current Condition

The present study focuses on the use of a modeling strategy for the simulation of annual forest NPP over peninsular Spain through the combination of a production efficiency model, which was used for the estimation of annual forest GPP, and the bio-geochemical model Biome-BGC, which was used to simulate GPP, RG, and RM. Both models were previously optimized for their application in the study area [8,19] and were driven by a combination of remote sensing observations and ground meteorological data. The modeling strategy described in [20] was developed to overcome one of Biome-BGC’s limitations, i.e., the assumption of the considered ecosystem to be in equilibrium with its surroundings (climax conditions), through the correction of Biome-BGC’s outputs to account for actual conditions of the considered ecosystem and through the use of a production efficiency model for the estimation of GPP. However, other limitations of Biome-BGC and the production efficiency model must be taken into account to evaluate the obtained results. The space where the models are applied is divided into cells and the simulations take place individually in each cell without interaction with other cells during the functioning of the models. Therefore, the spatial relationships and variations observed in the outputs of the models are provided by the ones already present in their inputs. This is important when interpreting the spatial patterns presented in Figure 4. Furthermore, both models assign a single ecosystem functional type (in the case of the present study, a forest type) to each cell and assume that it is homogeneous independently of the spatial resolution at which the models work and that it does not change through time. This affects light use efficiency in the production efficiency model and the ecophysiological parameters in Biome-BGC, but it is also relevant to interpret the results of NPP simulation in future condition because the possibility of the current forest type being replaced by another ecosystem is not contemplated. Only the current forest type reactions to possible variations in air temperature and precipitation due to climate change were therefore analyzed.
The two different GPP sources, i.e., the one estimated through the production efficiency model (GPPPEM) and the one simulated using Biome-BGC (GPPBGC), that are involved in the calculation of NPP were validated against daily EC-derived GPP from different ecosystems (forest included) in a previous study [19]. Coefficients of correlation between 0.64 and 0.86 and between 0.52 and 0.75 were respectively obtained for GPPPEM and GPPBGC, while RMSE ranged from 0.7 g m−2 d−1 to 1.7 g m−2 d−1 and from 0.9 g m−2 d−1 to 1.8 g m−2 d−1. In the present study, new datasets (e.g., the respirations simulated using Biome-BGC and the GSV map) with their own uncertainty were incorporated into the calculation of NPP, which was assessed through an indirect validation against a reference NPP (NPPNFI3) calculated from CAI observations derived from forest inventory measurements.
The strength of this approach is to combine the directly derived GPPPEM with all respirations and allocations simulated by the bio-geochemical model Biome-BGC. However, actual NPP is difficult to predict due to the numerous factors and human-induced disturbances affecting it [48]. In the current case, i.e., in all applications over large areas, the modeling approach must be based on easily collectable and available information otherwise its applicability is missed. Additionally, the accuracy of the NPP estimates is also strictly dependent on the quality of the drivers utilized, among which GSV is one of the most important [49].
GSV from NFI3 was currently used instead of the GSV map to evaluate the modeling strategy in the best possible conditions. As expected, greater relative errors were obtained for NPP than for GPP. Some differences in the results obtained among the five forest types can be appreciated. The greatest R2 and smallest errors were obtained for EBF, while the smallest R2 and greatest errors were obtained for HENF. This is in agreement with the results obtained in Italy by [49], who tested a similar NPP modeling approach against reference CAI data collected by the Italian National Forest Inventory. They also found a clear NPP underestimation for Mediterranean pines, which can be partly attributed to a suboptimal parameterization of the coefficients reported in Table 3. However, there is no conclusive indication that one GPP should be preferred over the other. Therefore, the use of GPPBGC in the calculation of NPP in future condition is justified.

6.2. NPP in Future Condition

Some tests were firstly performed using the scenarios from AdapteCCa (https://www.adaptecca.es) data as inputs in Biome-BGC for some sites where station-level AdapteCCa data obtained with the same model were available. Daily GPP and respirations were simulated for 2006–2100. However, unexpected extreme (i.e., both very high and very low) carbon fluxes values were obtained. This, together with the abundance of different models that provide different climate projections and their high uncertainty, led to the decision of generating a simplified dataset of possible future changes in meteorological variables trying to cover most of the climatic projections (both from different models and for different scenarios). This simplified dataset was elaborated by applying variations to a reference dataset that was previously used in the simulation of GPP by both the production efficiency model and Biome-BGC [19], the daily precipitation and maximum and minimum air temperature for the 2005–2012 period, so different sets of 8-year periods were produced instead of a long time series as in the case of AdapteCCa. Thus, the results presented in Section 5.2 might approximately correspond to different 8-year periods over the present century depending on whether any actions are taken to mitigate climate change and how intense their application is.
Still, the generation of the abovementioned 8-year datasets presents some important limitations. First, there are some limitations regarding the direction of the modifications applied to the reference period. Both global and regional climatic models predict increment in air temperature in the study area for the different emission scenarios. Therefore, only increments in temperature were considered. Whereas in the case of precipitation, an increase is also expected globally, but regional models predict reductions in annual precipitation over the study area. Therefore, only reductions in precipitation were considered. Second, there are limitations regarding the temporal variation of the applied modifications. While climatic models usually provide, for example, different variations for winter and summer within the same year, the datasets generated in the present study were produced by applying the same modification to each value in the series, that is, a regular variation. Regarding the spatial variation, which is especially relevant in the case of precipitation, it is provided by the spatial variation of the reference series itself.
Other limitations of the current approach are due to the use of a bio-geochemical model, which can be suboptimal for the simulation of future climatic conditions. Biome-BGC respirations are actually simulated using a Q10 function of air temperature [50], which however does not account for the possible response and adaptation of vegetation to high air temperature [51]. The consideration of stable GSV for driving the NPP modeling strategy in future condition is also disputable but can be justified by the unpredictable dependence of this forest attribute on human activities (tree cutting and thinning operations, wildfires of arson origin, etc.).
When comparing the annual forest NPP shown in Figure 4b (reference period 2005–2012) and Figure 4c (ΔT = +2 °C and δPRE = –20%), it can be appreciated how NPP is greater in Figure 4c in the regions where it was already great in Figure 4b and smaller in Figure 4c in the regions where it was already small in Figure 4b, that is, it presents a more extreme variation in Figure 4c. This is supported by the statistics reported in Table 5 and Table 6: maxima and standard deviations of NPP are greater in Table 6 (corresponding to Figure 4c, ΔT = +2 °C, and δPRE = –20%). The increment of NPP in already very productive regions (mostly the north of the study area) may be caused by the increment in air temperature and the water not being a limiting factor. On the other hand, the reduction of NPP in already scant productive regions (mainly the SE of the study area) may be caused by the reduction in precipitation and water already being a limiting factor during 2005–2012. From Figure 5 it can be appreciated how, although no significant changes are found in annual NPP for EBF nor HDBF, precipitation is the main factor responsible for NPP reductions in the rest of the forest types (LDBF, LENF, and HENF).

7. Conclusions

The current research was aimed at investigating the NPP dynamics of Spanish forests both in the current and expected climate conditions. The analysis was conducted using an integrated modeling approach applied at high spatial and temporal resolutions (i.e., 1 km2 and daily) in order to grasp the most relevant variations that characterize this highly heterogeneous Mediterranean country.
The main conclusions that can be drawn from the experiment are:
  • When applied in current climate conditions, the integrated modeling approach is capable of reproducing the general NPP variability of broadleaved forests in Spain (R2 > 0.60) but underestimates the NPP of needleleaf forests.
  • When applied in future climate conditions, the approach produces results that can be explained considering the ecoclimatic features of the present forest types. NPP increases are predicted for more temperate-humid forests, while decreases are predicted for forest already subject to water limitation.
All these results contribute to the understanding of the behavior and the risks of Mediterranean ecosystems in the context of expected climate variability. However, they are affected by the limitations in model functioning and the datasets used. Therefore, future research will focus on improving the model capacity to simulate NPP in warmer and drier climate conditions and identifying more plausible future scenarios.

Author Contributions

Conceptualization, methodology, software, validation, formal analysis, investigation, resources, data curation, writing-original draft preparation, writing-review and editing, visualization, supervision, project administration, and funding acquisition: S.S.-R., F.M., M.C., L.F., B.M., 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 study was partially funded by the ESCENARIOS project from the Spanish Ministry of Science, Innovation and Universities (CGL2016-75239-R) co-funded by FEDER, and by the LSA SAF CDOP-3 project from the European Organization for the Exploitation of Meteorological Satellites (EUMETSAT).

Acknowledgments

Part of the data was kindly provided by the Spanish Meteorological Agency (AEMet). 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).

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Lindner, M.; Maroschek, M.; Netherer, S.; Kremer, A.; Barbati, A.; Garcia-Gonzalo, J.; Seidl, R.; Delzon, S.; Corona, P.; Kolström, M.; et al. Climate change impacts, adaptive capacity, and vulnerability of European forest ecosystems. For. Ecol. Manag. 2010, 259, 698–709. [Google Scholar] [CrossRef]
  2. Baldocchi, D.D. Assessing the eddy covariance technique for evaluating carbon dioxide exchange rates of ecosystems: Past, present and future. Glob. Chang. Biol. 2003, 9, 479–492. [Google Scholar] [CrossRef] [Green Version]
  3. McRoberts, R.E.; Tomppo, E.O. Remote sensing support for national forest inventories. Remote Sens. Environ. 2007, 110, 412–419. [Google Scholar] [CrossRef]
  4. Lu, D. The potential and challenge of remote sensing-based biomass estimation. Int. J. Remote Sens. 2006, 27, 1297–1328. [Google Scholar] [CrossRef]
  5. Dirección General de Conservación de la Naturaleza, Ministerio de Medio Ambiente. DGCN II Inventario Forestal Nacional 1986–1996; Dirección General de Conservación de la Naturaleza, Ministerio de Medio Ambiente: Madrid, Spain, 1996.
  6. 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]
  7. Kumar, L.; Mutanga, O. Remote sensing of above-ground biomass. Remote Sens. 2017, 9, 935. [Google Scholar] [CrossRef] [Green Version]
  8. 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]
  9. 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]
  10. 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]
  11. Monteith, J.L. Solar Radiation and Productivity in Tropical Ecosystems. J. Appl. Ecol. 1972, 9, 747–766. [Google Scholar] [CrossRef] [Green Version]
  12. Turner, D.P.; Ritts, W.D.; Cohen, W.B.; Maeirsperger, T.K.; Gower, S.T.; Kirschbaum, A.A.; Running, S.W.; Zhao, M.; Wofsy, S.C.; Dunn, A.L.; et al. Site-level evaluation of satellite-based global terrestrial gross primary production and net primary production monitoring. Glob. Chang. Biol. 2005, 11, 666–684. [Google Scholar] [CrossRef]
  13. Turner, D.P.; Ritts, W.D.; Cohen, W.B.; Gower, S.T.; Running, S.W.; Zhao, M.; Costa, M.H.; Kirschbaum, A.A.; Ham, J.M.; Saleska, S.R.; et al. Evaluation of MODIS NPP and GPP products across multiple biomes. Remote Sens. Environ. 2006, 102, 282–292. [Google Scholar] [CrossRef]
  14. Townsend, P.A. Estimating forest structure in wetlands using multitemporal SAR. Remote Sens. Environ. 2002, 79, 288–304. [Google Scholar] [CrossRef]
  15. 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; Academic Press Inc.: Cambridge, MA, USA, 1993; pp. 141–158. [Google Scholar]
  16. 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]
  17. Chiesi, M.; Maselli, F.; Moriondo, M.; Fibbi, L.; Bindi, M.; Running, S.W. Application of BIOME-BGC to simulate Mediterranean forest processes. Ecol. Model. 2007, 206, 179–190. [Google Scholar] [CrossRef]
  18. Chiesi, M.; Cherubini, P.; Maselli, F. Adaptation of a modelling strategy to predict the NPP of even-aged forest stands. Eur. J. For. Res. 2012, 131, 1175–1184. [Google Scholar] [CrossRef]
  19. Sánchez-Ruiz, S.; Chiesi, M.; Fibbi, L.; Carrara, A.; Maselli, F.; Gilabert, 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]
  20. 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. Model. 2009, 220, 330–342. [Google Scholar] [CrossRef]
  21. FAO. Global Forest Resources Assessment 2010. Terms and Definitions; FAO: Rome, Italy, 2010. [Google Scholar]
  22. 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]
  23. Moreno, A.; Gilabert, M.A.; Camacho, F.; Martínez, B. Validation of daily global solar irradiation images from MSG over Spain. Renew. Energy 2013, 60, 332–342. [Google Scholar] [CrossRef]
  24. 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]
  25. Ballabio, C.; Panagos, P.; Monatanarella, L. Mapping topsoil physical properties at European scale using the LUCAS database. Geoderma 2016, 261, 110–123. [Google Scholar] [CrossRef]
  26. IGN. Documento Técnico SIOSE 2011; IGN: San Francisco, CA, USA, 2011. [Google Scholar]
  27. 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. [Google Scholar] [CrossRef] [Green Version]
  28. 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]
  29. Iqbal, M. An Introduction to Solar Radiation; Academic Press: Toronto, ON, Canada, 1983. [Google Scholar]
  30. Roujean, J.; Breon, F. Estimating PAR Absorbed by Vegetation from Bidirectional Reflectance Measurements. Remote Sens. Environ. 1995, 51, 375–384. [Google Scholar] [CrossRef]
  31. 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]
  32. Jensen, M.E.; Haise, H.R. Estimating evapotranspiration from solar radiation. J. Irrig. Drain. Div. 1965, 89, 15–41. [Google Scholar]
  33. Running, S.W.; Zhao, M. User’s Guide: Daily GPP and annual NPP (MOD17A2/A3) Products NASA Earth Observing System MODIS Land Algorithm. Version 3.0. 2015. Available online: https://www.ntsg.umt.edu/files/modis/MOD17UsersGuide2015_v3.pdf (accessed on 24 April 2020).
  34. Odum, E.P.; Barrett, G. Fundamentals of Ecology, 5th ed.; Thompson Brooks/Cole: Belmont, CA, USA, 2005; ISBN 0534420664. [Google Scholar]
  35. Alcaraz, D.; Paruelo, J.; Cabello, J. Identification of current ecosystem functional types in the Iberian Peninsula. Glob. Ecol. Biogeogr. 2006, 15, 200–212. [Google Scholar] [CrossRef]
  36. Vicente-Serrano, S.M.; Pérez-Cabello, F.; Lasanta, T. Assessment of radiometric correction techniques in analyzing vegetation variability and change using time series of Landsat images. Remote Sens. Environ. 2008, 112, 3916–3934. [Google Scholar] [CrossRef]
  37. Valbuena-Carabaña, M.; de Heredia, U.L.; Fuentes-Utrilla, P.; González-Doncel, I.; Gil, L. Historical and recent changes in the Spanish forests: A socio-economic process. Rev. Palaeobot. Palynol. 2010, 162, 492–506. [Google Scholar] [CrossRef]
  38. Dirección General de Conservación de la Naturaleza, Ministerio de Medio Ambiente. DGCN III Inventario Forestal Nacional 1997–2007; Dirección General de Conservación de la Naturaleza, Ministerio de Medio Ambiente: Madrid, Spain, 2006.
  39. Valcarcel, N.; Villa, G.; Arozarena, A.; Garcia-Asensio, L.; Caballlero, M.E.E.; Porcuna, A.; Domenech, E.; Peces, J.J.J. SIOSE, a Successful Test Bench Towards Harmonization and Integration of Land Cover/Use Information As Environmental Reference Data. In Proceedings of the International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Beijing, China, 3–11 July 2008; Volume 37, pp. 1159–1164. [Google Scholar]
  40. LSA-SAF Product User Manual Down-Welling Surface Shortwave Flux (DSSF). 2011. Available online: http://landsaf.ipma.pt (accessed on 24 April 2020).
  41. Schaaf, C.B.; Gao, F.; Strahler, A.H.; Lucht, W.; Li, X.; Tsang, T.; Strugnell, N.C.; Zhang, X.; Jin, Y.; Muller, J.-P.; et al. Global albedo, BRDF and nadir BRDF-adjusted reflectance products from MODIS. Remote Sens. Environ. 2002, 83, 135–148. [Google Scholar] [CrossRef] [Green Version]
  42. 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]
  43. Deng, H.; Runger, G. Gene selection with guided regularized random forest. Pattern Recognit. 2013, 46, 3483–3489. [Google Scholar] [CrossRef] [Green Version]
  44. Izquierdo-Verdiguier, E.; Zurita-Milla, R.; Rolf, A. On the use of guided regularized random forests to identify crops in smallholder farm fields. In Proceedings of the 2017 9th International Workshop on the Analysis of Multitemporal Remote Sensing Images (MultiTemp), Bruges, Belgium, 27–29 June 2017; pp. 1–3. [Google Scholar]
  45. Gorelick, N.; Hancher, M.; Dixon, M.; Ilyushchenko, S.; Thau, D.; Moore, R. Google Earth Engine: Planetary-scale geospatial analysis for everyone. Remote Sens. Environ. 2017, 202, 18–27. [Google Scholar] [CrossRef]
  46. 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]
  47. Christensen, O.B.; Goodess, C.M.; Harris, I.; Watkiss, P. European and Global Climate Change Projections: Discussion of Climate Change Model Outputs, Scenarios and Uncertainty in the EC RTD ClimateCost Project. Clim. Proj. Final Rep. 2011, 1, 1–36. [Google Scholar]
  48. Chirici, G.; Chiesi, M.; Corona, P.; Puletti, N.; Mura, M.; Maselli, F. Prediction of forest NPP in Italy by the combination of ground and remote sensing data. Eur. J. For. Res. 2015, 134, 453–467. [Google Scholar] [CrossRef] [Green Version]
  49. Hasenauer, H.; Petritsch, R.; Zhao, M.; Boisvenue, C.; Running, S.W. Reconciling satellite with ground data to estimate forest productivity at national scales. For. Ecol. Manag. 2012, 276, 196–208. [Google Scholar] [CrossRef]
  50. Ryan, M.G. Effects of Climate Change on Plant Respiration Author (s): Michael G. Ryan Stable URL: http://www.jstor.org/stable/1941808. Eff. Clim. Chang. Plant Respir. 2015, 1, 157–167. [Google Scholar]
  51. Atkin, O.K.; Tjoelker, M.G. Thermal acclimation and the dynamic response of plant respiration to temperature. Trends Plant Sci. 2003, 8, 343–351. [Google Scholar] [CrossRef]
Figure 1. Flowchart of the modeling strategy to assess annual forest net primary production (NPP) described in Section 2. Input data (with yellow background) are described in Section 3.3. Input data with purple borders were derived from remotely sensed data. The quality assessment is described in Section 4.1.
Figure 1. Flowchart of the modeling strategy to assess annual forest net primary production (NPP) described in Section 2. Input data (with yellow background) are described in Section 3.3. Input data with purple borders were derived from remotely sensed data. The quality assessment is described in Section 4.1.
Remotesensing 12 01356 g001
Figure 2. (a) Forest type map of peninsular Spain. OSA, NC, EBF, LDBF, HDBF, LENF, and HENF refer respectively to out of the study area, non-classified, evergreen broadleaved forest, low altitude deciduous broadleaved forest, high altitude deciduous broadleaved forest, low altitude evergreen needleleaved forest, and high altitude evergreen needleleaved forest. (b) Average annual (period 2005–2012) forest gross primary production (GPP) obtained using Equation (1). (c) Forest GSV produced as summarized in Section 3.3 [28] and aggregated at 1-km spatial resolution.
Figure 2. (a) Forest type map of peninsular Spain. OSA, NC, EBF, LDBF, HDBF, LENF, and HENF refer respectively to out of the study area, non-classified, evergreen broadleaved forest, low altitude deciduous broadleaved forest, high altitude deciduous broadleaved forest, low altitude evergreen needleleaved forest, and high altitude evergreen needleleaved forest. (b) Average annual (period 2005–2012) forest gross primary production (GPP) obtained using Equation (1). (c) Forest GSV produced as summarized in Section 3.3 [28] and aggregated at 1-km spatial resolution.
Remotesensing 12 01356 g002aRemotesensing 12 01356 g002b
Figure 3. Scatterplots of estimated NPP (NPPEST) vs. reference NPP (NPPNFI3) per forest type. NPPEST was calculated with GPPPEM and NFI3 GSV. 1:1 line in blue. Regression line in red.
Figure 3. Scatterplots of estimated NPP (NPPEST) vs. reference NPP (NPPNFI3) per forest type. NPPEST was calculated with GPPPEM and NFI3 GSV. 1:1 line in blue. Regression line in red.
Remotesensing 12 01356 g003
Figure 4. Average annual forest NPP obtained through the modeling strategy presented in Section 2: (a) using GPP estimated using Equation (1); (b) using GPP simulated using Biome-BGC; (c) using GPP simulated using Biome-BGC, ΔT = +2 °C, and δPRE = –20%.
Figure 4. Average annual forest NPP obtained through the modeling strategy presented in Section 2: (a) using GPP estimated using Equation (1); (b) using GPP simulated using Biome-BGC; (c) using GPP simulated using Biome-BGC, ΔT = +2 °C, and δPRE = –20%.
Remotesensing 12 01356 g004aRemotesensing 12 01356 g004b
Figure 5. Variations of annual forest NPP due to possible changes in daily air temperature (ΔT) and precipitation (δPRE) caused by future climate change. One case representative of each forest type is shown.
Figure 5. Variations of annual forest NPP due to possible changes in daily air temperature (ΔT) and precipitation (δPRE) caused by future climate change. One case representative of each forest type is shown.
Remotesensing 12 01356 g005
Table 1. Summary information about the used third Spanish Forest Inventory (NFI3) plots per forest type. Growing stock volume (GSV) and current annual increment (CAI) are averaged for all plots.
Table 1. Summary information about the used third Spanish Forest Inventory (NFI3) plots per forest type. Growing stock volume (GSV) and current annual increment (CAI) are averaged for all plots.
Forest TypeSpeciesNumber of PlotsGSV (m3 ha−1)CAI (m3 ha−1 a−1)
EBFQuercus ilex266451.06
LDBFQuercus robur1821224.10
HDBFFagus sylvatica2551974.57
LENFPinus pinea122712.97
HENFPinus nigra3071003.52
Table 2. Descriptive information of the forest types considered in the present study 1.
Table 2. Descriptive information of the forest types considered in the present study 1.
Forest Typeh (m)T (°C)Rg (MJ m−2)PRE (mm)PET (mm)
EBF584
(335)
14.53
(1.97)
5725
(666)
711
(300)
1193
(235)
LDBF520
(186)
12.73
(0.97)
4863
(524)
1084
(307)
910
(139)
HDBF1099
(205)
11.48
(1.19)
5272
(420)
850
(267)
936
(118)
LENF484
(223)
14.38
(1.68)
5586
(635)
714
(360)
1153
(209)
HENF1198
(290)
12.53
(1.54)
5742
(348)
604
(183)
1089
(165)
1h is the elevation from sea level, T is the average mean annual air temperature, Rg is the average annual incoming global solar radiation, PRE is the average annual precipitation, and PET is the average annual potential evapotranspiration. Averages refer to the period considered in the present study (2005–2012) and were calculated from the datasets described in Section 3.3. Standard deviations are expressed between brackets.
Table 3. Basic density (ρ), biomass expansion factor (C), and stem carbon allocation ratio (D) for the calculation of potential GSV and CAI [46].
Table 3. Basic density (ρ), biomass expansion factor (C), and stem carbon allocation ratio (D) for the calculation of potential GSV and CAI [46].
Forest Typeρ (Mg m−3)C (m3 m−3)D
EBF0.701.450.47
LDBF0.691.330.45
HDBF0.611.360.45
LENF0.531.530.42
HENF0.381.310.42
Table 4. Statistics resulting from the comparison between estimated (NPPEST) and reference (NPPNFI3) NPP. NPPEST was calculated using Equation (4). NFI3 GSV was used in Equation (6). GPPPEM was estimated using Equation (1). GPPBGC was simulated using Biome-BGC.
Table 4. Statistics resulting from the comparison between estimated (NPPEST) and reference (NPPNFI3) NPP. NPPEST was calculated using Equation (4). NFI3 GSV was used in Equation (6). GPPPEM was estimated using Equation (1). GPPBGC was simulated using Biome-BGC.
Forest TypeR2MBE (g m−2 a−1)RMSE (g m−2 a−1)
with
GPPPEM
with
GPPBGC
with
GPPPEM
with
GPPBGC
with
GPPPEM
with
GPPBGC
EBF0.690.74604510078
LDBF0.600.566420150150
HDBF0.620.60120140160190
LENF0.590.49−170−180210230
HENF0.460.51−66−75140140
Table 5. Minimum (min), maximum (max), mean, and standard deviation (std) of the different average annual GPP products and their corresponding NPPs involved in the estimation of annual NPP per forest type. GPPPEM was estimated using Equation (1). GPPBGC was simulated using Biome-BGC. All values are expressed in g m−2 a−1.
Table 5. Minimum (min), maximum (max), mean, and standard deviation (std) of the different average annual GPP products and their corresponding NPPs involved in the estimation of annual NPP per forest type. GPPPEM was estimated using Equation (1). GPPBGC was simulated using Biome-BGC. All values are expressed in g m−2 a−1.
Forest TypeGPPPEMGPPBGC
MinMaxMeanStdMinMaxMeanStd
EBF1281778916308131495828243
LDBF3101762122920220116111156220
HDBF11785104925417417291003312
LENF7187686431101909784284
HENF415927742202051681773256
Forest TypeNPP with GPPPEMNPP with GPPBGC
MinMaxMeanStdMinMaxMeanStd
EBF887618116812470143110
LDBF3683241813817788385134
HDBF466123212814679218134
LENF1197149116749512387
HENF154912584958312387
Table 6. Minimum (min), maximum (max), mean, and standard deviation (std) of the average annual NPP obtained as explained in Section 4.2 with ΔT = +2 °C and δPRE = −20% per forest type. All values are expressed in g m−2 a−1.
Table 6. Minimum (min), maximum (max), mean, and standard deviation (std) of the average annual NPP obtained as explained in Section 4.2 with ΔT = +2 °C and δPRE = −20% per forest type. All values are expressed in g m−2 a−1.
Forest TypeMinMaxMeanStd
EBF14875197218
LDBF171136564204
HDBF15696220133
LENF7772215179
HENF862210484
Table 7. Descriptive information of the forest sites selected for Figure 5.
Table 7. Descriptive information of the forest sites selected for Figure 5.
Forest
Type
LAT
(°)
LON
(°)
h
(m)
T
(°C)
Rg
(MJ m−2)
PRE
(mm)
PET
(mm)
EBF43.6161−7.982122513.7545141366856
LDBF42.6964−2.473267912.434855743921
HDBF42.7946−2.0714113512.4750981018962
LENF40.3303−6.098257413.9159848111220
HENF41.8304−2.9821116010.375513676923
LAT is latitude, LON is longitude, h is the elevation from sea level, T is the average mean annual air temperature, Rg is the average annual incoming global solar radiation, PRE is the average annual precipitation, and PET is the average annual potential evapotranspiration. Averages refer to the period considered in the present study (2005–2012) and were calculated from the datasets described in Section 3.3.

Share and Cite

MDPI and ACS Style

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. https://doi.org/10.3390/rs12091356

AMA Style

Sánchez-Ruiz S, Maselli F, Chiesi M, Fibbi L, Martínez B, Campos-Taberner M, García-Haro FJ, Gilabert MA. Remote Sensing and Bio-Geochemical Modeling of Forest Carbon Storage in Spain. Remote Sensing. 2020; 12(9):1356. https://doi.org/10.3390/rs12091356

Chicago/Turabian Style

Sánchez-Ruiz, Sergio, Fabio Maselli, Marta Chiesi, Luca Fibbi, Beatriz Martínez, Manuel Campos-Taberner, Francisco Javier García-Haro, and María Amparo Gilabert. 2020. "Remote Sensing and Bio-Geochemical Modeling of Forest Carbon Storage in Spain" Remote Sensing 12, no. 9: 1356. https://doi.org/10.3390/rs12091356

APA Style

Sánchez-Ruiz, S., Maselli, F., Chiesi, M., Fibbi, L., Martínez, B., Campos-Taberner, M., García-Haro, F. J., & Gilabert, M. A. (2020). Remote Sensing and Bio-Geochemical Modeling of Forest Carbon Storage in Spain. Remote Sensing, 12(9), 1356. https://doi.org/10.3390/rs12091356

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