Climate Change Impacts on the Water Resources and Vegetation Dynamics of a Forested Sardinian Basin through a Distributed Ecohydrological Model

Abstract

Climate change is impacting Mediterranean basins, bringing warmer climate conditions. The Marganai forest is a natural forest protected under the European Site of Community Importance (Natura 2000), located in Sardinia, an island in the western Mediterranean basin, which is part of the Fluminimaggiore basin. Recent droughts have strained the forest′s resilience. A long-term hydrological database collected from 1922 to 2021 shows that the Sardinian forested basin has been affected by climate change since the middle of the last century, associated with a decrease in winter precipitation and annual runoff, reduced by half in the last century, and an increase of ~1 °C in the mean annual air temperature. A simplified model that couples a hydrological model and a vegetation dynamics model for long-term ecohydrological predictions in water-limited basins is proposed. The model well predicted almost one century of runoff observations. Trees have suffered from the recent warmer climate conditions, with a tree leaf area index (LAI) decreasing systematically due to the air temperature and a vapor pressure deficit (VPD) rise at a rate of 0.1 hPa per decade. Future climate scenarios of the HadGEM2-AO climate model are predicting even warmer conditions in the Sardinian forested basin, with less annual precipitation and higher air temperatures and VPD. Using these climate scenarios, we predicted a further decrease in runoff and tree transpiration and LAI in the basin, with a reduction of tree LAI by half in the next century. Although the annual runoff decreases drastically in the worst scenarios (up to 26%), runoff extremes will increase in severity, outlining future scenarios that are drier and warmer but, at the same time, with an increased flood frequency. The future climate conditions undermine the forest’s sustainability and need to be properly considered in water resources and forest management plans.

1. Introduction

Over the past century, climate change has been altering precipitation and air temperature regimes across the world [1,2,3,4,5,6,7]. Precipitation and runoff are decreasing in Mediterranean regions, with impacts on water resources and ecosystem sustainability [8,9,10]. In Sardinia, an island in the Mediterranean Sea, Montaldo and Sarigu (2017) [10] estimated a persistent negative trend of winter precipitation (Mann−Kendall τ up to −0.40) in the 1922−2010 period, which produced an even more dramatic negative trend of runoff (τ up to −0.45), with runoff more than halved in recent decades, and a decrease in the mean annual basin runoff coefficient (Mann-Kendall τ value of −0.32 [11]. Runoff is a key term of basin water resources and a good indicator of water availability [12,13]. The impact on water resources has been dramatic in Sardinia, with frequent water supply restrictions, including for domestic consumption [14]. A further negative precipitation trend has been predicted for the near future [15,16], affecting the central Mediterranean basin and speeding up the desertification process [16,17,18,19,20,21].

At the same time, the air temperature is increasing around the world, and the positive trend of air temperature impacts the vapor pressure deficit (VPD) trend [22,23] and consequently the vegetation dynamics, which are strictly related to temperature and VPD (e.g., [24,25]).

Forests are important for carbon sequestration, erosion control [26], and biodiversity conservation [27] and need to be preserved [28,29]. Forests make up 30.74% of the planet’s surface [30], and in Mediterranean countries, the forest area is 88 million ha [31]. The declining trend of precipitation related to climate change can affect forest sustainability [32,33] due to the direct relation between precipitation and tree cover development in water-limited ecosystems. Sankaran et al. (2005) [34] and Axelsson and Hanan (2017) [35] estimated a mean annual precipitation (MAP) threshold of around 650−700 mm, below which rain restricted tree cover growth in African savannas. Understanding of forest response to decreased precipitation is still undetermined because tree strategies for reducing vulnerability to drought can differ (e.g., resource partitioning, tree species selection) related to tree species diversity and richness [36].

In Mediterranean forests, recent investigations on holm oaks showed that drought reduced their growth and increased tree mortality [37,38]. Typically, to survive dry conditions, trees optimize the use of water by extending their roots downward, where soils are wet for longer [39,40,41,42]. In Mediterranean basins, especially under water-restricted conditions, upstream hillsides are often characterized by trees growing on shallow soils above fractured bedrocks or cemented horizons, where tree roots can still enter through cracks, taking up water from the rock storage (up to 70–90% of transpiration) through the hydraulic lift process [40,43,44,45,46,47,48,49]. Drier conditions can affect tree development, making past survival and water-use optimization strategies insufficient, and these need to be investigated further. Furthermore, an increase in VPD due to the recognized increase in air temperature may affect vegetation growth negatively [23,50,51], and subsequently impact land surface fluxes and water resources [24,52,53]. Forest sustainability under climate change warming is still strongly questioned [54,55] and needs to be investigated using accurate modeling approaches.

Currently, long-term databases of land-observed hydrologic data (e.g., rain, runoff and air temperature) have been available for most of the last century, and the effects of climate change on water resources and vegetation dynamics can be evaluated [3,32,37]. At the same time, the availability of future climate scenarios may allow for the prediction of the states and the preservation of water and vegetation resources in the future [56,57,58]. In general, the effects of climate and land use changes may be overlapped and need to be separated [59]. A first step is to investigate climate change effects only, looking at areas with low human impact and few anthropogenic land cover changes in the past. In this sense, having relatively low urbanization and human activity (with ~50% of the total area covered by forest; [60]), Sardinia is an excellent reference for ecohydrological studies on past and future climate change effects [8,11,18,47,61,62] representing typical conditions of the western Mediterranean Sea basin, and endowed with an extensive and long-term hydrological database dating back almost a century [10,47].

Remote sensing observations of the vegetation index, like the normalized difference vegetation index (NDVI), provide spatial monitoring of forest dynamics at high time resolution [63,64]. Long NDVI databases are starting to become available, although they are not as extensive as the hydrological databases. Time series of the Advanced Very High Resolution Radiometer (AVHRR) optical sensor are, for instance, greatly extended (from 1984), but the coarse spatial resolution (1.1 km) is inappropriate for Mediterranean hydrologic basins characterized by strong spatial heterogeneity. Instead, the time series of NDVI from the MODIS sensor on the Aqua and Terra satellite platforms are robust and still long enough (from 2000), with a smaller spatial resolution (~250 m) [64,65], to offer appropriate support for ecohydrological and climate change studies.

Ecohydrological models are a key tool for predicting climate change impacts on water and vegetation resources, through the coupling of hydrologic and vegetation dynamics models [52,66,67,68,69,70,71]. With the growth of computational power, ecohydrological models have increased in complexity through the inclusion of more processes, controlling variables and parameterization [52,72,73]. These parameter additions are justified on the basis that the inclusion of additional controlling mechanisms should both improve the predictive skill and make the parameter values easier to estimate based on physiological characteristics or measurements, but the uncertainty of the parameterization has increased [52,73,74,75]. Hence, there is a need for models with parsimonious parametrization that are still able to accurately predict the main drivers of the processes, with a low computational load, becoming key requirements when long databases are simulated. We propose an ecohydrological model for both hydrological and vegetation dynamics modeling under water-limited conditions, which respond to the need for: (1) long-term (e.g., hundreds of years) hydrologic balance predictions, (2) the inclusion of vegetation dynamics, (3) the inclusion of the hydraulic redistribution by tree roots from deep to surface soil layers, an important contribution in water-limited forested ecosystems, (4) a spatially distributed approach that is able to capture the spatial variability of ecohydrological processes, physiographic properties of the basin and meteorological forcing, and (5) parsimonious model parameterization and a low computational load. With these objectives, daily or monthly model time scales can be adequate for water balance and vegetation dynamics predictions [76,77,78]. We started from the basin-scale distributed hydrologic model of Montaldo et al. (2007) [79] and the plot-scale ecohydrological model of Montaldo et al. (2008) [52] and Montaldo et al. (2013) [61], which couples a land surface model and a vegetation dynamic model. At a daily time scale, the Hortonian overland flow [80,81,82], typical of Mediterranean basins with steep hillslopes in semi-arid climate conditions, can be predicted using the Soil Conservation Service (SCS) method [81,83,84], which is widely used for surface runoff estimation (e.g., [85,86]). For the recharge of the shallow soil by tree roots through the hydraulic redistribution, recently, Montaldo et al. (2021) [47] successfully developed a simplified approach, which we adapted in our proposed model. We used the vegetation dynamic model of Montaldo et al. (2008) [68], running at a daily time scale and already tested in a Sardinian ecosystem for LAI predictions, which is parsimonious in model parameterization [52].

In Sardinia, the Marganai forest (mainly Quercus ilex trees) is a European Site of Community Importance (Natura 2000), part of the long-term ecological research network [87]), and is mainly in the Rio Fluminimaggiore basin, located on the southwest side of the island. In this part of Sardinia, a significant past decrease in rainfall was estimated by Montaldo and Sarigu (2017) [10], due to its strong correlation with the North Atlantic Oscillation (NAO) index. The Marganai forest experienced a recent drought in summer 2017 that affected the trees (Figure 1b), sounding a warning signal on the sustainability of this protected forest. Montaldo and Oren (2018) [11] showed a mean decreasing trend of evapotranspiration (ET) in Sardinian basins (~−5% for the period 1951−1975 compared to 1976−2000), and that ET was mostly affected by the decreasing of precipitation. The proposed ecohydrological model is used for predicting runoff, soil water balance and forest tree behavior during the long 1924−2021 period at the Fluminimaggiore basin, for which hydrological observations are available. Remote optical sensor observations of MODIS were used for estimating the NDVI of the basin and testing model predictions of vegetation dynamics. The ecohydrological model is also used for future scenario predictions. We used the future climate scenarios predicted by global circulation models (GCM) in the Fifth Assessment report of the Intergovernmental Panel on Climate Change (IPCC) [88], and selected the most reliable models after testing the past GCM predictions against observed data [11]. In summary, the objectives of this study are: (1) to develop a parsimonious spatially distributed model for long-term ecohydrological studies in water-limited basins; (2) to evaluate the past (almost one century) climate change effects on the soil water balance terms (i.e., runoff, and evapotranspiration), and tree LAI of the forested Fluminimaggiore basin, through the use of the proposed ecohydrological model, and (3) to predict future ecohydrological scenarios of the forested Sardinian basin, evaluating the impact on water resources and the hydrological state of the protected Marganai forest.

2. Methods

2.1. The Rio Fluminimaggiore Basin and the Marganai Forest

The case study is in the Rio Fluminimaggiore basin (area 83 km²), which is located in the southwest of Sardinia (Italy) (Figure 2), and includes the Marganai Forest, a Long–Term Ecosystem Research (LTER) Italian site and a European Site of Community Importance (Natura 2000).

The basin is predominantly mountainous (89% of the basin with height > 200 m a.s.l.), with a mean slope of ~24°. The mountains are granitic and limestone rocks covered mainly by thin soils (mostly silty-clay) ranging from 10 to 50 cm deep. The basin is mainly natural (~97% of the basin), and the predominant land cover is forest (~55% of the basin), with half being dense Quercus ilex forest, and the rest typical Mediterranean maquis (Arbutus unedo L., Phillyrea latifolia L., Euphorbia dendronides L., Cistus monspeliensis L., Cistus salvifolius L., Rosmarinus officinalis L.). Typically Mediterranean basins in water-limited conditions, tree roots expand below the thin soil layer into the fractured rocky layers, looking for rock moisture [47].

The climate is typical marine Mediterranean climate type, characterized by cool and wet winters and warm and dry summers. The mean historical annual precipitation is 782 mm (1922−2021) with mean monthly precipitation ranging from 125 mm (in December) to 3 mm (in July). The mean annual air temperature is 16.36 °C, with mean monthly values ranging from 9.45 °C (in January) to 24.29 °C (in August).

Daily rain data from 8 stations of the Sardinian Regional Hydrographic Service are available in the area of the basin (Figure 2), with data from 1922 to 2021. Daily data of the hydrometric station of the Fluminimaggiore river are available from 1924 to 2021 (a record length of 97 years and 22.68% of missing data). In the following, yearly estimates of hydrologic variables are computed for the hydrologic year, which in Sardinia starts in October (when precipitation usually starts) and ends in September of the next year (when the soil is dry) [11]. Data from the 8 thermometric stations of the Sardinian Regional Hydrographic Service are available, with at least 97 complete years of data (Figure 2), providing daily maximum and minimum air temperatures (T_max and T_min, respectively). We also used wind velocity data from a nearby anemometric station at Capo Frasca, and relative humidity and incoming solar radiation data from the ERA20C (data from 1924 to 1978, [89]) and the ERA5 (data from 1979 to 2021, [90,91] European Centre for Medium-Range Weather Forecasts (ECMWF) atmospheric reanalysis datasets, which we tested successfully using two years of data of incoming solar radiation from two stations (Punta Gennarta and Rio Leni) close to the basin (RMSE = 90.46 W/m²), and 15 years of data of relative humidity at a Sardinian station in Orroli (distance of 65 km; RMSE = 19.02%).

A digital elevation model was available (Figure 2), from which basin topographic properties are derived at 200 m spatial resolution. Soil characteristics (soil depth, texture, and hydrologic properties) and land cover maps are also available at the same spatial resolution. We used the soil depth map made by the Sardinian authority for agricultural scientific research and technological innovation (AGRIS) [92]. It was at 100 m spatial resolution, and we up-scaled at 200 m by averaging the 100 m data. The AGRIS soil depth map was also tested successfully by independent measurements of the forest Sardinian authority (FORESTAS).

NDVI is estimated from red and near-infrared spectral reflectance measurements of remote sensor data. We used MODIS data at 250 m spatial resolution and 16 days of temporal resolution, available from 2000 to 2021 (Figure 1). The LAI is derived from NDVI MODIS data through a non-linear relationship (LAI = 7658 NDVI^33.06 + 2.97 for the forest; [93]. Sentinel 2 images, which were characterized by higher spatial resolution (10 m), were also acquired, but they are available from 2016 only.

2.2. The Distributed Ecohydrological Model

The proposed numerical model is a spatially distributed ecohydrological model (Figure 2 and Figure 3), which couples a hydrologic model and a vegetation dynamic model. The model runs at a daily timescale and 200 m spatial scale. Three principal components can be identified in the model (Figure 3): (1) the soil water balance model (SWBM) that computes soil moisture, water balance, surface runoff and subsurface water flux, (2) the vegetation dynamic model (VDM) that computes biomass budgets and LAI, and (3) the simplified overland flow and base flow computations.

In the first component, the soil water budget is computed at each cell of the basin. From this model component, the surface runoff and drainage are predicted, and these are used by the third component of the model for computing runoff at the basin scale. The second component, the vegetation dynamic model, provides the LAI at each cell, which is used by the SWBM for computing the evapotranspiration (ET) and rain interception (Figure 3).

Meteorological model inputs are: precipitation (P), air temperature (T_air), wind speed (WS), relative humidity (RH), incoming short wave solar radiation (Rsw) and photosynthetically active radiation (PAR). We used the Thiessen method for spatial interpolation of the meteorological model inputs. The parameters of the model are in

Table 1. Model parameters of the vegetation components.

Parameter	Description	Grass	Forest
rs,min (s m−1)	Minimum stomatal resistance	150	290
θwp (-)	Wilting point	0.08	0.05
θlim (-)	Limiting soil moisture for vegetation	0.20	0.18
Tmin (°K)	Minimum temperature	279.15	268.15
Topt (°K)	Optimal temperature	293.15	288.15
Tmax (°K)	Maximum temperature	299.15	304.15
ca (m2 gDM−1)	Specific leaf areas of the green biomass	0.01	0.0065
cd (m2 gDM−1)	Specific leaf areas of the dead biomass	0.01	0.0062
ke (-)	PAR extinction coefficient	0.5	0.5
ξa (-)	Parameter controlling allocation to leaves	0.6	0.55
ξs (-)	Parameter controlling allocation to stem	0.1	0.1
ξr (-)	Parameter controlling allocation to roots	0.4	0.35
Ω (-)	Allocation parameter	0.8	0.1
ma (d−1)	Maintenance respiration coefficients for aboveground biomass	0.032	0.0001
ga (-)	Growth respiration coefficients for aboveground biomass	0.32	0.85
mr (d−1)	Maintenance respiration coefficients for root biomass	0.007	0.0003
gr (-)	Growth respiration coefficients for root biomass	0.1	0.1
Q10 (-)	Temperature coefficient in the respiration process	2.5	3
da (d−1)	Death rate of aboveground biomass	0.023	0.0019
dr (d−1)	Death rate of root biomass	0.005	0.0001
ka (d−1)	Rate of standing biomass pushed down	0.23	0.35
QN (-)	Soil respiration coefficient related to temperature	1.2
R10 (mmol CO2/m2s)	Reference respiration at 10 °C	2.54
zom,v (m)	Vegetation momentum roughness length	0.05	0.5
zov,v (m)	Vegetation water vapor roughness length	zom/7.4	zom/2.5
zom,bs (m)	Bare soil momentum roughness length	0.015
zov,bs (m)	Bare soil water vapor roughness length	zom/10

.

2.2.1. The Soil Water Balance Model

The SWBM predicts land surface fluxes and soil moisture evolutions for each grid cell. It is derived from the models of Montaldo et al. (2013) [78] and Montaldo et al. (2021) [47], and includes three land surface components for each cell: bare soil, grass and trees. Because the basin is mainly characterized by thin soils above fractured rocky layers, the model includes two layers vertically, a surficial soil of variable depth (d_s) according to the soil depth of the cell, and an underlying active fractured rock depth (d_r) for modeling root zone soil moisture. The model simulates the water balance of the two layers for each cell:

$$d_s\frac{\partial {\theta }_s}{\partial t}=P-E_w-Q_{sup}-f_{bs}{E}_{bs}-f_g{E}_g-{\xi }_t f_t{E}_{t,s}+f_d-D_r$$

(1)

$$d_{r }\frac{\partial {\theta }_r}{\partial t}=D_r-f_d-\left(1-{\xi }_t\right)f_t{E}_{t,r}-L_e$$

(2)

where θ_s is the moisture of the upper soil layer, θ_r is the moisture of the underlying rocky layer, P is the precipitation, Q_sup is the surface runoff, E_w is the wet evaporation, which is equal to the rainfall interception, characterized by a storage capacity of 0.2 LAI [94], E_bs is the bare soil evaporation, E_g is the grass transpiration, E_t,s is the tree transpiration from the surface soil layer, E_t,r is the tree transpiration from the fractured rock layer (with tree transpiration, E_t = E_t,s + (1 − ${\xi }_t$) E_t,r), ${\xi }_t$ is the percentage of tree root water uptake from the surface soil layer, f_d is the daily hydraulic redistribution flux between the surface soil and the fractured rock through the tree roots, D_r is the vertical drainage, L_e is the leakage, f_bs, f_g and f_t are the fractions of bare soil, grass cover and tree cover in each cell, respectively, with f_bs + f_g + f_t = 1. Q_sup is estimated using the Soil Conservation Service Curve Number (SCS-CN) method [81,83,84] from P, and computing the antecedent moisture conditions from the rain in the previous five days (Appendix A). We used the original seasonal rainfall limits of the SCS-CN method for distinguishing the antecedent moisture conditions (Appendix A) and calibrated the CN map (for properly predicting runoff). D_r and L_e are estimated using the unit gradient assumption of gravity drainage [95] and related to soil moisture following Clapp and Hornberger (1978) [96]:

$$D_r=k_{sat,s}{\left(\frac{{\theta }_s}{{\theta }_{sat,s}}\right)}^{2b_s+3}$$

(3)

$$L_e=k_{sat,r}{\left(\frac{{\theta }_r}{{\theta }_{sat,r}}\right)}^{2b_r+3}$$

(4)

where k_sat,s and k_sat,r are the saturated hydraulic conductivity of the surficial soil layer and the deeper layer, respectively, b_s and b_r are the slopes of the soil retention curve of the two layers, and θ_sat,s and θ_sat,r are the saturated soil moisture of the two layers.

The f_d contribution is estimated as a function of the soil moisture gradient between the surface and underlying rocky sublayer through [47,97]:

$$f_d=a_f{\left({\theta }_r-{\theta }_s\right)}^{{\mathrm{b}}_f},$$

(5)

where a_f and b_f are tree root parameters. Evapotranspiration components (E_g, E_t,s, and E_t,r) are estimated by the Penman−Monteith equation ([98] Equation 10.34), with the canopy resistance (r_c) given by [99]:

$$r_c=\frac{r_{s,min}}{\mathrm{LAI}}{\left[f_1\left(\mathsf{\theta }\right)f_2\left(T_a\right)f_3\left(VPD\right)\right]}^{-1}$$

(6)

where f₁, f₂ and f₃ are stress functions of soil/rock moisture, air temperature and VPD, calculated as [11,68]:

$$f_1\left(\theta \right)=\{\begin{array}{cc}0& for \theta \le {\theta }_{wp}\\ \frac{\theta -{\theta }_{wp}}{{\theta }_{\lim }-{\theta }_{wp}}& for {\theta }_{wp}<\theta \le {\theta }_{\lim }\\ 1& for \theta >{\theta }_{\lim }\end{array}$$

(7)

$$f_2\left(T_a\right)=\{\begin{array}{cc}0& for T_a\le T_{min} | T_a>T_{max}\\ 1-\frac{T_{opt}-T_a}{T_{opt}-T_{min}}& for T_{min}<T_a<T_{opt}\\ 1& for T_{opt}\le T_a\le T_{max}\end{array}$$

(8)

$$f_3=1-0.6 \log \left(VPD\right)$$

(9)

where θ_wp is the wilting point, θ_lim is the limiting soil moisture for vegetation, and T_min, T_opt and T_max are the minimum, optimal and maximum temperature for vegetation. Note that θ soil moisture contributors are θ_s for E_g and E_t,s and θ_r for E_t,r.

Bare soil evaporation is estimated as a function of soil moisture through E_bs = α PE [100], and PE is the potential evaporation based on (Penman−Monteith [98], Equation 10.15). Total evapotranspiration is given by:

$$ET=f_{bs}{E}_{bs}+f_g{E}_g+{\xi }_t f_t{E}_{t,s}+\left(1-{\xi }_t\right)f_t{E}_{t,r}+E_w$$

(10)

Note that the soil water balance model for the root zone can include a soil sublayer instead of the root-accessible fractured rock sublayer based on the same approach, but it would require different parameter values.

2.2.2. The Vegetation Dynamic Model

The VDM computes the change in biomass over time from the difference between the rates of biomass production (photosynthesis) and loss, such as what occurs through respiration and senescence (e.g., [101,102]). The VDM distinguishes trees and grass components and is derived from the models of Nouvellon et al. (2000) [103], Montaldo et al. (2005) [52] and Montaldo et al. (2008) [68]. In the VDM of trees, four separate biomass states (compartments) are tracked: green leaves (B_g), stem (B_s), living root (B_r), and standing dead (B_d). The biomass (g DM m⁻²) components are simulated by ordinary differential equations integrated numerically at a daily time step [52,101,103,104]:

$$\frac{d{B}_a}{dt}=a_{as}{P}_{g,s} {\xi }_t+a_{ar}{P}_{g,r}\left(1-{\xi }_t\right)-R_a-S_a$$

(11)

$$\frac{d{B}_s}{dt}=a_{ss}{P}_{g,s} {\xi }_t+a_{sr}{P}_{g,r}\left(1-{\xi }_t\right)-R_s-S_s$$

(12)

$$\frac{d{B}_r}{dt}=a_{rs}{P}_{g,s} {\xi }_t+a_{rr}{P}_{g,r}\left(1-{\xi }_t\right)-R_r-S_r$$

(13)

$$\frac{d{B}_d}{dt}=S_a-L_a$$

(14)

where P_g is the gross photosynthesis, a_as, a_ar, a_ss, a_sr, a_rs, and a_rr are allocation (partitioning) coefficients to leaves, stem and root compartments, R_a, R_s and R_r are the respiration rates from leaves, stem and root biomass, respectively, S_a, S_s and S_r are the senescence rates of leaves, stem and root biomass, respectively, and L_a is the litter fall. The equations of the terms in (11)−(14) are in

Table 2. Equations of the vegetation dynamic model components. Parameters are defined in Table 1.

Ecophysiological Term	Equation
Photosynthesis	$P_g={\epsilon }_P\left(PAR\right)f_{PAR}PAR\frac{1.37r_a+1.6r_{c,min}}{1.37r_a+1.6r_c}$ ${\epsilon }_P\left(PAR\right)=a_{0}PAR/\left(1+a_{1}PAR\right)$ $f_{PAR}=1-e^{-k_e\mathrm{LAI}}$
Allocation	For woody vegetation
	$a_{aj}=\frac{{\xi }_a}{1+\mathrm{\Omega }[2-\lambda -f_1\left({\theta }_j\right)]}$
	$a_{sj}=\frac{{\xi }_s+\mathrm{\Omega }\lambda }{1+\mathrm{\Omega }[2-\lambda -f_1\left({\theta }_j\right)]}$
	$a_{rj}=\frac{{\xi }_r+\mathrm{\Omega }(1-f_1\left({\theta }_j\right)}{1+\mathrm{\Omega }[2-\lambda -f_1\left({\theta }_j\right)]}$ $\lambda =e^{-k_{e}LAI}$ ${\xi }_a+{\xi }_s+{\xi }_r=1$
	For grass
	$a_a=\frac{{\xi }_a+\mathrm{\Omega }\lambda }{1+\mathrm{\Omega }[1+\lambda -f_1\left({\theta }_s\right)]}$
	$a_r=\frac{{\xi }_s+\mathrm{\Omega }(1-f_1\left({\theta }_s\right)}{1+\mathrm{\Omega }[1+\lambda -f_1\left({\theta }_s\right)]}$
	${\xi }_a+{\xi }_r=1$
Respiration	$R_a=m_a{f}_4(T)B_a+g_a\left[a_{as}{P}_{g,s}{\xi }_t+a_{ar}{P}_{g,r}\left(1-{\xi }_t\right)\right]$
	$R_s=m_s{f}_4(T)B_s+g_s[a_{ss}{P}_{g,s}{\xi }_t+a_{sr}{P}_{g,r}(1-{\xi }_t)]$
	$R_r=m_r{f}_4(T)B_r+g_r[a_{rs}{P}_{g,r}{\xi }_t+a_{rr}{P}_{g,r}(1-{\xi }_t)]$ $f_4(T)=Q_{10}^{Ta/10}$
Senescence	$S_a={\delta }_a{B}_a$
	$S_s={\delta }_S{B}_S$
	$S_r={\delta }_r{B}_r$
Litterfall	$L_a=k_a{B}_d$

. The key term of the VDM, P_g, is computed using the approach of Montaldo et al. (2005) [52], which started from a simplified form of Fick’s law applied to gas exchange in plants [102,105] and the Nouvellon et al. (2000) model [103], deriving a simplified expression that estimates P_g mainly by the photosynthetically active radiation PAR, soil moisture, and other routinely monitored variables (wind velocity, and air humidity and temperature). Note that in Equations (11)−(13), the θ soil moisture contributors are θ_s for P_g,s and θ_r for P_g,r.

The VDM of grass distinguishes only three biomass compartments (green leaves, roots and standing dead), and the biomass components are simulated using Equations (11), (13) and (14), respectively.

Leaf area index values are estimated from the biomass through linear relationships [52,103,105,106]:

$$\mathrm{LAI}=c_{a }{B}_a$$

(15)

where c_a is the specific leaf area of the green biomass. The other VDM equations are reported in Table 2. The SWBM is coupled with the VDM, which provides the LAI values of trees and grass daily, which are then used by the SWBM for computing the evapotranspiration partitioning, and the soil water content. The SWBM provides soil moisture and aerodynamic resistances to the VDM.

2.2.3. The Surface Runoff and Base Flow

In each cell of the basin, the Q_sup surface runoff is computed daily by the SWBM, and the total surface runoff at the basin scale is computed by summing the surface runoff contribution of the cells of the basin daily. This model simplification is reasonable for small basins such as the Fluminimaggiore basin, which are characterized by a small time of concentration (<4 h). Following Montaldo et al. (2007) [79], the Q_b subsurface flow at the basin outlet is computed at the j + 1 time step through a lumped conceptual approach, namely a linear reservoir method [107]:

$$Q_b^{j+1}=Q_b^j e^{-\Delta t/\kappa }+L_e^{j+1}\left(1-e^{-\Delta t/\kappa }\right)$$

(16)

where Δt is the time step (=daily), and k is a linear reservoir parameter.

At each daily time step, the total runoff at the basin outlet (Q) is given by the sum of total surface runoff and Q_b.

2.3. Other Field Experimental Sites for Testing Model Components

In the Fluminimaggiore basin, measurements of some key outputs of the SWBM, such as evapotranspiration of the main tree species, Quercus ilex, and the soil water balance for shallow soils with trees, a typical feature of the Marganai forest, are not available. For testing these model components, we used the data of two other case studies, the Castelporziano and Orroli sites, which are still in Mediterranean climate conditions, and for which the required measurements are available. We chose these two sites because the trees in Castelporziano are of the same species as the Marganai forest, while the Orroli site is characterized by a long dataset of soil moisture and micrometeorological data [108,109], a similar morphological soil type with a key role of the hydraulic redistribution and the water uptake by the tree roots from the underlying fractured rocks, and is located in Sardinia, 65 km from the Fluminimaggiore basin.

The Castelporziano Estate is located in Lazio, central Italy (41.70427°, 12.3572°, 19 m elevation), and is a Long-Term Ecological Research site with data available on the (FLUXNET2015 Dataset [110], site ID: IT-Cp2). The site is located in a coastal area and covers about 6000 ha [111]. The climate is typical Mediterranean semi-arid with a mean annual precipitation of 805 mm, similar to the Fluminimaggiore basin, and a mean annual temperature of 15.2 °C [111]. The soil is a sandy-loam with a mean depth of 45 cm. Vegetation is mainly composed of Quercus ilex with a percentage cover ~60% [111,112]. This site is equipped with an eddy covariance tower for measuring evapotranspiration, CO₂, and energy fluxes (e.g., [113]).

The Orroli site is located in east-central Sardinia (39°41′12.57″ N, 9°16′30.34″, 500 m elevation), with thin soils (mean soil depth of 15 cm). The climate at this site is maritime Mediterranean, with a mean annual precipitation of 643 mm and mean annual air temperature of 14.6 °C. The vegetation is a patchy mixture of wild olive tree (Olea Sylvestris) clumps with canopy cover of ~33% and herbaceous seasonal species. The site has been monitored since 2003 by a micrometeorological tower for the estimation of sensible heat flux, evapotranspiration and CO₂ exchanges with the standard eddy covariance method. Seven frequency domain reflectometer probes (FDR, Campbell Scientific Model CS-616) were inserted in the soil close to the tower (3.3−5.5 m away) to estimate moisture at half-hourly intervals in the thin soil layer. FDR calibration and the other instrumentation of the field site are described in Montaldo et al. (2020) [109].

2.4. IPCC Future Projections

For future climate scenarios, we have considered the IPCC scenarios in the Fifth Assessment report, focusing on the Atmosphere–Ocean General Circulation Models (AOGCMs) of Flato et al. (2014) [89]. We selected the climate model for future scenarios after comparing the model outputs of air temperature and precipitation with the observed data in the Fluminimaggiore basin. Following Montaldo and Oren (2018) [11], from the 12 AOGCMs (MPI-ESM-LR, CMCC-CMS, CNRM-CM5, GFDL-ESM2G, GFDL-ESM2M, HadCM3, HadGEM2-CC, HadGEM2-ES, MPI-ESM-MR, NorESM1-M, NorESM1-ME, HadGEM2-AO) of Flato et al. (2014) [88], we selected the HadGEM2-AO, because it predicted better than the others the observed past changes of precipitation and air temperature in the Sardinian basin. We compared the changes of precipitation and air temperature in the period 1976−2000 with 1951−1975. The HadGEM2-AO model was able to capture past climate changes, with differences in predicted changes of winter precipitation of less than 2% of the observed changes (less than 4% for yearly data), and with low differences (~0.1 °C) in yearly changes of air temperature.

Future climate scenarios are generated up to 2100 for the Fluminimaggiore basin from HadGEM2-AO predictions through the multivariate bias correction technique [114], which uses the land-observed meteorological data of the past for estimating the required statistics. Future climate scenarios are generated up to 2100 for the Fluminimaggiore basin from HadGEM2-AO predictions through the multivariate bias correction technique [114], which uses the land-observed meteorological data of the past for estimating the required statistics. The reference period is from 1950 to 2005. The stations with land-observed data are in Figure 2, and basin areal averages are computed using the Thiessen method. The multivariate bias correction technique uses an iterative algorithm that consists in three main steps: (1) constructing a uniformly distributed random orthogonal rotation matrix and applying that to the land-observed meteorological target data and to the AOGCM meteorological source data of the past, (2) correcting the marginal distributions of the rotated sources data by using an empirical quintile map, (3) applying the inverse rotation to the resulting data. These steps are repeated until the multivariate distribution converges to the target distribution. The multivariate bias correction statistically downscaled the AOGCM data at a resolution of 1.25° latitude and 1.87° longitude (~150 km) to the basin spatial scale. The temporal resolution of the HadGEM2-AO predictions is the same as the model time step (daily). Four representative concentration pathways (RCP) of HadGEM2-AO are considered (RCP 26, RCP 45, RCP 6, RCP 85).

2.5. Comparisons and Statistical Data Analysis

Data were analyzed on daily, monthly, seasonal and yearly time scales. Model goodness of fit was evaluated by comparing the modeling results with observations and using the following statistics: mean (μ), standard deviation (SD), coefficient of variation (CV), mean error (me), root mean square error (RMSE), correlation coefficient (ρ), Nash-Sutcliffe efficiency (NSE) [115] and mean ratio (R_μ).

For evaluating trends in the data series, we used the Mann−Kendall (1938) method [116], which allows us to estimate the τ index of the trend, and the Theil−-Sen method (1950,1968) [117,118] for estimating the β trend line slope. The Mann−Kendall and Theil−Sen methods are widely used in ecohydrological studies of climate change [10,119,120].

Finally, we used the Generalized Extreme Value (GEV) distribution for the frequency analysis of the annual maximum runoff and Kolmogorov–Smirnov and Chi-squared tests for evaluating the model goodness. The GEV distribution is widely used in hydrological applications (e.g., [121,122] and is well suited for the Sardinian basins [123].

3. Results

A long drought affected the Marganai forest in summer 2017, due to an unusual dry spring (48.37 mm of rain compared with the average 109.37 mm) that preceded a common dry summer. At the end of August 2017, the south-facing Quercus ilex trees turned red due to the critical state (Figure 1b), while usually, they were still green at the end of the summer (Figure 1a). The state of the trees was clearly identified by Sentinel 2 images (Figure 1c,d) due to their high spatial resolution (10 m). However, the use of the coarse MODIS images still allowed us to capture the lower NDVI of the tree cover at the end of August 2017 (mean NDVI ~0.7), compared to the end of August 2016 (mean NDVI ~0.8; Figure 1c,d). Looking at the long MODIS NDVI time series of the Marganai forest, in summer 2017, the NDVI actually reached the minimum value (~0.7) in the last 13 years; however, in the first 8 years of the investigated 2000−2021 period, the NDVI decreased to 0.7 another five times (Figure 1g); highlighting it was more common for trees to reach critical conditions in the first decade of the twenty-first century. However, 2017 was characterized by a longer and unusual drought (from spring) that caused a longer series of very low NDVI values (from April; Figure 1g).

The Fluminimaggiore basin is characterized by low urbanization conditions and low human activities, which has not changed from the past, as can be seen from the comparison of historical (1954) and recent (2016) aerial photography (Figure S1). The basin is still mostly natural (~97% of the basin), forested (55%), and the increase in urban areas is negligible (from 0.45% to 0.79% of the basin area), while the mine area (1.4%) is still present, although no longer active.

The analysis of the historical precipitation and temperature data at the basin scale pointed out that the yearly precipitation trend is negative but not significant (Mann–Kendall τ of −0.065 with p = 0.345, Theil Sen β of −0.596 mm/year; Figure 4a), while the yearly air temperature is increasing significantly (Mann–Kendall τ of 0.300 with p < 0.001, Theil–Sen β of 0.010 °C/year; Figure 4b). Due to the precipitation regime of the basin (Figure 5a), precipitation is mainly concentrated in the autumn and winter months, with mean monthly precipitation ranging from 72 mm in March to 125 mm in December (Figure 5a), and the decrease in precipitation was significant in winter (Mann–Kendall τ of −0.15 with p < 0.05, and Theil–Sen β of −0.26 mm/year; Figure 5b,c). Conversely, the air temperature increase was higher in the spring and summer months (Mann Kendall τ up to 0.26 in summer; Figure 5b,c), when the air temperature was already much higher (up to 25 °C on average in summer; Figure 5a). Due to the positive air temperature trend, the mean annual VPD increased significantly (Mann–Kendall τ of 0.13 with p = 0.06, Theil–Sen β of 0.01 hPa/year).

The negative trend of observed annual runoff was significant and higher than the negative trend of annual precipitation, reaching a Mann–Kendall τ of −0.26 (p = 0.001) and a Theil−Sen β of −2.206 mm/year (Figure 4c). In the 1990s, annual runoff was very low, with values lower than 100 mm/year for three years (Figure 4c).

3.1. The Ecohydrological Model Testing and Its Predictions for the Past Period

Firstly, the model components were tested locally at the two experimental sites without considering surface and subsurface flows. The ET model components were tested using the data of the Castelporziano site. The model parameter values are in Table 1 and

Table 3. Soil-related parameters of the model.

Parameter	Description	Mean	Range
CN	Curve Number of the Soil Conservation Service method	89.7	80–99
θsat,s (-)	Saturated soil moisture	0.44	0.41−0.46
bs (-)	Slope of the retention curve in the surface soil	10.28	7.75−11.40
ksat,s (m/s)	Saturated hydraulic conductivity	2.82 × 10−7	10−8−10−6
ds (m)	Soil layer depth	0.36	0.20−0.85
qsat,r (-)	Saturated rock moisture	0.48	0.45−0.50
br (-)	Slope of the retention curve in the fracture rock layer	7
ksat,r (m/s)	Saturated hydraulic conductivity in the fracture rock layer	1.41 × 10−7	5 × 10−9−5 × 10−7
dr (m)	Fractured rock depth	1.5

. The model component for the evapotranspiration estimate of Quercus ilex was successfully tested using the 2-year Castelporziano ET observations (RMSE = 0.039 mm/d, R² = 0.73 and p < 0.05; Figure 6a,b). The model well predicted the seasonal variability of ET, which reached 4 mm/d in summer (Figure 6a). Soil moisture and f_d (a_f = 169.83, b_f = 8.61 in Equation (6); ξ_t = 0.3 if θ_s > = 0.2, otherwise ξ_t = 0.1) are successfully tested using the long data base (15 years) of the Orroli site (RMSE = 0.067 R²= 0.78 and p <0.05 for soil moisture prediction; RMSE = 0.48; R² = 0.74 and p <0.05 for seasonal f_d; Figure 6c,d). The model well predicted the strong seasonal variability of soil moisture, ranging from dry (~0.1 in summer) to saturated (~0.5 in winter) conditions (Figure 6c).

The distributed version of the ecohydrological model was then applied at the Fluminimaggiore basin comparing the modeled and observed runoff for the long 1924−2020 historical period (Figure 7). We first used the 1924−1975 period for calibration (Table 1 and Table 3; k = 20.8 d) and last the 1976−2020 period for validation. Both monthly and yearly runoff predictions were tested successfully (Figure 7; for monthly runoff: RMSE = 14.85 mm/month and R²= 0.94 and p < 0.001 in calibration, and RMSE = 17.76 mm/month and R² = 0.84 and p < 0.001 in validation; for yearly runoff: RMSE = 36.38 mm/year and R² = 0.98 and p < 0.001 in calibration, and RMSE = 48.51 mm/year and R² = 0.90 and p < 0.001 in validation), and the total predicted runoff of the whole observed period (~28 m) matched well to the total observed runoff (difference of 0.02%; Figure 7c). A sensitivity analysis of the model to the soil depths was also performed. We varied the basin soil depths uniformly with respect to the reference soil depth map (made by AGRIS) in a wide range (increasing up to +80% and decreasing up to −80%). The runoff predictions were mainly sensitive to the decrease in the soil depths (for values lower than −40%), due to the high increase in predicted runoff (which increased up to 58%, and RMSE increased up to 210 mm/y). Using more realistic scenarios of soil depths (variations from −40% to 40% respect to the reference soil depth map), the model results have low sensitivity to the soil depth variations (RMSE always lower than 53 mm/y), confirming the robustness of the model results.

The VDM predicted the annual LAI estimated from MODIS NDVI quite well. We successfully tested the basin-averaged LAI predictions (RMSE = 0.10, R²= 0.69, mean ratio = 0.99), and the mean LAI predictions of the forest (RMSE = 0.18, R²= 0.72, mean ratio = 0.99, Figure 8c).

We also tested the basin average yearly ET predictions with ET estimates derived from runoff and precipitation observations. In semiarid regions at an annual scale, the runoff is roughly equal to the precipitation minus ET from the basin water balance, assuming negligible changes of the basin soil water storage [11,124], so that annual ET estimates can be derived as the difference between observed runoff and observed precipitation. The annual ET model predictions well matched the annual ET estimates (RMSE = 3.93 mm/year; R² = 0.98, p = <0.05), confirming the model’s reliability at the basin scale.

At the basin scale, a nonsignificant negative trend of annual ET was predicted for the past (Mann−Kendall τ of −0.08 with p = 0.22, Theil−Sen β of −0.19 mm/year, Figure 8a), while the negative trend is significant for tree transpiration (Mann−Kendall τ of −0.16 with p = 0.02, Theil−Sen β of −0.26 mm/year, Figure 8b), and even more for tree LAI (Mann−Kendall τ of −0.24 with p < 10⁻³, Theil−Sen β of −0.005, Figure 8c)

Finally, we compared the changes in the past of runoff, ET and tree LAI predicted by the model using the actual observations of the Sardinian meteorological network, with those predicted using the HadGEM2-AO data, opportunely downscaled, as model atmospheric inputs. We compared the 1976−2000 period with the 1951−1975 period as proposed by Montaldo and Oren (2018) [11] in Sardinia. The use of the HadGEM2-AO atmospheric data of the past allowed us to sufficiently well predict the past decrease in the annual runoff (~20%), ET (~2%) and tree LAI (

Table 4. Past mean changes (1976−2000 with 1951−1975) of annual precipitation (P_y), air temperature (T_y), and predicted annual runoff (Q_y), evapotranspiration (ET_y), and tree LAI (LAI_t,y) using observations of the Sardinian meteorological network and the HadGEM2-AO meteorological data as model input.

	ΔPy (%)	ΔTy (%)	ΔQy (%)	ΔETy (%)	ΔLAIt,y (%)
From land-observations	−8.47	1.36	−20.73	−1.85	−1.39
HadGEM2-AO	−9.42	1.67	−16.18	−2.57	−2.07

).

3.2. Future Ecohydrological Scenarios

Future scenarios of HadGEM2-AO predicted a decrease in annual precipitation for all the representative concentration pathways, reaching a decrease of 14% for the 2061−2080 period and RCP 45, and 12% for the farther scenario (2081−2100) and the worst RCP 85 (Figure 9a) when compared with the recent 2000−2020 period. Annual air temperature is predicted to increase in future scenarios by up to 4 °C in RCP 85 and the 2081−2100 period (Figure 9b), impacting VPD, which is predicted to increase starting from 2041, by up to 20% in RCP 85 and in the 2081−2100 period (Figure 9c). Using the calibrated distributed ecohydrological model, the yearly runoff is predicted to further decrease in the future for most of the scenarios, mainly amplifying the rain changes (Figure 9d), and reaching the highest decrease of 26% for the RCP 45 configuration and the 2061−2080 period.

The warmer future scenarios (2022−2100) impacted the evapotranspiration, the trend of which is predicted to remain negative in the future, although not significantly (τ = −0.03, p = 0.68, β = −0.04 mm/year; Figure 8a), and even more so for tree transpiration, the decrease of which is predicted to be higher and significant in the future (τ= −0.22, p = 0.004, β = −0.23 mm/year; Figure 8b). The most alarming decreasing trend is for tree LAI (τ = −0.62, p <0.001, β = −0.01mm/year; Figure 8c), for which the average reduction is predicted to be 31% at the end of the predicted period when compared to its value at the start of the investigated period (Figure 8c).

4. Discussion

In several areas of the world, such as the Mediterranean basin, warmer climate conditions due to climate change are impacting water resources and forests [9,33,37,125], highlighting the need to carefully evaluate this impact using observations and accurate ecohydrological modeling, an essential tool for the prediction of the effects of future climate change scenarios. The Fluminimaggiore basin, in which the Marganai forest is located, is a meaningful case study because it is mainly natural (~97% of the basin) with low urbanization and has experienced negligible land use changes in the last century (Figure S1), providing the opportunity to evaluate the climate change impacts alone, unaltered by anthropogenic land use change effects [59]. The Marganai forest is a European Site of Community Importance (Natura 2000), whose habitats need to be preserved for the conservation of biodiversity in Europe, and a recent drought, in Summer 2017, stressed and affected the predominant tree species, the Quercus ilex, a common evergreen oak species in the Mediterranean region (e.g., [126,127]). Actually, the Fluminimaggiore basin was affected by changes in the main climatic factors, rainfall and air temperature in the past. The long database (almost one century) of the Sardinian forested basin allowed us to evaluate past climate changes. While annual rain decreased, but not significantly in the past, the winter precipitation decreased significantly with a Mann−Kendall τ of −0.15 (p < 0.05). Because the winter precipitation is the key contribution to the runoff in the Sardinian basins [8], the annual runoff trend followed the winter precipitation trend with even higher negative values, reaching a Mann−Kendall τ of −0.26 (p = 0.001), in agreement with the findings of Montaldo and Sarigu (2017) [10] for several Sardinian basins. Actually, runoff was very low in the 1990s, when water supply restrictions, including for domestic consumption, became frequent in Sardinia [14]. The decrease in the winter rain produced an alarming decrease in runoff and water resources supply, with potential consequences on forest and ecosystem sustainability, due to the decrease in the water reservoir in the deeper soil layers, as revealed by the significant negative trend of rock moisture (Mann−Kendall τ of −0.15, p = 0.03). In these water-limited environments, the rock moisture is a key water resource used by trees in dry seasons [47,49], and its decline can undermine the resilience of the trees to droughts.

At the same time, in the Sardinian basin, there was a significant positive trend of air temperature, with an increase of ~1 °C of mean annual air temperature in the last century (Figure 4b), which will be further amplified according to the future scenarios of the HadGEM2-AO climate model, which predicts an increase in air temperature by up to 4 °C in the worst scenario (RCP 85) and in the furthest period (2081–2100). We selected the HadGEM2-AO climate model from the 12 AOGCM2 of Flato et al. (2013) [88] as being the best when precipitation and air temperature predictions of the past are compared with observations of the Sardinian hydrographic service network, as proposed by Montaldo and Oren (2018) [11]. Then, we downscaled the original GCM to match the observed data statistical distribution to increase the correlation between the observed data and the simulated one taking into account the inter-variable relationships. Trees are suffering and will suffer even more in the warmer conditions in the future, with the predicted tree transpiration decreasing by 2.3 mm for the decade. The warmer conditions affected tree LAI, which decreased in the past, and is predicted to decrease still more in the future, reaching low values up to 2 (on average ~2.5) at the end of the twenty-first century, while tree LAI is ~4 currently (Figure 9c). The decrease in tree LAI is related more to the rise in air temperature (by 89%) than to the decline in precipitation, due to the impact of air temperature on VPD, which is a main climate control of the vegetation growth in these Mediterranean ecosystems [24,50]. Future scenarios predict a decrease in the mean annual precipitation (MAP) from ~810 mm in the recent 2001−2020 period to ~690−790 mm in the future (varying according to the representative concentration pathways scenarios of HadGEM2-AO), which could still be sufficient MAP values for sustaining the actual tree cover density in the basin. Sankaran et al. (2005) [34], Yang et al. (2016) [128], Axelsson and Hanan (2017) [35] and Corona and Montaldo (2020) [129] identified a MAP threshold of ~650−700 mm under which woody cover development becomes directly limited by rain in semi-arid ecosystems. In this Sardinian basin, future hydrological scenarios seem to still satisfy tree water needs to the limit: in the future scenarios, the rock moisture is predicted to further decrease on average by up to 10%, highlighting the achievement of a dangerous borderline condition, below which the underlying root-zone reservoir could no longer sustain tree water needs. In this respect, other Sardinian basins are characterized by lower MAP (e.g., ~40% of the Sardinian basins analyzed by Montaldo and Sarigu, 2017 [10], such as most of the basins in the Mediterranean area [130], and these basins could be in even more critical conditions.

Forest stress under climate change has been evaluated worldwide [32,33,36,37,125]. The negative impact of VDP increase on vegetation growth in the future was predicted on a global scale by Yuan et al. (2019) [23], and Williams et al. (2013) [51] predicted a strong increase in drought stress on the forests of the southwestern United States due to the VPD increase. Here, we predicted a reduction by almost half of the tree LAI in a protected natural forest in the Mediterranean basin mainly due to the VPD increase. The Marganai forest is suffering and will suffer even more from a reduction of tree cover that jeopardizes the forest’s sustainability and survival in the future, and its resilience to climate change. Forest management needs to deal with the impact of climate change on the forests [131], defining strategies for increasing the forest resilience and sustainability [36].

The proposed ecohydrological model well predicted the long runoff data series, including all the main hydrological and vegetation dynamics processes. The model was lightweight computationally, and it was possible to write the model code in MATLAB, although the spatially distributed approach is typically heavy for numerical model calculations. The model answered the need for parsimonious parametrization [52,74,75,76], mainly using a simplified approach for runoff estimation based on the Soil Conservation Service (SCS) method [81,83,84], which is widely used by other hydrological models (e.g., [86]) such as the SWAT model (e.g., [18,85], and the simplified approach of Montaldo et al. (2008) [68]) for vegetation dynamics modeling. At the same time, the model included the hydraulic lift effects due to tree root expansion in the rocky sublayer typical of the Mediterranean water-limited basins, which is a key term (on average ~50% of E_t) for properly predicting tree transpiration in semi-arid ecosystems [47,48,49]. The proposed ecohydrological model proved to be an effective tool for predicting climate change effects, which needed to be evaluated with long data series to avoid erroneous evaluations due to the limitation of the data extension [62,132].

In the Sardinian basin, runoff, a representative term for the water resources status, decreased in the past with a reduction by half in the last century (Figure 7a). The model predicted a further decrease in the runoff in the future (using the HadGEM2-AO scenarios) by up to an additional 26% (Figure 9d), highlighting a serious concern about the water resources supply in the future of the Sardinian forested basin. Sardinian basins have suffered historically from water scarcity [133], but the predicted future scenarios will be even worse. These results are consistent with those of Piras et al. (2014) [62] in another Sardinian basin, who used a short period (3 years) for model calibration, and the results of Fonseca and Santos (2019) [134] in a Portuguese Mediterranean basin. Based on HadGEM2-AO climate model scenarios, the future of the Fluminimaggiore basin will be characterized by less annual runoff but at the same time more extreme runoff events, as suggested by a frequency analysis of daily runoff extremes (Figure 10). The increases in the extreme precipitation in the future scenarios of the Sardinia region (by up to 10%, Figure S2) and the regions in the north Mediterranean Sea basin were also predicted by Trambay and Somot (2018) [135] and Marras et al. (2021) [18]. Rainfall extremes in the reference period were well estimated using the downscaled HadGEM2-AO data (Figure S2). Using the HadGEM2-AO (in the RCP 85 configuration) climate forcing, an increase in the frequency of maximum runoff was predicted, as can be deduced by the frequency analysis of the annual maximum daily runoff using a Generalized Extreme Value distribution, distinguishing between the past (1924−2021) and future (2022−2100) periods, and return periods from 5 to 100 years (Figure 10; the statistics of the Kolmogorov–Smirnov and Chi-squared tests are in

Table 5. Results of the statistical tests (Chi-square and Kolmogorov−Smirnov tests) of the Generalized Extreme Value distribution fitting to the annual maximum daily runoff time series of the Fluminimaggiore basin for the past (1924–2021) and the future climate scenario (2022–2100), predicted using HadGEM2-AO data in the RCP 85 configuration. Statistics of the Chi-square and Kolmogorov−Smirnov tests for normality, applied to the series, have 5% significance level.

Period	Kolmogorov–Smirnov	Χ2
1924−2021	0.0721 (0.1354)	4.2951 (5.991)
2022−2100	0.0783 (0.1505)	5.1175 (5.991)

). Daily runoff increased by up to 34% for return periods of 100 years for the future climate scenarios. We predicted that in the future, despite the reduction trend of the annual runoff, the increase in precipitation extremes will produce an increase in the possibility of extreme runoff events such as floods, with possible consequences for human life, structural damages, and river erosion, which may also be amplified due to the predicted reduction of tree cover in the basin [136,137]. These future scenarios will need to be carefully considered and included in the water resources management and planning of the Sardinian forested basins, using robust modeling approaches, such as our proposed distributed ecohydrological model.

5. Conclusions

A parsimonious model that couples a hydrological model and a vegetation dynamic model for long-term ecohydrological predictions in forested basins under water-limited conditions was successfully proposed. The model predicted the soil water balance in two root zone layers, one of which is an underlying rocky layer where tree roots penetrate, the runoff contribution using the simplified Soil Conservation Service (SCS) method, and the vegetation dynamics using a biomass balance approach for grass and tree LAI computations. Despite the model’s parsimony, it well predicted almost one century of runoff observations in a forested Sardinian basin, the Fluminimaggiore basin, characterized by low human impact and urbanization, which has been left in almost the same land use conditions over the past century.

The Sardinian basin is affected by climate change with a decrease in winter precipitation and, consequently, in annual runoff, which has reduced by half in the last century, and an increase of ~1 °C of the mean annual air temperature. In the Fluminimaggiore basin, the Marganai forest, which is a protected forest and a European Site of Community Importance (Natura 2000), has experienced alarming droughts due to the warmer climate conditions. Forest trees have suffered from these conditions, with tree LAI decreasing systematically mainly due to the air temperature and VPD rise at a rate of 0.1 hPa per decade.

Future scenarios (up to 2100) of the HadGEM2-AO climate model predict even warmer climate conditions, with less annual precipitation and higher air temperature and VPD, impacting both runoff and tree transpiration and LAI, with a predicted reduction of tree LAI by half in the next century. The annual runoff is predicted to decrease drastically in the worst scenarios (by up to 26%), and, at the same time, runoff extremes are predicted to increase severity, outlining future scenarios that will be drier and warmer with an increasing flood frequency. The future climate conditions will undermine the water resources and forest sustainability of the forested Sardinian basin, and their resilience to climate change.

The climate change impacts on water resources and forests need to be carefully considered in forest and water resources planning of the Sardinian forested basin, and, in general, of Mediterranean basins under water-limited conditions.