On the Impacts of Historical and Future Climate Changes to the Sustainability of the Main Sardinian Forests

Abstract

The Mediterranean Basin is affected by climate changes that may have negative effects on forests. This study aimed to evaluate the ability of 17 forests located in the Island of Sardinia to resist or adapt to the past and future climate. Sardinia is experiencing a decreasing anthropic pressure on forests, but drought-triggered dieback in trees was recently observed and confirmed by the analysis of 20 years of satellite tree-cover data (MOD44B). Significant negative trends in yearly tree cover have affected the broad-leaved vegetation, while significative positive trends were found in the bushy sclerophyllous vegetation. Vegetation behavior resulted in being related to the mean annual precipitation (MAP); for MAP < 700 mm, we found a decline in the tall broad-leaved stands and an increase in the short ones, and the opposite was found for bushy sclerophyllous vegetations. In forests with MAP > 700 mm, both stands are stable, regardless of the growing trends in the vapor-pressure deficit (VPD) and temperature. No significative correlation between bushy sclerophyllous tree cover and the climate drivers was found, while broad-leaved tree cover is positively related to MAP1990–2019 and negatively related to the growing annual VPD. We modeled those relationships, and then we used them to coarsely predict the effects of twelve future scenarios (derived from HADGEM2-AO (CMIP5) and HadGEM3-GC31-LL (CMIP6) models) on forest tree covers. All scenarios show an annual VPD increase, and the higher its increase, the higher the trees-cover loss. The future changes in precipitation were contrasting. SC6, in line with past precipitation trends, predicts a further drop in the mean annual precipitation (−7.6%), which would correspond to an average 2.1-times-greater reduction in the tree cover (−16.09%). The future changes in precipitation for CMIP6 scenarios agree on a precipitation reduction in the range of −3.4% (SC7) to −14.29% (S12). However, although the reduction in precipitation predicted in SC12 is almost double that predicted in SC6, the consequent average reduction in TC is comparable and stands at −16%. On the contrary, SC2 predicts a turnaround with an abrupt increase of precipitation (+21.5%) in the upcoming years, with a reduction in the number of forests in water-limited areas and an increase in the percentage of tree cover in almost all forests.

1. Introduction

In recent years, a growing awareness of climate changes and related consequences has developed, and many countries are promoting mitigation and adaptation actions. However, despite this, climate changes are still taking place at an increasing rate across the globe [1,2] In the last few decades, Mediterranean countries, for example, have experienced increasing drought frequency and severity [3,4,5,6,7,8]. These findings are in line with the comparative analysis by Giorgi (2006) [9] that places the Mediterranean Basin among the most responsive regions to global climate change [10]. Future climate predictions are not encouraging, since they foresee a further air-temperature increase, particularly in the summer season, and a precipitation decrease, especially in winter [11,12].

Differences in temperature or precipitation patterns have produced a direct impact on the distribution of the terrestrial ecoregions and biomes [13]. The Mediterranean vegetation biome, for example, is characterized by the presence of shrubs (maquis) and forests, and their importance is first and foremost linked to the goods and services they provide. The Mediterranean Basin is the world’s second largest biodiversity hotspot; it is rich in plant diversity [14,15]; and it is recognized that its forests make a vital contribution to rural development, to poverty alleviation, and to food security, as well as to carbon up-take and accumulation, to controlling floods and erosion, and to sustaining biodiversity [14]. Although it has been recognized that Mediterranean forests are affected by climate change [16,17], their response to current and projected climate changes is not yet understood [1].

In many parts of the world, forested ecosystems are being transformed by both climate and land-use changes [18]. The impact of the former is more difficult to estimate when both coexist [19]. The relatively low urbanization and human activity makes the island of Sardinia, located in the middle of the Western Mediterranean Basin, an excellent reference condition for hydrological studies on past and future climate change [20,21,22,23]. Approximately half of Sardinia is covered by broadleaf forests and Mediterranean shrubs and garrigue [24], and it is the Italian region with the largest forest area. Furthermore, an extended and robust hydrological database is available in Sardinia, with almost one century of data, providing the opportunity to investigate climate change trends and their effects on forest ecosystems.

Regarding the air temperature in the Mediterranean Basin, a clear warming trend over the last four decades, for both annual and seasonal averages, is observed [25,26]. Based on the MAR1 report [27], the temperatures over land are projected to rise in the range of 0.9 to 1.5 °C or 3.7 to 5.6 °C during the 21st century for low (RCP2.6) or high greenhouse gas emissions (RCP 8.5), respectively, exceeding global mean temperature changes by 20% on an annual basis and 50% in summer, and with an increasing frequency and intensity of hot extremes [10,28]. Caloiero and Guagliardi (2020) [29] estimated an increasing trend of the Sardinian maximum monthly temperature, in the spring and summer months, for the period 1982–2011. Regarding precipitation, changes in the amounts and timing of rainfall events have been estimated in the Mediterranean regions. For instance, Feidas et al. (2007) [30], over the Greek territory, estimated a significant downward trend in annual precipitation determined by a reduction in winter precipitation, as confirmed more recently by Varlas et al. (2022) [31]. In Serbia and Montenegro, Tošić (2004) [32] analyzed the spatial and temporal variability of winter and summer precipitation at 30 rainfall stations for the period 1951–2000, finding again a decreasing trend in winter precipitation. A decreasing trend of annual precipitation was also found in Turkey [33], in the Northern Syrian Coastal Region [34], in Egypt [35], in Algeria, in Tunisia, and in Morocco [36]. A general negative trend in winter and autumn precipitation, as opposed to a positive one in spring and summer, has been found in Italy’s southern regions [37]. Montaldo and Sarigu (2017) [23] demonstrated that, in the period 1922–2011, winter precipitation in Sardinia was decreasing alarmingly, with a strong link between the winter rainfall trend and the North Atlantic Oscillation (NAO) index. The decreasing winter precipitation produced a drastic reduction of runoff, which diminished by more than 40% between 1975 and 2010 compared to the previous 1922–1974 period [23], confirming the sensitivity of island of Sardinia to climate-change effects.

Future climate scenarios by Global Climate Models (GCMs) are predicting an increase of droughts and warmer conditions in the Mediterranean Basin [11] and on the island of Sardinia [21].

From an ecological point of view, the variation in precipitation regimes and the increasing frequency of drought influence the coexistence of different biomes. Forests and savannas, for example, may be considered alternative biome states [38]. Sankaran et al. (2005) [39] found that, in Africa, a MAP of 650–700 mm is a threshold for woody cover density, which linearly increases with MAP up to 650 mm, then ceasing to be related to rainfall increase above 650, and then ceasing to 700 mm [40]. Most of the studies on this topic are focused on Africa, America, and Australia [38,39,41,42,43], while analogous research in Mediterranean regions, with a similar precipitation gradient, has not yet been seen.

Trees do not all respond in the same way to changing climates and environments, and some species are naturally resistant and resilient to climate-influenced disturbances, whereas others are sensitive and vulnerable [44,45,46]. The importance of past climate-change conditions in forest ecosystems is still under investigation [28,44,45,47], and accumulating evidence indicates that warmer temperatures and the increasing vapor-pressure deficit (VPD) are raising the mortality risks for forests during drought [48] and are negatively impacting vegetation growth [49]. Recently, Gazol and Camarero (2022) analyzed the compound effect of high summer VPD and low soil moisture on defoliation and tree mortality across European forests, highlighting that the steady increase in VPD is raising the risk of increased drought-triggered mortality rates of forests in Southern and Eastern Europe particularly.

Nowadays, long time series of optical remote-sensor observations that allow the monitoring of tree-cover dynamics are starting to be available [50,51], although not as extended as the hydrological database, and can be an appropriate support for climate-change studies. Although the observations of the Advanced Very High Resolution Radiometer (AVHRR) optical sensor started early (from 1984), they are at a coarse spatial resolution (1.1 km) that is inappropriate for heterogeneous Mediterranean ecosystems. Instead, the observations of the MODIS sensor on the Aqua and Terra satellite platforms are robust and still extended enough (from 2000), and with a more appropriate spatial resolution (250 m) [9,52,53]. In particular, a MODIS product, the MOD44B Version 6 Vegetation Continuous Fields—Collection 6 provides robust estimates of yearly tree-cover dynamics [43,54,55,56,57,58,59]

In the last decades, Sardinia has experienced long droughts that impacted the forest cover and tree density in several forests (see Supplementary Figures S2 and S3), showing alarming vulnerable conditions and raising the question of the actual and future hydrologic sustainability of the Sardinian forest ecosystems. Given this, and considering that it is well documented that MAP has a downward trend [60,61], while air temperature has a positive trend [27], in the Mediterranean Basin, we investigate the effects of climate change on the hydrological sustainability of the Sardinian forests. The specific objectives are to (1) quantify past trends in the amount, type, and distribution of tree cover in the main Sardinian forests; (2) evaluate the historical trends of the precipitation, air temperature, and vapor-pressure deficit in the Sardinian forests; and (3) evaluate the impacts of past and also future changing climate conditions on the Sardinian forest tree cover. The outcome of this study is expected to identify the key climate factors controlling tree growth and mortality in the main Sardinian forests, as well as the responses of different tree species, and therefore provide suggestions for forest management practice.

2. Methodology

2.1. Study Area

The island of Sardinia is a region of Italy located between 38°510′N and 41°150′N latitude and 8°80′E and 9°500′E longitude, with a surface of about 24,100 km², and is second in size only to Sicily among the islands of the Mediterranean Basin (Figure 1). Sardinia has an ancient geoformation and is not earthquake prone. Its granitic basement rocks date, in fact, from the Palaeozoic and the pre-Palaeozoic Era, and they are overlain by discontinuous sedimentary and volcanic units deposited in marine or continental environments between the Permo-Carboniferous and the Quaternary [23,62]. Due to long erosion processes, the island’s highlands, formed of granite, schist, trachyte, basalt, sandstone, and dolomite limestone, range between 300 m and 1000 m. The highest peak is Punta La Marmora (1834 m), part of the Gennargentu mountain range in the center of the island. Although the island’s ranges and plateau are separated by wide alluvial valleys and flatlands (e.g., Campidano, Nurra, etc.), and although the island has an average elevation of 334 m, the slopes are so steep that most of the island is considered mountainous [62]. The climate is maritime Mediterranean, with mild, rainy winters contrasting with hot, dry summers. Two macrobioclimates (Mediterranean pluviseasonal oceanic and Temperate oceanic), one macrobioclimatic variant (Sub-Mediterranean), and four classes of continentality (from weak semihyperoceanic to weak semicontinental) can be found in Sardinia [63].

2.2. Tree Cover of the Main Sardinian Forests

Sardinia’s wildlands were selected based on the Corine Land Cover (CLC) classification produced within the framework of the Copernicus Land Monitoring Service, referring to the land cover/land use status of years 2018 and 2000 (the source dataset is described and accessible from https://land.copernicus.eu/pan-european/corine-land-cover (accessed on 1 September 2022)). This database relies on standard methodology and nomenclature with the following base parameters: 44 classes in the hierarchical 3-level CLC nomenclature, a minimum mapping unit for status layers of 25 hectares, and a minimum width of linear elements of 100 m. Among all the available classes, only the broad-leaved forests (Class ID of CLC: 3.1.1 broad-leaved forests, BLFs), the coniferous forests (Class ID of CLC 3.1.2 coniferous forests, CFs), and the bushy sclerophyllous vegetation (Class ID of CLC 3.2.3 sclerophyllous vegetation, SV) were used in this study (

Table 1. Corine Land Cover classification and description of the land cover analyzed by this study.

CLC Class ID	Description
3.1.1.	Broad-leaved forest (BLF) Vegetation formation composed principally of trees, including shrub and bush understory, where broad-leaved species predominate. This class is applicable for mature forests of natural or anthropogenic origin. In Sardinia, it mainly includes evergreen broad-leaved woodlands composed of sclerophyllous trees (mainly Quercus Ilex, Quercus suber, Quercus Rotondifolia). Broad-leaved area must cover at least 75% of the forest tree component.
3.1.2.	Coniferous forest (CF) Woodlands consisting mainly of trees, but also of bushes and shrubs, in which coniferous forest species dominate. This class is applicable for mature coniferous (needle-leaved) forests of natural or anthropogenic origin pure or mixed stands of fir (Abies), pine (Pinus), spruce (Picea), cedar (Cedrus), cypress (Cupressus), juniper (Juniperus), yew (Taxus), and species, among others. The coniferous surface must constitute at least 75% of the forest tree component.
3.2.3.	Bushy sclerophyllous vegetation (BSV) Bushy sclerophyllous vegetation in a climax stage of development, including maquis, matorral, garrigue, and thermo-Mediterranean brushes, characterized by shrubs of sclerophyllous character (Quercus, P. latifolia, Pistacia lentiscus, Arbutus unedo, Mirtus, Juniperus, Cistus, Lavandula, Rosmarinus, Euphorbia, Genista, Erica, Laurus, etc., species), with hard leaves and short internodes, usually evergreen.

). These classes include only land-cover types that have not been significantly modified by human activity, and they exclude urbanized and agriculture lands. Table 1 describes the main properties of the three land-cover classes.

Tree cover was estimated from satellite information, using a product of the widely used optical MODIS sensor, the MOD44B Version 6 Vegetation Continuous Fields—Collection 6 [42,53,54,55,56,57,58]. MOD44B provides a yearly global representation of the percentage of three ground-cover components at a 250 m spatial resolution: tree-canopy coverage (TC), non-tree vegetation (NTV), and non-vegetated (bare) coverage (BS). Moreover, trees are defined as woody vegetation higher than 5 m. Data for the 2000–2019 period were used. CLC information for the years 2000 and 2018 was resampled at the 250 m spatial resolution of the MODIS product. In this study, we did not consider any pixels falling into one of the following categories: (i) nonhomogeneous land cover, (ii) land-cover change between CLC 2000 and CLC 2018, and (iii) affected by forest fires between 2000 and 2020. We selected 17 forests (Figure 1) that have been stable since 1932 [64] and included mainly tall trees with a sparse understory layer of Q. ilex subsp. Ilex, dense forests with a garrigue or maquis formation, and open woodlands of Q. ilex subsp. ballota (frequently described as savannah-like ecosystems, with low canopy cover, accompanied by grass or shrub understory).

Table 2. Properties of the 17 Sardinian forests indicated in Figure 1: code; name; surface; mean altitude; and percentage of the surface covered by broad-leaved vegetation (BLF), bushy sclerophyllous vegetation (BSV), coniferous forest (CF), and other vegetation.

Code	Name	Surface	Mean Altitude	Surface Covered by BLF	Surface Covered by BSV	Surface Covered by CF	Surface Covered by Other Vegetation
	(-)	(km2)	(m asl)	(%)	(%)	(%)	(%)
F1	Barigadu	36.23	639.8	43%	35%	5%	17%
F2	Sos Littos—Sas Tumbas (Bitti)	19.27	378.2	73%	7%	2%	18%
F3	Monti del Gennargentu	447.53	1047.3	24%	8%	3%	65%
F4	Supramonte di Oliena, Orgosolo e Urzulei—Su Sercone	235.08	765.9	26%	17%	2%	54%
F5	Foresta di Monte Arcosu	303.8	449.6	52%	24%	-	24%
F6	Monte dei Sette Fratelli e Sarrabus	92.98	531.5	72%	11%	-	17%
F7	Monte Linas—Marganai	236.46	551.8	26%	28%	-	46%
F8	Monte Arci	101.82	414	32%	35%	5%	28%
F9	Monte Ferru	32.05	336.7	28%	36%	0%	37%
F10	Villasalto	22.79	490.5	2%	73%	15%	10%
F11	Fiorentini	15.88	796.4	85%	-	-	15%
F12	Monte Lerno	99.26	612.9	22%	26%	25%	27%
F13	Sorilis	49.92	548.2	-	50%	12%	38%
F14	Siamanna	37.93	343	0%	56%	13%	31%
F15	Tacchixeddu	29.07	461.5	19%	55%	6%	20%
F16	Su Sartu	11.03	555.5	-	41%	-	59%
F17	Nuoro	14.2	630.5	54%	16%	-	30%

summarizes the main properties of the forests under investigation. The forest areas range from 11.03 km² of Su Sartu (F16) to the 447.5 km² of the Monti del Gennargentu (F3), which is also the forest at the highest average altitude (1147 m asl).

2.3. Meteorologic Data

The rain-gauge network over Sardinia consists of 414 stations [23], and data are available for almost one century (1922–2019). We further selected 179 rain gages with more than 80 complete years of data. In more detail, 51.4% of them have ≥80 complete years of data, 86% have ≥70 complete years of data, 91.62% have ≥60 complete years of data, and 100% have more than 59 years of data.

Daily maximum and minimum air temperatures at 112 stations are available for the period 1981–2010 (39 years), but only a few of them have remained available in the last ten years, for which tree-cover data are instead available. To overcome the reduced spatial availability and density of the observed temperature data, and to extend the data set to 2019, the ERA5 monthly reanalysis dataset for the 1980–2019 period for both air temperature and dew point was used. We successfully tested the accuracy of the reanalysis dataset (correlation coefficient > 0.62, RMSE < 1.14 °C) in 6 of the observing Sardinian stations, with almost-complete datasets (Carbonia C.ra Flumentepido, Tertenia, Senorbi, Macomer, Villa Verde, Alghero, Italy).

VPD was calculated by using the ERA5 monthly near-surface air temperature (T) and dew point (Td), using the following equation [65]:

$$VPD=c_1\ast \exp (\frac{c_2\ast T}{c_3+T})-c_1\ast \exp (\frac{c_2\ast Td}{c_3+Td})$$

(1)

where c₁ = 0.611 kPa, c₂ = 17.5, and c₃ = 240.978 °C. Td and T are in °C, and the resulting VPD is in kPa.

2.4. Statistical Analysis

The mean air temperature, VPD, and precipitation data were analyzed at the annual (T_y, VPD_y, and P_y) and seasonal scales. We used the Sardinian hydrologic year that begins on October 1 (when runoff usually resumes) and ends on September 30 of the following year (when the dry season usually ends) [21]. We distinguished the seasons as (i) ‘extended winter period’ (December through March, T_ew, VPD_ew, and P_ew); (ii) ‘shortened spring period’ (April and May, T_ss, VPD_ss, and P_ss); (iii) ‘extended summer period’ (June through September, T_es, VPD_es, and P_es); and (iv) the ‘shortened fall period’ (October and November, T_sf, VPD_sf, and P_sf). Note that the use of the extended seasons is common in climate-change studies [21,23,66].

Trends were estimated by using the Mann–Kendall non-parametric test [67], while the slope of the linear trend was estimated with the Theil–Sen method [68,69], both of which are widely used in hydrologic applications [23,70,71,72]. To reduce the impact of autocorrelation on the statistical significance of the Mann–Kendall test and on the Theil–Sen method, different prewhitening methods are proposed by the literature [73,74,75,76]. The algorithm, proposed by Collaud et al. (2020), was used here.

The Pettitt test, a method widely used in hydrology and climate analysis [23,77,78], was used to detect the presence of a change point in the time series of the average annual precipitation over Sardinia. The spatial distribution of the MAP values observed in the rain-gauge stations over the Sardinia region was estimated through an ordinary kriging model (linear semi-variogram, 250 m resolution, and 16 neighbors). Finally, the Pearson correlation coefficient, ρ, was used for estimating the linear correlation between two datasets.

2.5. Future Projections

For future climate scenarios of precipitation and air temperature, we considered the HadGEM2-AO climate model [79] predictions from the Coupled Model Intercomparison Project Phase 5 (CMIP5), and a model coming from the following generation of HadGEM configuration, the HadGEM3-GC31-LL [80,81], from the Phase 6 (CMIP6). Montaldo and Oren (2018) showed that the HadGEM2-AO model is the one that best represents the past variations in precipitation and air temperature in Sardinia. Future climate scenarios up to 2090 were used to predict the effect of future climate change on tree cover in the forests. Following Montaldo and Oren (2018), we considered two main representative concentration pathways (RCPs) of the HadGEM2-AO model (atmospheric horizontal resolution of 1.875° × 1.25°), i.e., RCP 45 and RCP 85. In addition, two Shared Socioeconomic Pathways (SSPs) of the HadGEM3-GC31-LL model (atmospheric horizontal resolution of 1.875° × 1.25°) were considered, i.e., SSP 245 and SSP 585). The multivariate bias-correction algorithm [82], which transfers all statistical characteristics of the land-observed meteorological data to the corresponding multivariate distribution of variables simulated by the GCMs, was applied for the downscaling and bias correcting of both HadGEM2-AO and HadGEM3-GC31-LL for Sardinian forests.

3. Results

3.1. Long-Term Trends in Precipitation

Rainfall is not abundant, and on the coast, the MAP typically ranges from 415 mm (Cagliari station at 7 m asl) to 934 mm (Baunei C.ra station at 480 m asl) mm per year (Figure 2a) and follows the Mediterranean pattern, so that rainfall is concentrated in the autumn (mean value of 180 mm) and winter (mean value of 305 mm), gradually decreases during the spring (mean value of 94 mm), and hits a low in the summer (mean value of 74 mm). In inland areas, rainfall locally exceeds 700 mm per year in hilly areas and 1000 mm in the mountains (Figure 2a). Considering the full rainfall dataset (1922–2019), the mean of the island’s average annual rainfall is 729.08 mm (SD of 143.8 mm), and it varied between 469.5 mm (year 1999) and 1062 mm (year 1933). Nine of the top ten driest years that took place in this century occurred after 1980. Pettitt’s test detects a changing point in the average annual precipitation over the island in 1980. A total of 51% of the observing stations showed a significant negative trend in annual precipitation, but the spatial variability is high (Figure 3a). Limiting the rainfall dataset to the stations within close proximity of the forests analyzed, it emerges that the annual precipitation decreases significantly in almost all forests, except for F3, F4, F11, and F17, and that the β_Py < −3 mm per decade in the rainfall stations is close to F5 (−3.6 mm/y at Is Cannoneris Caserma), F7 (−3.72 mm/y at Monti Mannu Caserma), F8 (−3.1 mm/y at Marrubiu, and −3.31 at Baradili), and F12 and F13 (−3.12 mm/y at Ala dei Sardi). Seasonally, precipitation decreases were largest in the winter (Figure 3b), while little and no significant changes are observed for the other seasons. The island’s average winter precipitation (Pew) was 344.47 mm (standard deviation, SD, of 98.18 mm), and it varied between 127.66 mm (year 1999) and 546.14 mm (year 1927). Even in the extended winter average precipitation, we can observe an important spatial variability, with β_Pew < −25 mm per decades in the rainfall stations within close proximity of F5 (−3.03 mm/y at Is Cannoneris Caserma), F7 (−2.70 mm/y at Monti Mannu Caserma), F8 (−2.43 mm/y at Marrubiu), and F13 (−2.41 mm/y at Ala dei Sardi). Variations in the Theil–Sen slope estimator of precipitation during the short spring (β_Pss) and the extended summer (β_Pes) across the island are in the range ± 0.65 mm/y, and ±1.00 mm/y for the short fall (β_Psf).

The past change of the annual precipitation was well depicted in Figure 2, where the spatial variability of MAP in the island during the whole observation period (Figure 2a) and before and after the change point in 1980 was compared (Figure 2b,c). By comparing the spatial distribution of MAP in the 1922–1979 period (MAP_1922–1979, 57 years; Figure 2b) with the spatial distribution of MAP in the 1980–2019 period (MAP_1980–2019, 39 years; Figure 2c), it emerges that MAP is decreasing drastically in many parts of the island. In the period 1922–1979, areas with a MAP value of less than 700 mm occupied 26% of the island’s surface. This value rises to 62% in the period 1980–2019 and even includes some areas such as those in the southern part of the island and the Campidano Plain, with a MAP lower than 450 mm. At the 17 selected forests, we observed that MAP is decreasing everywhere by an average of −122.86 mm. F13 and F14 are the forests that recorded the highest (−243 mm) and the lowest (−84 mm) reduction values, respectively.

Looking at a more recent period (1990–2019, 28 hydrological years; Figure 2d)), aridity increased even more. The Theil–Sen estimator applied to the annual precipitation is negative for 54 stations, peaking above 5 mm/y for 11 stations and even above 8 mm/y for 4 stations (Marrubiu C.ra (β_Py = −12.9 mm/y), Lanusei (β_Py = −10.5 mm/y), Thiesi (β_Py = −9.7 mm/y), and Perdasdefogu (β_Py = −8.7 mm/y)). The spatial distribution of MAP in the 1990–2019 period (MAP_1990–2019, 28 years; Figure 2d) showed that areas with less than 700 mm amounted to 78% of the island’s surface. In regard to the forests, two have a MAP_1990–2019 <600 mm (F8 and F14), seven have a 600 mm ≤ MAP1990–2019 < 700 mm (F5, F6, F10, F11, F12, and F17), seven have a 700 mm ≤ MAP1990–2019 < 800 mm (F1, F2, F4, F7, F9, F13, and F15), and only F3 has a MAP_1990–2019 of about 900 mm/y.

3.2. Trends in Tree Cover

A total of 155,518 MODIS pixels were sampled across the study area for monitoring Sardinian forests’ tree cover (Figure 1), corresponding to about 28% of the entire island. The dominant class (63.59%) is bushy sclerophyllous vegetation (BSV), while broad-leaved vegetation (BLF) accounts for 30.97%, and coniferous vegetation (CF) for just 4.51%. Regarding the 17 forests (Figure 1), in five of them (F22, F21, F20, F6, and F5) the percentage of broad-leaved vegetation is >50%, while seven forests (F23, F17, F19, F18, F16, F14, and F13) do not contain any broad-leaved vegetation (Figure 1). Bushy-sclerophyllous vegetation is present at all sites, and it dominates the other land uses in F11 (93%), F10 (73%), F14 (56%), F15 (55%), and F13 (50%). Coniferous vegetation is present in 13 of the 17 forests, but only in F12 and F18 does it exceed 20%. Due to the small number and the sparse distribution of coniferous compared to other vegetation categories, CF is not further analyzed in this paper.

The long-term trends of the percentage of tree-canopy coverage, non-tree vegetation, and bare soil over the period 2000–2019 were calculated based on MOD44B data. Figure 4 provides the spatial distribution of the Mann-Kendall τ values for both the BLF (Figure 4a–c) and BSV (Figure 4d–f) vegetation types, while

Table 3. The mean Mann-Kendall τ and the Theil-Sen slope β of the tree canopy coverage (TC), non-tree vegetation (NTV), and non-vegetated (bare) coverage (BS) for all the broad-leaved forests (BLFs) and bushy sclerophyllous vegetation (BSV) observed pixels falling within each of the 17 sites analyzed.

Code	τTCBLF	τTCBSV	τNTVBLF	τNTVBSV	τBSBLF	τBSBSV	βTCBLF	βTCBSV	βNTVBLF	βNTVBSV	βBSBLF	βBSBSV
	(-)	(-)	(-)	(-)	(-)	(-)	(%TC/y)	(%TC/y)	(%NTC/y)	(%NTC/y)	(%BS/y)	(%BS/y)
F1	−0.15	0.06	0.08	−0.07	0.10	0.01	−0.16	0.14	0.07	−0.16	0.08	0.01
F2	−0.12	0.04	0.10	−0.03	0.02	−0.07	−0.15	0.09	0.14	−0.07	0.07	−0.02
F3	0.03	−0.11	0.02	0.05	−0.12	−0.24	0.12	−0.24	0.05	0.35	−0.14	−0.20
F4	0.00	0.06	0.06	−0.05	−0.14	−0.02	0.05	0.12	0.11	−0.12	−0.14	0.01
F5	−0.30	0.09	0.18	0.00	0.00	−0.16	−0.30	0.20	0.30	−0.03	0.00	−0.09
F6	−0.12	0.11	0.03	−0.13	0.11	0.11	−0.09	0.19	0.02	−0.29	0.10	0.06
F7	−0.10	0.22	0.03	−0.15	0.11	−0.19	−0.11	0.35	0.01	−0.27	0.09	−0.09
F8	−0.19	0.14	0.06	−0.15	0.25	0.05	−0.21	0.20	0.05	−0.32	0.34	0.13
F9	−0.14	0.27	0.02	−0.19	0.17	−0.07	−0.15	0.50	0.01	−0.41	0.10	−0.01
F10	−0.21	−0.03	-	−0.01	0.08	0.04	−0.19	−0.06	NaN	−0.01	0.02	0.07
F11	−0.09	0.09	0.05	0.02	−0.03	−0.15	−0.09	0.17	0.09	−0.01	0.01	−0.15
F12	−0.09	0.15	0.03	0.12	0.06	−0.30	−0.10	0.22	0.04	0.17	0.04	−0.35
F13	-	0.01	-	−0.01	-	0.01	-	0.01	-	−0.04	-	0.01
F14	-	0.07	-	−0.10	-	0.08	-	0.13	-	−0.20	-	0.05
F15	−0.29	0.08	0.19	−0.07	0.10	−0.03	−0.22	0.14	0.22	−0.16	0.06	0.01
F16	-	0.01	-	−0.04	-	0.10	-	0.07	-	−0.13	-	0.07
F17	−0.07	0.04	0.03	−0.07	0.07	0.02	−0.05	0.11	0.02	−0.14	0.04	0.04

provides the average Mann-Kendall τ and the Theil-Sen slope β values of the tree-canopy coverage, non-tree vegetation, and bare soil within each of the 17 forests.

Broad-leaved forests with significant negative trends in the Mann-Kendall τ values for tree-canopy coverage are clustered mainly in the forested areas covered by this study and particularly in F5, F8, F10, and F15. Ten of the thirteen forests that contain broad-leaved vegetation show a Theil-Sen slope β ranging between −0.30% TC per year and +0.12% TC per year. Forest F5, which is the second largest forest on the island, is the one with the lowest value of Mann-Kendall τ for the tree cover canopy (τTC), equal to −0.3 (Theil-Sen slope β = −0.38%TC per year). In this forest, as in F15 the decrease in τTC contrasts with an increase in the τ of the vegetation of less than 5 m (τNTV_BLF), while, for example, in F8 and F9, there was a growth in the value of bare soil (τBS). The Mann-Kendall τ of the tree-cover canopy of forests F3 and F4, which are in high elevation mountains in the central east of Sardinia, is, instead, not significant and moderately positive. The spatial distribution of trends in bushy sclerophyllous vegetation (Figure 4d–f) highlights an overall increase in treecover coverage, and a relative reduction in non-tree vegetation and non-vegetated areas. Bushy sclerophyllous vegetation in the analyzed forests shows a sensitive negative trend in τTC only in F3, and, on the contrary, 15 sites show a positive trend in τTC_BSV, with βTC_BSV ranging from +0.01 to +0.5% per year.

The relationship between tree canopy coverage and MAP_1922–2019 for BLF and BSV was investigated year by year. As an example, Figure 5 shows the relationship of MAP_1922–2019 versus tree-canopy coverage for the year 2008 for both BLF and BSV. We evaluated the breakpoint (MAP_cp, identified by the regression analysis in Figure 5) at which the maximum tree cover is attained. The change in tree-canopy coverage as a function of MAP has been determined by using a 99th quantile piecewise linear regression with the identification of the breakpoint (the rainfall at which the maximum tree cover is attained), and of the maximum upper bound on tree cover (TC_max). We showed the scatter plot for the year 2008 in Figure 5 as an example because it was the year where the changing point was closest to the average value. The value of MAP_cp during the 2000–2019 period was, on average, 695.3 mm (SD of 51.65 mm) for broad-leaved vegetation (between 607.47 mm (year 2002) and 759 mm (year 2019)) and 736.76 mm (SD of 48.30 mm) for bushy sclerophyllous vegetation (between 630.00 mm in the year 2004, and 810 mm in the year 2014). The maximum percentage of tree-cover canopy (TC_max) ranges from 67.08% (year 2017) to 80.30% (year 2013), with an SD of 3.38%, for broad-leaved pixels; and from 56.75 (year 2003) to 69.84% (year 2013), with a SD of 3.53%, for bushy-sclerophyllous pixels (Figure 6). There is a significant negative trend of TC_max for broad-leaved vegetation (τ = −0.47, p-value < 0.05), while for bushy sclerophyllous vegetation, the trend is not significant (τ = −0.02 and p-value = 0.92).

3.3. Long-Term Trends in Air Temperature and Vapor-Pressure Deficit

The average yearly air temperatures (Ty) in Sardinia have increased significantly by nearly 0.38 °C per decade (τ = +0.33 with p-value < 0.05) since 1990. Warming has been occurring significantly during all seasons but has been greatest during the short spring (τ = +0.33 with p-value < 0.05) and extended summer seasons (τ = +0.40 with p-value < 0.05), with a magnitude of increase of 0.48 and 0.51 °C per decade, respectively. The spatial variability of the Mann-Kendall τ values of the annual air temperature (τTy) on the island is displayed in Figure 7a. It shows an increasing trend over the entire island, with τ values ranging from 0.34 in the hilly areas of the island up to more than 0.5 in the northern and eastern part of the island and particularly close to the coastline. The Theil-Sen slope estimator (Figure 7b) reveals that Sardinia is warming on average by 0.4 °C per decade, with a gradient of increase from west (on average 0.035 °C/y) to east (on average 0.04 °C/y). Focusing on the 17 forests analyzed, all of them are subject to increasing annual and seasonal average temperatures. The maximum rises in average annual temperature are recorded in the forests located in the eastern side of the island (F9, F15, and F6, with a Theil-Sen slope estimator value of 0.52, 0.49, and 0.44 °C per decade, respectively), while by contrast, those located on the central western part of the island recorded the lowest rises (F2, F1, and F11, with a Theil-Sen slope estimator value of 0.52, 0.44, and 0.44 °C per decade, respectively). The Theil-Sen slopes (β) of the trend lines of air temperature during the short spring and extended summer are positive for all the forests, and the rise in T_ss and T_es is greater than 0.05 °C/y in 15 and 7 forests, respectively.

Figure 8a,b show the spatial variability of the Theil-Sen slope estimator β on the island of the annual VPD (VPD_y) and in the ‘extended summer season’ (VPD_es). Overall, an increase of VPD was observed for most of the study area (Figure 8a). Increases in annual VPD have risen at a rate of approximately 0.003–0.007 kPa/y for 75% of the island, and this extrapolates to a total of 0.12 kPa for the 30-year period. The largest variations in annual VPD were in the western and northern portions of the island and the lower midsouth, while the lowest ones were in the inland hilly areas. During the extended summer season, the VPD increased over the largest areal extent of land, for which, 80% of the surface showed a Theil-Sen slope estimator in the range of 0.010–0.017 kPa/y. The most significant increases were in the mid-northern part of the island. VPD is increasing in all 17 of the forests analyzed, both annually and in every season. Increases in annual VPD have risen at a rate of approximately 0.005–0.009 kPa/y in seven forests, while for the remainder, the rate of growth is slightly lower.

3.4. Evaluation of the Climate Impact on Forest Cover Changes

We found interesting relationships between the trends of the forest tree cover and some climate drivers, in particular, the MAP_1990–2019 and the Theil-Sen slope estimator of the average annual VPD (Figure 9). Correlation coefficients (

Table 4. Pearson correlation coefficients (ρ), and relative p-value between the Theil-Sen slope estimators of the percentage of tree cover and bare soil in broad–leaved (βTC_BLF and βBS_BLF) and in bushy sclerophyllous vegetation (βTC_BSV and βBS_BSV) with some of the main climate variables.

	βTCBLF		βTCBSV		βBSBLF		βBSBSV
	ρ	p-val	ρ	p-val	ρ	p-val	ρ	p-val
MAP1990–2019	0.71	0.004	0.13	0.46	−0.68	0.01	−0.44	0.07
βTes	−0.22	0.45	−0.13	0.61	0.02	0.24	0.18	0.49
βTew	−0.22	0.45	−0.09	0.73	0.05	0.70	0.14	0.60
βTsf	−0.27	0.35	−0.01	0.98	−0.03	0.35	0.13	0.62
βTss	0.04	0.91	0.16	0.53	0.00	0.89	−0.23	0.38
βTy	−0.33	0.25	−0.07	0.78	−0.09	0.25	0.35	0.17
βVPDes	−0.12	0.68	−0.42	0.09	0.16	0.52	0.27	0.30
βVPDew	−0.31	0.28	0.01	0.97	−0.13	0.39	0.26	0.32
βVPDsf	−0.32	0.26	−0.34	0.18	0.04	0.08	0.32	0.21
βVPDss	−0.28	0.33	−0.39	0.12	0.00	0.10	0.46	0.06
βVPDy	−0.61	0.02	−0.34	0.18	−0.62	0.02	0.44	0.08

) whose magnitudes are between 0.5 and 0.7 indicate variables that can be considered moderately correlated, while correlation coefficients whose magnitudes are between 0.3 and 0.5 indicate variables that can be considered only slightly correlated. The Theil-Sen slope estimator of the percentage of tree cover of broad-leaved vegetation (βTC_BLF) was found to be positively correlated with MAP_1990–2019 (ρ = 0.71, p-value 0.004) and negatively with β_VPDy (ρ = −0.61, p-value 0.02). Opposite to this is the behavior of the Theil-Sen slope estimator of the percentage of bare soil in the analyzed broad-leaved vegetation (βBS_BLF), which can be considered a measure of the forest mortality rate. βBS_BLF shows a positive correlation with the annual change of annual VPD (β_VPDy, ρ = 0.62, p-value 0.02) and a negative correlation with the mean annual precipitation (MAP_1990–2019, ρ = −0.68, p-value 0.006). The correlations between the climatic variables analyzed and the variation in the percentage of tree cover and bare soil within the bushy sclerophyllous vegetation pixels of the forests analyzed are all low. Overall, therefore, the broad-leaved vegetation seems to be more affected by a variation in climate drivers than the bushy sclerophyllous vegetation. Figure 9 shows the scatterplots between the Theil–Sen slope estimators of the tree canopy (βTC_BLF) and of the bare soil (βBS_BLF) coverages and of the annual VPD(β_VPDy) and MAP_1990–2019 climate drivers. We fit (a) a linear function to describe the variability of βTC_BLF to β_VPDy (βTC_BLF = −41.60 β_VPDy + 0.092 (2), R² = 0.37; Figure 8a); (b) a sigmoidal function to describe the relationship between βTC_BLF and MAP_1990–2019 (βTC_BLF = −0.162+ 0.283/(1 + exp(−(MAP_1990–2019- 750.68)/14.985); R² = 0.81; Figure 8b); (c) a polynomial function to describe the variability of βBS_BLF to β_VPDy (βBS_BLF= −0.136+ 24.82*β_VPDy + 2124.17*β_VPDy²; R² = 0.62; Figure 8c); and (d) a linear function to describe the relationship between βBS_BLF and MAP_1990–2019 (βBS_BLF = −0.00095 MAP_1990–2019 + 0.706; R² = 0.47; Figure 8d). All of these four relationships are significant highlighting that the warmer conditions, expressed by both an increase of vapor-pressure deficit and a decrease of precipitation, are affecting tree cover in the Sardinian forests. In Figure 9a, most of the data are close to the βTC_BLF (β_VPDy) linear regression line, except the F3 data, for which βTC_BLF is moderately positive (but not significantly), and β_VPDy is positive. Note that the F3 climate conditions are different from the other forests’ conditions because F3 is the forest at the highest elevation (Figure 1), with the lowest annual mean air temperature (~13 °C) and vapor-pressure deficit (~1.7 kPa) and the highest MAP (~900 mm). We used a sigmoidal function for βTC_BLF (MAP_1990–2019) relationship for fitting the F3 and F4 βTC_BLF positive values, clearly related to the particularly wet conditions (MAP > 750 mm).

We used the estimated βTC_BLF(β_VPDy) and βTC_BLF (MAP_1990–2019) relationships of Figure 9 in the past for making a coarse prediction of the effects of future scenarios of precipitation and air temperature on forest tree covers.

Using the six future scenarios of the HadGEM2-AO of the CMIP5 and the six future scenarios of the HadGEM3-GC31-LL of the CMIP6 (

Table 5. The twelve future scenarios of the CMIP5 and CMIP6 [79,81].

Future Scenario	Coupled Model Intercomparison Project Phase	Representative Concentration Pathways	Shared Socioeconomic Pathways	Period
SC1	CMIP5	RCP 45		2020–2030
SC2	CMIP5	RCP 85	-	2020–2030
SC3	CMIP5	RCP 45	-	2055–2075
SC4	CMIP5	RCP 85	-	2055–2075
SC5	CMIP5	RCP 45	-	2081–2090
SC6	CMIP5	RCP 85	-	2081–2090
SC7	CMIP6	-	SSP 245	2020–2030
SC8	CMIP6	-	SSP 585	2020–2030
SC9	CMIP6	-	SSP 245	2055–2075
SC10	CMIP6	-	SSP 585	2055–2075
SC11	CMIP6	-	SSP 245	2081–2090
SC12	CMIP6		SSP 585	2081–2090

), we estimated the future changes of tree canopy coverage (with respect to 2000–2019) in broad-leaved forests (Figure 10). All scenarios show a positive change in annual vapor-pressure deficit (Figure 10a), ranging from 4.6% in SC7 to 38.64% in SC12. A reduction in the percentage of tree cover is predicted by using the proposed βTC_BLF (β_VPDy) relationship for all the scenarios, ranging from 4.43% in SC7 to 38.52% in SC12 (peaking above −40% in the five forests, F3, F12, F4, F11, and F2; Figure 10c). Future changes of precipitation are predicted as both positive and negative for CMIP5 scenarios (SC1–SC6), and only negative for CMIP6 scenarios (SC7–SC12). SC6, S11, and S12 were the worst scenarios, with a reduction of 7.6%,12%, and 14.3%, which produced the maximum reduction of tree cover of 16.09%, 16.46%, and 17%, respectively (using the proposed βTC_BLF (MAP_1990–2019) relationship).

4. Discussion

Sardinian forests are subject to the effect of climate factors and are suffering, as demonstrated by the decrease of the canopy cover almost everywhere in recent decades. There is currently no tool to map the large-scale spatial variability and intensity of this phenomenon. This prompted us to investigate the possibility of using satellite data, and in fact, the analysis of the historical MODIS satellite data over Sardinia reveals the presence of hotspots of significant negative trends in the tree-cover percentage (Figure 4). These are mainly located in some of the forests analyzed by this study that are protected forests, not affected by fires during the period 2000–2020, and where the human impact is insignificant (no change in CLC between 2000 and 2018). The significant negative trends in the tree cover particularly affect the broad-leaved vegetation. A total of 86% of forests show a reduction in broad-leaved tree cover, and for seven forests, the Theil–Sen slope of the percentage of tree cover in broad-leaved vegetation is (βTC_BLF) < 0.12%/y. The minimum value occurs in the forest of Monte Arcosu (F5), with a value of −0.30%/y, while the maximum one occurs in F3, equal to 0.12%/y (Table 4). In some forests, such as F5 and F15, the decrease in βTC_BLF is contrasted by an increase in shorter vegetation, and this is confirmed by the positive values of the Theil–Sen slope of the percentage of non-tree cover in broad-leaved vegetation (βNTV_BLF), while in other forests, such as F8 and F9, there is an increase in the percentage of bare soil (βBS_BLF). These results agree with the drought-triggered dieback observed in some of the analyzed forests in the last decade, visible also, for instance, from a comparison of aerial photography of the F5 forest in 2009 and 2019 (Figure 2) and, in general, after the long drought of 2017 (Supplementary Figure S3). The behavior of the bushy sclerophyllous vegetation is the opposite: in 15 out of the 17 forests the Theil–Sen slope of the percentage of tree cover (βTC_BSV) is positive, and it reaches the maximum value in F12, F07, and F09. These results are in line with the findings of others, who reported a reduction in woody cover along all the Mediterranean Basin [83], or in specific regions such as Spain and Italy (see Gentilesca et al. (2017) [84] for a review) and in Portugal [85]. Definitely, a main reason for which the forest is suffering is the strong decrease of precipitation in the last century in Sardinia. Using the extended Sardinian rainfall database, we performed an accurate trend analysis of Sardinian rainfall that mainly highlighted a change point (1980) in rainfall after which there was an important reduction in the percentage of surfaces with an average value of rainfall of more than 700 mm, which decreases significantly from 74% in 1922–1979 to 38% in 1980–2019, and an increase from 1 to 8 (out of 17) of the number of forests that are below the same threshold. The estimated mean annual precipitation (MAP) threshold (Sankaran et al., 2008) that limits the tree cover in semi-arid forests is 700 mm. Hence, we estimated the first reason of the forest suffering, which is the drastic decrease of the rainfall in the last century in Sardinia. Starting from this new drier condition of the Sardinian forests, we analyzed the recent period, for which more data are hopefully available, and we investigated the climate factors that are affecting forests the most. The analysis of the 20 years of historical MODIS satellite data over Sardinia reveals the presence of hotspots of significant negative trends in tree-canopy coverage. We analyzed the historical precipitation pattern over Sardinia, highlighting a significant average negative trend in the average annual precipitation (the average Mann-Kendall τ over the island is −0.2, with a p-value ≤ 0.05, and the Theil–Sen slope (β_Py) is −13.81 mm per decade), which is due mainly to the negative trend of winter precipitation (average Mann-Kendall τ over the island is −0.21, with a p-value < 0.005). These alterations of the rainfall regime are a widespread phenomenon in the Mediterranean Basin [23,30,31,32,33,35,37], due to climate-change effects [35]. In Sardinia, 61% of the rainfall stations showed a significant negative trend for the average annual precipitation, but the phenomenon has a marked variability over space (Figure 3a), with a decrease of the negative trend toward the mid-eastern coast (e.g., F3 and F4 forests), as anticipated by Montaldo and Sarigu (2017), who explained that this may be caused by the combined effect of the orography and of North Atlantic Oscillation (NAO) [86]. The reduction on annual total annual precipitation caused a reduction of MAP (Figure 2), and this has been recognized as one of the factors that most influence the spatial distributions of vegetation [87]. Sankaran et al. (2008) demonstrated that, when the MAP value falls below 700 mm, the upper bound of the percentage of tree cover may be limited mainly by water availability, while above this threshold, TC has little to no dependence on rainfall. The MAP analysis over Sardinia (Figure 2), performed by this study for before and after 1980 (change point in annual precipitation estimated with Petitt’s method), highlighted an important reduction in the percentage of surfaces with an average value of rainfall of more than 700 mm, which decreases significantly from 74% in 1922–1979 to 38% in 1980–2019, and an increase from 1 to 8 (out of 17) of the number of forests that are below the same threshold (Table 3). More recently, in the period 1990–2019, the reduction in precipitation is even more marked, and in fact, in this period, 78% of the island records less than 700 mm of rainfall per year. A total of 16 out of 17 of the forests were affected by a rainfall reduction, peaking below −8.9% in F14, F13, and F8, and with F14 and F8 that dropped below the threshold of 600 mm of rainfall per year.

Having established that MAP plays a key role in determining the spatial distributions of vegetation, which is decreasing drastically in most of the island, and that this causes an increase in the number of historical forests recently located in water-limited areas (MAP < 700 mm), an attempt was made to assess the presence of any possible relationship with the spatiotemporal pattern of vegetation. The tree cover versus MAP relationships, which were found by other authors [40,42,56,88] in savanna-like ecosystems, are in line with those found in the Mediterranean context proposed by this study (Figure 5). Indeed, our data indicate a strong dependence of tree cover on rainfall below the thresholds of 695.3 ± 51.65 mm and of 736.76 ± 48.3 mm for broad-leaved and bushy sclerophyllous vegetation, respectively, while above this, there is a higher probability of witnessing a woody plant encroachment. However, looking at the maximum tree cover reached by the vegetation over time (Figure 6), we noticed that the broad-leaved vegetation is showing a 3.4% per decade reduction in TC_max, while no trend was found for bushy sclerophyllous vegetation. This means that something is reducing the ability of broad-leaved vegetation to return to its original state, and this has been recognized as an important indicator of forest degradation [89,90]. Although the analysis of the data (Figure 5) suggested that precipitation is a key driver of tree cover stability, we investigated the further main climate factors that may explain the variability of the tree cover, such as global warming and the rising vapor-pressure deficit [91].

The average yearly air temperature (Ty) is increasing alarmingly in the entire island (with Mann-Kendall τ in the range 0.23–0.47) by about 0.39 °C per decade (see Figure 7), confirming climate-change-induced temperature-rise projections [92,93]. The annual mean temperature is growing in all forests; the smallest increases in βTy are found in the forests of the mountainous areas of the hinterland (F1 and F17) and along the west coast (F7 and F8), while forests located along the east coast record the greatest positive trends (F6, F9, and F15). Summer global warming is affecting the whole island (Figure 7b) and, thus, all forests (>0.4 °C per decade), and this exacerbates the deleterious effects of leaf overheating and may affect a wide number of metabolic processes, including photosynthesis, respiration, water transport, and phenology [16]. Moreover, the increasing trend in the air temperature may aggravate the effects of drought, because temperature modulates plant water use through its effects on VPD [25,94]. A positive trend in the average annual vapor-pressure deficit has been found in the entire island (average β_VPDy of about 0.05 kPa per decade), except for inland mountain areas, where the increase is moderate (i.e., F4 and F17 forests, Figure 8). The greatest increases in VPD_y (β_VPDy > 0.07 kPa/yr) affected the forests located in the proximity of the Campidano Plain (F14 and F8) and of the east coast of the island (F9). Plants respond to the rising vapor-pressure deficit regulating the gas exchange with the atmosphere through the closure of the stomata [94,95], and the stomatal sensitivity to the vapor-pressure deficit is highly variable across and within species [96]. The latter, such as Quercus ilex, which is the dominant species in the broad-leaved vegetation type (Figure 1), have a variable midday leaf water potential and keep their stomata open and a high photosynthetic rate for long periods, even in the presence of decreasing leaf water potential. This behavior might be beneficial when water is abundant, but under conditions of intense drought, it might endanger the plant [97]. The increase in the vapor-pressure deficit during the summer season is even greater in all forests (βVPD_es = 0.007–0.015 kPa/y), making the phenomenon just described more acute and impacting plant transpiration and photosynthesis [95].

To sum up, the increasing vapor-pressure deficit substantially influences vegetation growth and forest mortality [91,98], and if high VPD conditions are combined with reduced precipitation, the risk of experiencing drought stress quickly increases [98]. Sardinia is experiencing a situation of this type. In the last three decades, on the one hand we are witnessing an increase in both average temperatures and the vapor-pressure deficit; and on the other, there is an important variation in MAP caused mainly by a drastic reduction in the winter average precipitation (see Figure 2 and Figure 3), which is important for soil water recharge and is positively correlated with tree growth [20]. The combination of those climate drivers, as found by other authors [15,99,100,101,102], may be the main cause of the shift from tall tree to shorter or no vegetation (Figure 4), and to the reduced ability of the vegetation to return to its original state (Figure 6), particularly for BLF vegetation. In closed-canopy forests, Bennett et al. (2015) [103] explain that taller trees face a more challenging hydraulic environment because they are exposed to higher solar radiation and leaf-to-air vapor-pressure deficits than those in the relatively buffered understory. This is what is most likely happening in F5, which is characterized by a MAP_1922–2019 of 697 mm, which dropped to 659 mm MAP_1990–2019, and growing values of temperature (τT_y = 0.33) and vapor-pressure deficit (τVPD_y = 0.3) and where there is a decreasing trend in broad-leaved tree canopy coverage (τTC_BLF = −0.3) and a positive trend of broad-leaved non-tree vegetation coverage (τNTV_BLF = +0.18). In F8, instead, the reduction in broad-leaved tree canopy coverage (τTC_BLF = −0.19) is countered by an increase in bare soil (τBS_BLF = +0.25). No significant trend was found in these forests for the percentage of tree cover in bushy sclerophyllous vegetation, demonstrating the greater resistance of Mediterranean maquis to climatic stresses, which is a clear step toward replacing forests with shrubland [15,104]. Areas with positive trend in the percentage of broad-leaved tree canopy coverage (F4 and F3) are mainly located in the inland mountain areas of the island, where water availability is not a problem, and the broad-leaved vegetation might benefit from the warmer weather.

The decreasing trend of broad-leaved tree canopy coverage (βTC_BLF) is positively related to MAP_1990–2019 (ρ = 0.71, p-value 0.004) and negatively related to the Theil–Sen slope of annual VPD (ρ = −0.61, p-value 0.02). On the contrary, the Theil–Sen slope of non-vegetated areas within broad-leaved plots (βBS_BLF) shows a moderate positive correlation with the Theil–Sen slope of annual VPD (ρ = −0.62, p-value 0.01) and a negative one with MAP_1990–2019 (ρ = −0.68, p-value 0.006). This means that, where we have a predominance of broad-leaved vegetation, as the average annual precipitation decreases or the annual VPD increases, the percentage of bare soil increases, while the tree cover decreases. Low or no correlation was found between bushy sclerophyllous vegetation, and the variables analyzed demonstrated low sensitivity at climate drivers. This is likely due to the fact that, drought- and heat-resistant species such as P. latifolia and Pistacia lentiscus, which are common in bushy sclerophyllous vegetation, tolerate the increase of temperature and dry conditions better than more mesic and less cold-sensitive species, such as Q. ilex [105,106].

The use of the proposed simplified βTC_BLF (β_VPDy) and βTC_BLF (MAP_1990–2019) relationships (Figure 9) estimated from past data allowed us to coarsely predict the percentage of tree cover in the future climate scenarios of the GCMs and account for variations in average annual VPD and MAP_1990–2019. Six future air-temperature and precipitation scenarios were derived from the HadGEM2-AO [79] predictions of the CMIP5, identified by Montaldo and Oren (2018) as the most suitable for Sardinia, and a further six scenarios from the following generation of HadGEM configuration of CMIP6, the HadGEM3-GC31-LL predictions. We considered three time periods of increasing distance in the future and two different atmospheric CO₂ concentrations (RCP 45 and RCP 85) for CMIP5, and two Shared Socioeconomic Pathways (SSP 245 and SSP 585) for CMIP6. All scenarios agree on the negative effect of the annual average VPD increase, and F2, F3, F4, F11, and F12 are the forests most affected. As expected, the higher the average annual VPD increase, the higher the tree-canopy-coverage decrease, with SC12 being the worst scenario. Future changes in precipitation were contrasting for CMIP5 scenarios; for example, SC6, in line with past precipitation trends, predicts a further drop in mean annual precipitation (−7.6%), which, based on the proposed model, would correspond to an average 2.1 times greater reduction in the tree cover (−16.09.5%). On the contrary, SC2 predicts a turnaround with an abrupt increase of precipitation (+21.5%) in the upcoming years, with a reduction in the number of forests in water-limited areas and an increase in the percentage of tree cover in almost all forests. Future changes in precipitation for CMIP6 scenarios agree on a precipitation reduction in the range of −3.4% (SC7) to −14.29% (S12). However, although the reduction in precipitation predicted in SC12 is almost double that predicted in SC6, the consequent average reduction in TC is comparable and stands at about −16%. Using the ‘Space-for-time’ substitution approach, a widely used technique to infer future trajectories of ecological systems from contemporary spatial patterns [107], we can predict a change of forest patterns in Sardinia, with more resistant tree species, such as those that compose the bushy sclerophyllous vegetation (P. latifolia, Pistacia lentiscus, Arbutus unedo, Mirtus, etc.). This already applied, for example, in a holm oak forest located in the Montserrat Mountain (Spain) [101], where, after drought-induced events of canopy die-off, there was a change in the composition and structure of the community, involving the increase in species, mostly small shrubs and herbaceous plants, as well as the arrival of non-dominant shrubs and short trees to the canopy. This approach, which follows the principles of adaptive forest management [84], can represent a valid alternative to the conservative management, based on the close-to-natural principles, that is currently required by forestry law in protected areas.

5. Conclusions

This study added evidence about the power of using satellite data for understanding the long-term tree-cover trend in forests and its relationships with climate drivers. We proposed an economic, user-friendly, and expeditious approach, which can be used by the forest management authorities in both developed and developing countries. The method allows users to map the current status of vegetation and historical trends and patterns and to identify future trajectories based on climate drivers. The results show that vegetation is responding to a warming and drier climate in different ways, as a function of species, mean annual precipitation (MAP), and vapor-pressure deficit (VPD) trends. Within broad-leaved vegetation (BLF) with a MAP lower than 700 mm, we are generally witnessing the decline of tall stands in favor of the increase of short stands, while under the same conditions, bushy sclerophyllous vegetation (BSV) showed an opposite behavior. On the contrary, within broad-leaved vegetation with a MAP higher than 700 mm, stands still have the potential to be stable even under the observed positive trends in temperature and VPD. Moreover, regardless of MAP, broad-leaved vegetation showed a drastic negative trend in TC_max (τ = −0.61 and p-value < 0.05), demonstrating that it is not able to recover its original status after the climate disturbances of the last decades. We thus believe that Sardinian broad-leaved forests, as well as others Western Mediterranean oak forests, are declining [84,85,101,108]. In general, the most negative trends in tree-canopy coverage are observed in the broad-leaved vegetation of forests located in the southern-west part of the island, where there are important historical forests composed mainly of tall trees of Quercus ilex. This finding is in line with that of other authors who demonstrate that tall trees, and in particular holm oak, are at high risk of mortality during this changing climate [99,102,103,105]. The medium- and long-term forecasts, based on future climate scenarios and on the simplified models proposed, show that the situation is destined to worsen and deserves attention. It would be useful to demonstrate the compound effects of climate variables and the uncertainties of complex climate effects to vegetation; however, such an analysis is beyond the scope of this paper, which aims only to analyze historical precipitation data in Sardinia, show the actual state of the vegetation and the connection with historical climate data, and to provide simple models that may be easily used by managers and end-users of the forests. The observed overall tree-cover trends may help to identify necessary forest-management activities within the island, whether they are close-to-natural or adaptive. If management of the second type is chosen, the positive trends in tree-canopy coverage observed in most of the bushy sclerophyllous vegetation of the forests analyzed, rather than the absence of correlation with the climate drivers investigated, would lead us to encourage the use of species falling within this category of vegetation. Moreover, it is well-known that forest-cover dynamics have a substantial impact on local water resources [58], and given that Sardinia is experiencing a drastic reduction in runoff [23], future studies should consider the compound effects of those phenomena, even under climate-change scenarios.