Signiﬁcant Loss of Ecosystem Services by Environmental Changes in the Mediterranean Coastal Area

: Mediterranean coastal areas are among the most threated forest ecosystems in the northern hemisphere due to concurrent biotic and abiotic stresses. These may affect plants functionality and, consequently, their capacity to provide ecosystem services. In this study, we integrated ground-level and satellite-level measurements to estimate the capacity of a 46.3 km 2 Estate to sequestrate air pollutants from the atmosphere, transported to the study site from the city of Rome. By means of a multi-layer canopy model, we also evaluated forest capacity to provide regulatory ecosystem services. Due to a signiﬁcant loss in forest cover, estimated by satellite data as − 6.8% between 2014 and 2020, we found that the carbon sink capacity decreased by 34% during the considered period. Furthermore, pollutant deposition on tree crowns has reduced by 39%, 46% and 35% for PM, NO 2 and O 3 , respectively. Our results highlight the importance of developing an integrated approach combining ground measurements, modelling and satellite data to link air quality and plant functionality as key elements to improve the effectiveness of estimate of ecosystem services.


Introduction
The provisioning of ecosystem services (ES) by urban and peri-urban forests has direct effect on the social, economic and environmental benefits. Regulating ecosystem Services (RES) represent the benefits to mankind deriving from the regulation of ecosystem processes [1]. The role of forests as nature-based solutions to mitigate the effects of climate change has been largely discussed [2], becoming a fundamental part of the proposal formulated at international level during important summits as the G20 (Rome) and COP26 (Glasgow) that fix the ambitious goals of planting 1 trillion trees (Rome G20) and stop deforestation (COP26 Glasgow) by 2030.
Estimating the contribution of "process-driven" RES such as air quality improvement, water purification, and climate regulation to human well-being is typically made by models [3], often challenged by a lack of input data, especially for estimations at broad scales [4].
Canopy models help unravel the factors that influence the response of the forest ecosystems in different environmental conditions by testing scenarios [5] or evaluating the impact of different stressors [6]. Such models need climatic and ecophysiological inputs that are often unavailable, and direct measurements needs to be performed to obtain reliable estimates [7]. Meteorological and air quality data for big cities are often available at hourly to daily resolution [8][9][10], and when missing (rural areas), high-resolution gridded

Study Area
The Presidential Estate of Castelporziano is a natural reserve of about 60 km 2 (of which 85% are natural forests), 25 km away from the south-southwest area of Rome on the Tyrrhenian Sea coast. The Estate represents a hotspot for biodiversity in the Mediterranean area, which hosts more than 1000 plant species [38,39].
Oak forests and mixed deciduous broad-leaved woods predominate, occupying 23.6 km 2 (40% of the woods), a surface higher than the evergreen oaks, and are represented by holm oaks (4.3 km 2 ) and cork oaks (2.4 km 2 ), by pine forests (9.1 km 2 ) and by Mediterranean maquis favored by anthropic land use changes in past centuries [40]. The most significant aspect that makes Castelporziano a territory now unique in the Mediterranean Basin is given by the presence of a specific and therefore very high genetic biodiversity, with 5037 species recorded, representing an exceptional wealth of flora and fauna on an area of just 60 square kilometers [41], (Figure 1). However, this valuable Green Infrastructure being located in the Metropolitan area of Rome is subjected to wide types of stressors, both natural and anthropogenic [42,43]. The intense urban sprawls nearby the Estate have caused an overexploitation of the water table and enhanced the seawater intrusion that exacerbate the natural drought that typically occurs during the summer. Moreover, high concentrations of pollutants that are transported over Castelporziano from the city center of Rome have been recorded. Indeed, the Estate receives plumes of air masses all day long from the sea and from the city of Rome because of the wind circulation, that follows a sea-land breeze regime, the dominant wind direction is S-SW during the morning and N-NE during the afternoon. The site can experience high summer levels of tropospheric ozone up to 90 ppb [44]. Soil presents a flat topography, and it can be divided in two main soil types ( Figure S1): the coastal area, characterized by a sandy texture (77% sand, 14% silt, 9% clay) and low water-holding capacity, and the inland However, this valuable Green Infrastructure being located in the Metropolitan area of Rome is subjected to wide types of stressors, both natural and anthropogenic [42,43]. The intense urban sprawls nearby the Estate have caused an overexploitation of the water table and enhanced the seawater intrusion that exacerbate the natural drought that typically occurs during the summer. Moreover, high concentrations of pollutants that are transported over Castelporziano from the city center of Rome have been recorded. Indeed, the Estate receives plumes of air masses all day long from the sea and from the city of Rome because of the wind circulation, that follows a sea-land breeze regime, the dominant wind direction is S-SW during the morning and N-NE during the afternoon. The site can experience high summer levels of tropospheric ozone up to 90 ppb [44]. Soil presents a flat topography, and it can be divided in two main soil types ( Figure S1): the coastal area, characterized by a sandy texture (77% sand, 14% silt, 9% clay) and low water-holding capacity, and the inland area, with a loamy-sandy soil type (47% sand, 29% silt, 24% clay), as shown in Figure S1 [45]. The study area is characterized by the typical Mediterranean climate, with pronounced seasonality. Summers are hot and dry, and winters are moderately cold. Precipitations occur prevalently during winter, spring and autumn. Annual precipitation up to 1019 ± 105.5, 854 ± 130.8 and 507 ± 95.6 mm y −1 was observed in the 3 years of study in 2014, 2019 and 2020, respectively, considering averaged values for all the meteorological stations in the Estate. During the three years of study, mean annual temperature was of 16 • C with seasonal peaks in summer 2019 ( Figure S2).

The AIRTREE Model
The Aggregated Interpretation of the Energy balance and water dynamics for Ecosystem services assessment (AIRTREE) model [14] is a one-dimensional multi-layer model that couples soil, plant and atmospheric processes to predict exchanges of CO 2 , water, ozone and particulate matter (PM) between leaves and the atmosphere and integrates them through five layers to obtain fluxes at the canopy level. The model allows to estimate stomatal conductance (gs) and photosynthesis (A) through the Farquhar-Von Caemmerer-Berry model of photosynthesis (FvCB model) and the Ball, Woodrow and Berry stomatal conductance (BWB) model at different levels from the top to the bottom of the canopy [46]. The model accounts for oxidative limitations (i.e., drought and ozone stress) to gas exchanges (A and gs) through linear functions of soil water content [47] and stomatal ozone [48][49][50].
In this study NO 2 , SO 2 and CO fluxes were calculated as: where F is the pollutant flux (µmol m 2 s −1 ), V d is the deposition velocity (m s −1 ) and air pollutant concentration (µmol m −3 ). V d was estimated through a resistance scheme [51,52]: where R a is the atmospheric resistance (see Fares et al. [14] for details), R b is the boundary layer resistances (s m −1 ) calculated according Pederson et al. [53] and R c is the canopy resistance calculated according to Bidwell & Fraser [54] as: where Sc is the Schmidt number, Pr is the Prandtl number (0.72), k is the von Karman constant (0.41), u is wind speed (m s −1 ), r s (s m −1 ) is the inverse of g s , the formulation of which is extensively described in Fares et al. [14]. Soil resistance r soil (s m −1 ) was set equal to 2941 and 2000 (m s −1 ) during vegetative and non-vegetative periods, respectively. Mesophyll resistance r m (s m −1 ) was set to zero for SO 2 [55] and 600 for NO 2 [56]. Cuticular resistance r t (s m −1 ) was set to 8 × 10 3 and 2 × 10 4 for SO 2 and NO 2 , respectively [56]. Rc for CO was set equal to 5 × 10 4 and 1 × 10 6 during vegetative and non-vegetative periods, respectively [54]. Concerning PM deposition, the model incorporates traditional atmospheric models to predict particle deposition extensively described in Fares et al. [10,14].

Acquisition of Meteorological Variables
For the entire Estate, we performed a Voronoi tessellation using data from a network of 6 meteorological stations where continuous monitoring of local climatic conditions have been carried out ( Figure S3). Each station is equipped with sensors measuring air temperature and relative humidity, wind speed and direction, precipitation and solar radiation; for details on instrumentation, see Aromolo et al. [57]. In addition, the Estate hosts a flux tower for continuously monitoring gas exchanges (i.e., Carbon dioxide CO 2 , ozone O 3 , nitrogen dioxide NO 2 , particulate matter PM) between the vegetation and the atmosphere in the coastal area, 1.5 km away from the Thyrrenean sea (see [30,44] for details). Just outside the Estate, within 1.5 km, there are 4 air quality monitoring stations maintained by the regional agency for the environmental protection (ARPA Lazio) that continuously monitor concentration of PM, NO 2 , SO 2 , CO, O 3 and other pollutants. Those air quality data were used as input for the AIRTREE model. The data were assumed to be exemplificative of Castelporziano site since a comparison between PM 2.5 concentration measured in 2014 and 2015 at the ICOS It-Cp2 flux site and average values measured at the 5 ARPA stations showed a good correlation (Pearson's r > 0.7, data not shown). The Voronoi tessellation method allowed us to generate a spatial-metric decomposition determined by distances to discrete set of elements in space (meteorological stations), and therefore, 6 Voronoi polygons were created for the whole Estate.

Acquisition of Vegetation Map
The vegetation map of the Castelporziano Estate ( Figure 1) used in this study [58] is characterized by 7 groups of forested land cover type (LCT) (Other broadleaves; Low shrubs; Holm oak forests; Mediterranean maquis; Pine forests; Deciduous oak forests; Corks) and 9 non-forested LCT (Internal waters; Permanently non-productive areas; Permanently non-productive natural areas). To associate species-specific ecophysiological traits to each of the forested LCT group, we used the state-of-the-art of literature about the vegetation in Castelporziano to identify the most representative species included in each group [38]. Therefore, the values of carbon, PM and other pollutants deposition for each pixel (associated to a specific LCT) are the results of a single model simulation. For each model run, the AIRTREE model was parametrized for an ideal tree-type by synthetizing the ecophysiological characteristics of the dominant trees representative of the LCT (Table S1).

Retrieval of Biometric Vegetation Data
A raster map of the species associations was superimposed and aligned with the maps of Leaf Area Index (LAI), canopy cover and heights ( Figure S4). Concerning LAI and canopy cover, we used 300 m resolution satellite images from Copernicus Global Land Service, the Earth Observation programme of the European Commission [59]. For each year of study (i.e., 2014, 2019, 2020), the maximum values observed during the growing season were used. AIRTREE simulations were conducted for each pixel of the study area, and for each year of study, we extracted the LAI values filtering out pixels classified as buildings, urban, wetlands, meadows and artificial settlements. LAI is a critical value for the proper evaluation of ecosystem services via AIRTREE model [14]. Therefore, we compared satellite data with field data observations collected with multiple simultaneous measurements of transmittance using the LAI-2200C Plant Canopy Analyzer (Li-Cor, Lincoln, NE, USA) and with a method based on the use of hemispherical photos collected in the frame of ICOS activities [60]. The sampling points were selected through a stratified sampling procedure ( Figure S5) over the main four LCTs: Mediterranean maquis, Deciduous broadleaves, Pine forest and Holm oak forest (Table S2). For more details on the statistical design and the field campaigns for LAI measurement, please see the section "In situ LAI measurements" in the Supplementary File. To derive a function suitable to provide a realistic estimate of the annual dynamics of LAI for each LCT group, monthly values of LAI were fit to a set of parametric non-linear models (Gaussian, Exponential, Rational, Power, Sin, Weibull and Fourier). The function that best fitted data (lowest RMSE) was implemented in the model. The Fourier model was identified as the most appropriate model to reproduce LAI dynamics through the three years of study for the LCTs functional groups: where a, b, c and d were LCT-specific fixed values (i.e., the coefficients of the fitted model), and x is the day of the year ( Figure 2). Coefficients of the Fourier models used for each LCT and for each of the three years of study are summarized in Table S3. Non-linear model intercomparisons are shown in Tables S4-S24.
ence year in order to derive the daily LAI at each model time step. Canopy heights were only available for the year 2019 at a spatial resolution of 25 m [61]. These were resampled to 300 m in order to align the maps with those of LAI and canopy cover. All maps out on the boundaries of the Estate were cut, and pixels falling into the nonforested LCT were removed, so as to work only on the actual forested areas.  Table S1.

LAI and Canopy Cover Dynamics
Overall, average LAI for the entire Estate was 4.88, 4.31 and 3.90 in the 3 years of study, respectively ( Figures S4 and S6), indicating a decrease of 11.83% and 20.05% in 2019 and 2020 compared to 2014, respectively. LAI showed a decreasing trend for all LCTs (Figure 3), in particular for Pine forest (for which the representability of the total LAI changes from 14% in 2014 to 12% in 2020, Figure S5). In 2020, LAI increased for low shrubs by only 1% when compared to 2019. This is also confirmed by canopy cover data ( Figure S7) showing an average loss of 6.83% between 2014 and 2020. In particular, the Pine forest showed a loss of cover of 13.33%, of which 9.83% occurred in a year, between 2019 and 2020. Indeed, a positive line-  Table S1.
Such a function was finally implemented into the AIRTREE model for each reference year in order to derive the daily LAI at each model time step. Canopy heights were only available for the year 2019 at a spatial resolution of 25 m [61]. These were resampled to 300 m in order to align the maps with those of LAI and canopy cover. All maps out on the boundaries of the Estate were cut, and pixels falling into the non-forested LCT were removed, so as to work only on the actual forested areas.

LAI and Canopy Cover Dynamics
Overall, average LAI for the entire Estate was 4.88, 4.31 and 3.90 in the 3 years of study, respectively ( Figures S4 and S6), indicating a decrease of 11.83% and 20.05% in 2019 and 2020 compared to 2014, respectively. LAI showed a decreasing trend for all LCTs (Figure 3), in particular for Pine forest (for which the representability of the total LAI changes from 14% in 2014 to 12% in 2020, Figure S5). In 2020, LAI increased for low shrubs by only 1% when compared to 2019. This is also confirmed by canopy cover data ( Figure S7) showing an average loss of 6.83% between 2014 and 2020. In particular, the Pine forest showed a loss of cover of 13.33%, of which 9.83% occurred in a year, between 2019 and 2020. Indeed, a positive linear correlation between changes in LAI and changes in canopy cover was observed (Pearson's r = 0.63). Although an increase in forest cover of 2% has been observed between 2010 and 2015 in the Mediterranean basin [62], the loss of Cover (7%) and LAI (up to 20%) experienced by the vegetation of the Estate between 2014 and 2020 confirms the increasing vulnerability of forests due to the climatic stressors and aridity the vegetation is exposed to [63]. Indeed, from 2014 (previously described wet year [30]), warmer and drier years succeeded. This may have caused stress and loss of productivity for Mediterranean forests, even if they are adapted to the typical dry and hot summers [64][65][66][67]. All ecosystems except shrubs showed a reduction in their LAI up to 14% (Figure 3).
The environmental monitoring data concerning the health status of the Castelporziano forests, detected from 2015 to today through the diachronic interpretation of the NDVI values provided by Sentinel-2 data, showed widespread suffering in pine and deciduous oaks forests [68]. Regarding the pine forests, two insect infestations occurred in 2015 and 2017, causing a widespread decline in the photosynthetic capacity (−34% in the average NDVI values for the period 2015-2021) of the forest, which caused the death of plants over large area [68]. In Deciduous oak forests, the limiting factor seems to be strongly correlated with summer drought and springs characterized by low amount of rain, causing a mean decrease in NDVI values of −27% for the observation period [37,68].
Interestingly, there is high resilience from several species of Evergreen shrubs, capable of adapting to drought stress even more than Quercus ilex has been previously observed [69,70], and therefore, a transformation of forest structure into a shrubland under increasing heatwaves and drought is the most likely scenario [64,71]. Indeed, Mediterranean plants have developed various morphological and physiological strategies to adapt to drought [72], but the largest trees are frequently more sensitive and less resistant and resilient to increases in aridity [72,73]. In line with our observations, Lloret et al. [74] showed that shorter trees are more resistant and resilient to increases in drought than taller trees. Tall shrubs such as Phillyrea sp., which are abundant in the Estate of Castelporziano, have been shown to have a higher capacity than trees to adapt and resist intensive droughts [75,76] and have a higher capacity than sympatric trees to maintain their foliage and concentrations of non-structural carbohydrates after droughts [76].
The inter-comparison between non-linear models fitted to monthly values of LAI for each year of study and for each LCT is shown in Tables S4-S24. On the 21 fits performed (i.e., 3 years for 7 LCTs), the Fourier function ( Figure 3 and Table S3) was found to be the best model (highest R 2 adj and lowest RMSE) 12 times, although it produced some minor biases (i.e., increase in LAI for "other broadleaves" and "Corks" at the end of the year).
A general change in phenology was observed, in 2020 compared to 2014, with an anticipated peak in LAI that shifted from June to July in 2014 to May to June 2020 ( Figure 3). This is in line with previous findings showing that key Mediterranean species such as Quercus ilex display a longer period of growth by approximately 10 days by advancing the onset of spring by winter warming, with an early cessation of growth in spring and summer [77]. 2015 in the Mediterranean basin [62], the loss of Cover (7%) and LAI (up to 20%) experienced by the vegetation of the Estate between 2014 and 2020 confirms the increasing vulnerability of forests due to the climatic stressors and aridity the vegetation is exposed to [63]. Indeed, from 2014 (previously described wet year [30]), warmer and drier years succeeded. This may have caused stress and loss of productivity for Mediterranean forests, even if they are adapted to the typical dry and hot summers [64][65][66][67]. All ecosystems except shrubs showed a reduction in their LAI up to 14% (Figure 3).  Satellite LAI values used in this study were compared with field measurements with an overestimation of 20% considering all LCTs (Figure 4). Deciduous broadleaves and Holm oak forests showed the highest LAI values (4.12 and 4.41, respectively; statistics are shown in Table S2), in accordance with what was previously found for Deciduous oaks stands [78]. We measured lower LAI values for Pine forests and Mediterranean maquis (2.08 and 3.85, respectively), in line with previous findings [79,80].
Compared to 2014, the Estate showed a decrease in productivity by 31.6% and 34.3% in 2019 and 2020, respectively (Figures 5a and S10). Between 2019 and 2020, Deciduous Oak and Mediterranean maquis and Low shrubs showed an increase in carbon uptake by 5.7%, 2.6% and 22.5%, respectively (in line with our previous statements about highest tolerance to drought stress).
Compared to 2014, the Estate showed a decrease in productivity by 31.6% and 34.3% in 2019 and 2020, respectively (Figures 5a and S10). Between 2019 and 2020, Deciduous Oak and Mediterranean maquis and Low shrubs showed an increase in carbon uptake by 5.7%, 2.6% and 22.5%, respectively (in line with our previous statements about highest tolerance to drought stress). Considering that pollutants such as ozone have been found to reduce carbon assimilation in the Holm oak forest up to 10% [14], we do not exclude a possible beneficial effect of reducing atmospheric pollutants in 2020 due to lockdown status in the first part of the year. Other LCTs instead, showed a decrease in carbon uptake up to 17%, and such decrease can be highly appreciated by comparing 2014 and 2020 ( Figure S9).
We explain such behavior considering the different temperature recorded during the warm seasons in these two dry years. Otherwise, in 2019, the average temperature during summer was higher than in 2020 ( Figure S2), and accordingly, the high vapor pressure deficit may have caused stomatal closure and reduced photosynthesis, as previously discussed in Conte et al. [30]. Looking at single LCT, the data show ( Table 1) that broadleaf species and Holm oaks are the ecosystems with the higher NPP, with values up to 2.46 kg CO2 m 2 y −1 in the wet year (2014), with abrupt decreases in productivity in the dry years. These model results, especially for Holm oak, are in line with data measured with Eddy Covariance at the ICOS Holm oak It-Cp2 site at the Estate (data not shown). Like Italy, warmer and drier conditions linked to increases in Atlantic multidecadal oscillations are associated with the increase in post-1990 defoliation in the forests of Spain [81]. Our results are in agreement with previous studies that have clearly indicated that the decreases such as dieback, defoliation and lower growth in Mediterranean oaks (Quercus spp.) and pines (Pinus spp.) in southern Europe are mainly due to more frequent drought, often interacting with higher temperatures (higher water demand) and pathogenic attack [82][83][84][85][86]. Such a cumulative effect of stressors has been associated with pathogen infestations in Mediterranean forests [86][87][88][89]. Considering that pollutants such as ozone have been found to reduce carbon assimilation in the Holm oak forest up to 10% [14], we do not exclude a possible beneficial effect of reducing atmospheric pollutants in 2020 due to lockdown status in the first part of the year. Other LCTs instead, showed a decrease in carbon uptake up to 17%, and such decrease can be highly appreciated by comparing 2014 and 2020 ( Figure S9).
We explain such behavior considering the different temperature recorded during the warm seasons in these two dry years. Otherwise, in 2019, the average temperature during summer was higher than in 2020 ( Figure S2), and accordingly, the high vapor pressure deficit may have caused stomatal closure and reduced photosynthesis, as previously discussed in Conte et al. [30]. Looking at single LCT, the data show ( Table 1) that broadleaf species and Holm oaks are the ecosystems with the higher NPP, with values up to 2.46 kg CO 2 m 2 y −1 in the wet year (2014), with abrupt decreases in productivity in the dry years. These model results, especially for Holm oak, are in line with data measured with Eddy Covariance at the ICOS Holm oak It-Cp2 site at the Estate (data not shown). Like Italy, warmer and drier conditions linked to increases in Atlantic multidecadal oscillations are associated with the increase in post-1990 defoliation in the forests of Spain [81]. Our results are in agreement with previous studies that have clearly indicated that the decreases such as dieback, defoliation and lower growth in Mediterranean oaks (Quercus spp.) and pines (Pinus spp.) in southern Europe are mainly due to more frequent drought, often interacting with higher temperatures (higher water demand) and pathogenic attack [82][83][84][85][86].

Sequestration of Pollutants
Total yearly estimates of ozone deposition (O 3 ) indicate that the Estate sequestrated 159, 145 and 103 t O 3 y −1 for the years 2014, 2019 and 2020, respectively (Figures 2 and S8). Deciduous Oaks, Pine forests and Holm oak forests were the LCTs that contributed most. Deciduous Oaks contributed up to 46% to total ozone deposition (70 t y −1 ) and Pine forest up to 33% (51 t y −1 ). The contribution to total ozone deposition by the other ecosystems (Holm oak forest up to 9%, Mediterranean maquis up to 7%, Other broadleaves up to 3%, corks up to 4% and low shrubs <1%) is below 25% (Figures S11 and S12). The modelled ozone deposition on Holm oak agrees with our previous measured with the Eddy Covariance carried out in 2013 and 2014 [43].
Compared to 2014, the Estate showed a decrease in ozone deposition by 8.8% and 35.2% in 2019 and 2020, respectively. In particular, Pine forest, Deciduous oaks, Other broadleaves and Corks showed a decrease in ozone uptake by 32.7%, 25.2%, 29.4% and 35.3%, respectively. Such a decrease in deposition (Figure 3), in part due to the loss of canopy cover and in part due to physiological limitations as demonstrated by NPP decreases (Table 1), suggests that land use changes and environmental stressors heavily compromise the capacity of the Estate and, in general, Mediterranean peri-urban forests. We must point out, however, that most of the decreases in 2020 compared with 2014 are due to a strong reduction in atmospheric ozone concentration due to the lockdown status in the first part of the year. Nevertheless, ozone concentrations were similar between 2014 and 2019, highlighting evidence of stomatal limitation to ozone deposition in the dry year as stressed in our previous work [26,90].
Total yearly estimates of particulate matter (PM) indicate that the Estate sequestrated 357, 241 and 215 t PM 10 y −1 and 27, 21, and16 t PM 2.5 y −1 for the years 2014, 2019 and 2020, respectively ( Figures S8 and S13a,b). Deciduous oaks and Pine forests were the LCTs that contributed equally to PM 10 deposition (up to 41% Pine forest in 2014 and Deciduous oaks in 2020, Figure S13), while a higher contribution (up to 52%) of Pine forest was found for PM 2.5 ( Figure S13). The contribution to total PM deposition by the other ecosystems is below 25% (Other broadleaves up to 3%, Corks up to 4%, Mediterranean maquis up to 6%, Holm oak forest 8% and low shrubs <1%) ( Figure S13).
Compared to 2014, the Estate showed a decrease in PM 10 deposition by 32.5% in 2019 and 39.6% in 2020 (Figure 5b), while there was a decrease in PM 2.5 deposition of 25% in 2019 and 40% in 2020. Between 2019 and 2020, a decrease of 10.5% and 19.9% was observed for PM 10 and PM 2.5 , respectively. In particular, the Pine forest showed a reduction of 19.8% and 26% for PM 10 and PM 2.5 in 2020, respectively (Figures 2, S14 and S15). The same conclusions drawn for ozone work here, with a significant reduction in PM emissions during the lockdown in 2020 ( Figure S2). To note, several hectares (0.7 km 2 in 2017 and 2022) of forests have been cut in 2017 and 2020 (and another 1 km 2 is going to be cut in the near future) probably due to a strong infection of Tomicus destruens and Toumayella parvicornis on Pinus pinea stands, as visible in the white pixels in Figure S13. Moreover, as shown in Table 1 (and also found elsewhere by Fares et al. [10]), individual performances in PM removal are higher in Pine forests; therefore, a loss of these forest stands represents a major loss in RES provision. This loss is particularly relevant for cities such as Rome, where Pinus pinea is a key tree of the urban landscape.
Concerning Nitrogen Dioxide (NO 2 ), we found values of 92, 76 and 49 t NO 2 y −1 sequestration for the years 2014, 2019 and 2020, respectively (Figures 2 and S8). Deciduous Oaks removed up to 45%, and Pine forests removed up to 34%. Mediterranean maquis and holm oak forest contribution to total fluxes were constant (up to 8%) in the three years of study ( Figure S16). The contribution to total NO 2 deposition by the other ecosystems (Other broadleaves up to 3%, Corks up to 4% and Low shrubs < 1%) was below 10% ( Figure S16). Compared to 2014, the Estate showed a decrease in NO 2 deposition of 17% and 46.8% in 2019 and 2020, respectively. In particular, the highest decreases in NO 2 deposition were observed for Holm oak forest (up to 36.2%) and Corks (up to 38%), Pine forest (37.9%) and Other broadleaves (36.9%). Similar to ozone, the same conclusions can be drawn for NO 2 , since 2020 recorded much lower emissions due to lockdown status ( Figures S2 and S17).
As a minor entity, the Estate sequestrated other pollutants such as sulfur dioxide (SO 2 ) and carbon monoxide (CO). In agreement with the removal dynamics observed for the other pollutants, a strong decrease was observed, with values decreasing from 4.4 to 1.6 and 1.26 t SO 2 y −1 and from 0.26 to 0.25 and 0.18 t CO y −1 for the years 2014, 2019 and 2020, respectively ( Figure 3).

Conclusions
The Mediterranean Basin is a global hotspot of biological diversity and the most diverse biome in Europe. However, biotic and abiotic stressors can compromise survival of native ecosystems, especially in highly urbanized contexts. This motivated us to investigate which environmental factors affect Mediterranean forests' health status and whether natural and anthropogenic stressors can have an impact on RES. Our study is the result of the integration of remote sensing products with a mechanistic modelling approach to estimate plant functionality and stress response. Despite previous results showing that an increase in forest cover by 2% has been observed between 2010 and 2015 in the Mediterranean region [91], we highlight that in response to climatic changes, pollution and biotic stresses, not only the extensiveness of peri-urban forests can be reduced (we described a loss of canopy cover by 7% and of LAI up to 20% between 2014 and 2020) but also their capacity to deliver RES. Our results also warn that future forest composition may be altered with an increase in Mediterranean shrubs in place of forests stands populated by pines and oaks.
The obtained results can have important implications for future forest management, and to raise awareness on the high ecological and economic value of RES. In this context, the quantitative estimates of RES may play a key role in supporting decision makers and management planners to evaluate possible future impacts of different practices or environmental policies.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/f13050689/s1, Figure S1: Soil classification for the Castelporziano natural reserve derived from the Harmonized World Soil Database viewer (HWSD v.1.2) is reported. Figure S2: Meteorological inter-comparison between the three years of study seasonal values of mean temperatures, and cumulated precipitation (top left and top right, respectively) are shown together with average concentrations of PM, NO2, SO2 and CO. Figure S3: Spatial association of meteorological stations to homogeneous portions of the natural reserve. Yellow dots indicate the position of each one of the meteorological stations. Climatic conditions of the portions of the Estate (red shapes) were associated to each meteorological station. Figure S4: Superimposed maps of vegetation (colored shapes) to LAI (black and white pixels). Colored shapes represent homogeneous groups of vegetation (LCTs). Figure S5: Castelporziano map, showing the sampling points chosen randomly using QGis. Figure S6: Average values of LAI for the entire Estate of Castelporziano are reported for each of the three years of study (top). Relative contribution to overall LAI by each LCT is shown in pie charts (bottom), for 2014, 2019, and 2020, respectively. Figure S7: Bar plot showing the percent change of canopy cover values from satellite data between 2014 and 2019, in blue, and 2020 in yellow. Figure S8: Total estimates of LAI, NPP, and pollutants deposition for the whole Estate in the three years of study. Figure S9: Cumulated values of NPP for the entire Estate of Castelporziano are reported for each of the three years of study (top). Figure S10: Superimposed maps of vegetation (colored shapes) to NPP (black and white pixels). Figure S11: Superimposed maps of vegetation (colored shapes) to O3 fluxes (black and white pixels). Figure S12: Cumulated values of O3 deposition for the entire Estate of Castelporziano are reported for each of the three years of study (top). Figure S13: Superimposed maps of vegetation (colored shapes) to PM10 (left) and PM2.5 (right) deposition on canopies (black and white pixels). Figure S14: Cumulated values of PM10 deposition for the entire Estate of Castelporziano are reported for each of the three years of study (top). Figure S15: Cumulated values of PM2.5 deposition for the entire Estate of Castelporziano are reported for each of the three years of study (top). Figure S16: Cumulated values of NO2 deposition for the entire Estate of Castelporziano are reported for each of the three years of study (top). Figure S17: Superimposed maps of vegetation (colored shapes) to NO2 deposition on canopies (black and white pixels). Table S1: For each LCT, the corresponding association is reported. For each association, average ecophysiological parameters (e.g. Vcmax and ozone tolerance) of the dominant species are used as model input to characterize the LCT ideal tree-type. Table S2: LAI, (m2m-2) measured for the different LCTs in Castelporziano Estate. Table S3: Coefficient of the Fourier model used in this study to simulate LAI at each model time-step. For each LCT, yearly coefficient that best (lowest RMSE) fit LAI measured by satellite data are shown. Table S4: Intercomparison between nonlinear models used to fit the day of the year (doy) and LAI of Corks for the year 2014. Table S5: Intercomparison between nonlinear models used to fit the day of the year (doy) and LAI of Other Broadleaves for the year 2014. Table S6: Intercomparison between nonlinear models used to fit the day of the year (doy) and LAI of Deciduous Oaks for the year 2014. Table S7: Intercomparison between nonlinear models used to fit the day of the year (doy) and LAI of Holm oak forest for the year 2014. Table S8: Intercomparison between nonlinear models used to fit the day of the year (doy) and LAI of Mediterranean maquis for the year 2014. Table S9: Intercomparison between nonlinear models used to fit the day of the year (doy) and LAI of Pine forest for the year 2014. Table S10: Intercomparison between nonlinear models used to fit the day of the year (doy) and LAI of Low shrubs for the year 2014. Table S11: Intercomparison between nonlinear models used to fit the day of the year (doy) and LAI of Other broadleaves for the year 2019. Table S12: Intercomparison between nonlinear models used to fit the day of the year (doy) and LAI of Corks for the year 2019. Table S13: Intercomparison between nonlinear models used to fit the day of the year (doy) and LAI of Deciduous oaks for the year 2019. Table S14: Intercomparison between nonlinear models used to fit the day of the year (doy) and LAI of Holm oak forest for the year 2019. Table S15: Intercomparison between nonlinear models used to fit the day of the year (doy) and LAI of Mediterranean maquis for the year 2019. Table S16: Intercomparison between nonlinear models used to fit the day of the year (doy) and LAI of Pine forest for the year 2019. Table S17: Intercomparison between nonlinear models used to fit the day of the year (doy) and LAI of Low shrubs for the year 2019. Table S18: Intercomparison between nonlinear models used to fit the day of the year (doy) and LAI of Other broadleaves for the year 2020. Table S19: Intercomparison between nonlinear models used to fit the day of the year (doy) and LAI of Corks for the year 2020. Table S20: Intercomparison between nonlinear models used to fit the day of the year (doy) and LAI of Deciduous oaks for the year 2020. Table S21: Intercomparison between nonlinear models used to fit the day of the year (doy) and LAI of Holm oak forest for the year 2020. Table S22: Intercomparison between nonlinear models used to fit the day of the year (doy) and LAI of Mediterranean maquis for the year 2020. Table S23: Intercomparison between nonlinear models used to fit the day of the year (doy) and LAI of Pine forest for the year 2020. Table S24: Intercomparison between nonlinear models used to fit the day of the year (doy) and LAI of Low shrubs for the year 2020.