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

: 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 e ﬃ ciency 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 e ﬃ ciency model was used to estimate daily gross primary production (GPP), while Biome-BGC was used to simulate growth ( R G ) and maintenance ( R M ) respirations. To account for actual forest conditions, GPP, R G , and R M 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. Coe ﬃ cients of determination ranged from 0.46 to 0.74 depending on the forest type. A simpliﬁed 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 di ﬀ erent climate change scenarios reported in the literature. The modiﬁed meteorological series were used to simulate new GPP, R G , and R M through Biome-BGC and corrected using GSV. Precipitation was conﬁrmed 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.


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 (R A ), 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 R A [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 m 3 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.

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 (f APAR ), precipitation, and air temperature: 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]. f APAR 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 k 0 , k 1 , and k 2 for an optimal angular geometry in the solar principal plane. f APAR was calculated as a linear function of the renormalized difference vegetation index ε W is the water stress coefficient (C WS ) as used in [31] that accounts for short-term water stress 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]. tree cover and GSV and are calculated as follows: and v = 50 GSV ρ C / σ (6) where LAI (m 3 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 m 3 ha −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.

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 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.
NPP is the difference between GPP and R A , with R A equal to the sum of growth (R G ) and maintenance (R M ) 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]: 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: and v = 50 GSV ρ C/σ where LAI (m 3 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 m 3 ha −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).

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.

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).
Remote Sens. 2020, 12, x FOR PEER REVIEW 7 of 20 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.   (1). (c) Forest GSV produced as summarized in Section 3.3 [28] and aggregated at 1-km spatial resolution.

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. Table 3. Basic density (ρ), biomass expansion factor (C), and stem carbon allocation ratio (D) for the calculation of potential GSV and CAI [46].  Table 1. For each of these pixels, a reference annual NPP was calculated from NFI3 CAI as  (1). (c) Forest GSV produced as summarized in Section 3.3 [28] and aggregated at 1-km spatial resolution.
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. Table 2. Descriptive information of the forest types considered in the present study 1 . 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 k 0 , k 1 , and k 2 . f APAR 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 f APAR was linearly interpolated to obtain daily images.

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

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 (GPP PEM , Figure 2b) and the one simulated by Biome-BGC (GPP BGC ) [19], which also simulated R G and R M . In Equation (6), the GSV map and the data reported in Table 3 [46] were used. Table 3. Basic density (ρ), biomass expansion factor (C), and stem carbon allocation ratio (D) for the calculation of potential GSV and CAI [46].  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 NPP NFI3 = 50 CAI ρ C/D D being the stem carbon allocation ratio (Table 3), and the 50 factor accounting for both the unit conversion to m 3 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 NPP NFI3 by means of the following statistics calculated by forest type: coefficient of determination (R 2 ), mean bias error (MBE), and root mean square error (RMSE).

Simulation of Forest NPP in Future Condition
The simulation of GPP, R G , and R M 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 CO 2 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, R G , and R M 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.

NPP in Current Condition
The assessment of estimated NPP per forest type is illustrated by the scatterplots in Figure 3. The largest NPP NFI3 explained variance amongst forest types was obtained for EBF with almost 70% (R 2 = 0.69);~60% (R 2 = 0.60, R 2 = 0.62, and R 2 = 0.59 respectively) was obtained for LDBF, HDBF, and LENF; and slightly less than 50% (R 2 = 0.46) for HENF. Small scattering was observed for EBF except for a great NPP NFI3 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 R 2 = 0.59 with GPPPEM to R 2 = 0.49 with GPPBGC). Around half this difference was found for EBF, LDBF, and HENF (respectively from R 2 = 0.69, R 2 = 0.60, and R 2 = 0.46 with GPPPEM to R 2 = 0.74, R 2 = 0.56, and R 2 = 0.51 with GPPBGC), while only 2 pp were found for HDBF (from R 2 = 0.62 with GPPPEM to R 2 = 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). The accuracy statistics resulting from the NPP assessment are reported in Table 4. The greatest difference between the two NPP EST to explain the NPP NFI3 variance was found for LENF with 10 percentage points (pp) (from R 2 = 0.59 with GPP PEM to R 2 = 0.49 with GPP BGC ). Around half this difference was found for EBF, LDBF, and HENF (respectively from R 2 = 0.69, R 2 = 0.60, and R 2 = 0.46 with GPP PEM to R 2 = 0.74, R 2 = 0.56, and R 2 = 0.51 with GPP BGC ), while only 2 pp were found for HDBF (from R 2 = 0.62 with GPP PEM to R 2 = 0.60 with GPP BGC ). More of the NPP NFI3 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 NPP NFI3 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). Table 4. Statistics resulting from the comparison between estimated (NPP EST ) and reference (NPP NFI3 ) NPP. NPP EST was calculated using Equation (4). NFI3 GSV was used in Equation (6). GPP PEM was estimated using Equation (1). GPP BGC was simulated using Biome-BGC.  Figure 4a,b show, respectively, the average annual forest NPP calculated using Equation (4) with GPP PEM and with GPP BGC . 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 GPP BGC was generally smaller than NPP calculated with GPP PEM . Table 5 summarizes the basic statistics of the mentioned GPPs and their corresponding NPPs per forest type. GPP PEM maxima are greater than GPP BGC maxima for broadleaved forest, but the opposite relationship occurred for needleleaved forests. GPP PEM means were greater than GPP BGC 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.

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. 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 .  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 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%. 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. GPP PEM was estimated using Equation (1). GPP BGC was simulated using Biome-BGC. All values are expressed in g m −2 a −1 .

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 GPP BGC , ∆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 GPP BGC 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. 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 . 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.   LAT is latitude, LON is longitude, h is the elevation from sea level, T is the average mean annual air temperature, R g 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)(2006)(2007)(2008)(2009)(2010)(2011)(2012) and were calculated from the datasets described in Section 3.3.

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, R G , and R M . 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 (GPP PEM ) and the one simulated using Biome-BGC (GPP BGC ), 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 GPP PEM and GPP BGC , 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 (NPP NFI3 ) calculated from CAI observations derived from forest inventory measurements.
The strength of this approach is to combine the directly derived GPP PEM 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 R 2 and smallest errors were obtained for EBF, while the smallest R 2 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 GPP BGC in the calculation of NPP in future condition is justified.

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 Tables 5 and 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).

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 km 2 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: 1.
When applied in current climate conditions, the integrated modeling approach is capable of reproducing the general NPP variability of broadleaved forests in Spain (R 2 > 0.60) but underestimates the NPP of needleleaf forests.

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