Regional patterns and drivers of modelled water flows along environmental, functional, and stand structure gradients in Spanish forests

The study of the water cycle in the forest at large scales, such as countries, is challenging due to the difficulty of correctly estimating forest water flows. Hydrological models can be coupled with extensive forest data sources, such as national forest inventories, to estimate the water flow of forests over large extents, but so far the studies conducted have not analysed the role of stand structure variables or the functional traits of the forest on predicted blue and green water flows in detail. In this study, we modelled the water balance of Spanish forests using stand structure and species data from forest inventories to understand the effects of climate, stand structure, and functional groups on blue water flows. We calculated blue water and green water flows and expressed them relative to received precipitation. Relative blue water flow was mainly concentrated in the wetter regions (Atlantic and alpine biomes) of Spain (around 25 %) in comparison with the Mediterranean biomes (10 %–20 %) and during the autumn–winter season. The leaf area index (LAI) of the forest stand is the most important predictor of relative blue water, exhibiting a negative effect until it reaches a plateau at higher levels (around 2.5–3). Deciduous forests showed a greater relative blue water flow than evergreen functional groups (25 %–35 % and 10 %–25 %, respectively) primarily due to leaf fall during the autumn–winter season. This study highlights how green water is decoupled from blue water; namely, blue water depends on winter and autumn precipitation, while green water depends on the spring and summer water demand and how the species' functional traits (deciduous vs. evergreen) can influence blue water production.

1 Introduction

Forests are one of the most important ecosystems on the planet and constitute a supplier of carbon and water for humanity (MEA, 2005). Even though water is the main determinant of vegetation worldwide (Wang-Erlandsson et al., 2022), the water cycle is particularly challenging to study because it is difficult to measure its components. The water cycle in forests depends on the precipitation that falls in the forest, which can be partitioned into the green water and the blue water (Caldwell et al., 2016; Schlesinger and Jasechko, 2014). The green water is the evapotranspiration flow from the forest to the atmosphere formed by the water transpired by plants together with the evaporation of intercepted precipitation by vegetation and the evaporation from the soil surface (Llorens et al., 2011). Climate, specific traits, stand forest structure, and their interactions determine the partitioning of precipitation between green water and blue water. Forest evapotranspiration is determined by local climate conditions, such as air temperature, radiation, vapour pressure deficit, or soil water availability (Granier et al., 2000). Moreover, plant transpiration is different among species depending on their specific traits, as leaf area and xylem traits determine the efficiency of water transport or leaf phenology (evergreen and deciduous) (Ford et al., 2011). Stand forest structure is also key in rainfall partitioning through its relationship with stand leaf area index (LAI; total area of leaves of the canopy per unit horizontal ground area), which is a key variable determining transpiration (Granier et al., 2000) and/or the interception of vegetation. On the other hand, blue water can be split into two flows: the surface runoff and the surplus of groundwater by the system and available downstream.

Green water reduces the amount of blue water produced at a local scale but recycles the water through evapotranspiration at a global scale (Ellison et al., 2012). Green water flow is usually greater than blue water's in forests and grasslands in comparison with croplands or wetlands due to the greater transpiration of the former (Oki and Kanae, 2006). Green water represents around 40 %–70 % of water flow in temperate and boreal forests (Schlesinger and Jasechko, 2014), and even 90 %–100 % in drier environments or years, including Mediterranean forests (Ungar et al., 2013; Campos et al., 2016; Qubaja et al., 2020). The combination of different rainfall dynamics, transpiration, stand structure, and specific life traits can modify the partitioning between green and blue water. For instance, water yield can be affected by the temporal distribution of rainfall events (seasonality and torrentiality). Relative water interception by trees decreases when rainfall increases, being greater in arid environments compared to humid forests (Levia and Frost, 2003). In Mediterranean forests, the interception is greater than in temperate forests due to the smaller amount of rainfall and intensity and the higher evaporation caused by the higher temperature (Limousin et al., 2008). However, evaporation of intercepted water is lower during late-summer intense convective storms since the rainfall intensity is very high and lasts a short time (torrentiality) (Llorens and Domingo, 2007). Physiological and phenological differences among species also modulate water flows. Evergreen Mediterranean trees can reduce their transpiration in summer (water saver), but deciduous trees increase it due to higher water demand (water spender) (Baquedano and Castillo, 2006; Klein, 2014; Link et al., 2014; McDowell et al., 2008). Therefore, Mediterranean forests show transpiration seasonality, with a reduction in stomatal conductance during the summer drought, while the transpiration in temperate forests follows leaf phenology more closely (Llorens et al., 2011).

Water balance in forests has been studied mainly at local and landscape scales, like in watersheds (Caldwell et al., 2016; Guzha et al., 2018; Schwärzel et al., 2020) or forest stands (Benyon et al., 2017; Simonin et al., 2007). Analyses at regional scales have also been carried out (Hoek van Dijke et al., 2022; Mastrotheodoros et al., 2020; Sun et al., 2005), analysing the water flow variation by topoclimate or landscape vegetation cover with hydrological models. The study of water fluxes in forest ecosystems is an important challenge for forest management. Ecohydrological simulations have been shown to be a valuable tool to understanding water fluxes at different spatial scales, complementing the limited field data available (Hoek van Dijke et al., 2022; Mastrotheodoros et al., 2020). The effect of stand structure variables on the water flow have been studied at a local scale with field data (Benyon et al., 2017; Simonin et al., 2007). Nevertheless, landscape-level ecohydrological simulations have been most often carried out without incorporating the role of the stand structure variables or the differences in species and functional composition (Hoek van Dijke et al., 2022; Mastrotheodoros et al., 2020).

In our study, we analysed the spatial pattern and partitioning of blue and green water in Spain using ecohydrological simulations with the MEDFATE water balance model (De Cáceres et al., 2023, 2015), which uses detailed field stand structure and species composition derived from forest inventory data. We simulated water flows and their green and blue water components at a daily resolution, and we studied the spatial distribution of relative blue water and its seasonality. Spain is an adequate region to study the partitioning between blue and green water because it shows contrasting climate gradients in terms of temperature, amount of rainfall, and seasonal distribution. Importantly, it includes a contrasted climate between the Mediterranean area – with hot temperatures, low precipitation, and seasonal drought – and the temperate area – with warm temperatures and a higher amount of rainfall and without seasonal drought (Rivas-Martínez et al., 2011). In addition, peninsular Spain harbours both temperate and Mediterranean forests that span a broad gradient of forest types and stand structures that can modify the partitioning between blue and green water. The questions addressed by the study are as follows: (1) which are the spatio-temporal patterns of blue and green water throughout the Spanish Peninsula and among climate subregions? (2) How does this partitioning between blue and green water vary among contrasting forest functional groups? (3) Which are the main climatic and forest structural drivers of this partitioning between blue and green water?

2 Material and methods

2.1 Study area

The studied area spans the forested areas of Spain, including the Iberian Peninsula and the Balearic Islands, but excluding the Canary Islands. This encompasses two climatic domains: the temperate climate with Atlantic Sea influence in the north and west parts of the Iberian Peninsula and the Mediterranean climate on the rest of the territory. The temperate climate is wetter and colder, while the Mediterranean climate is hotter and drier. We have classified the Iberian territory into six biomes (see map in Fig. S1 in the Supplement) according to the Iberian climate classification of Allué Andrade et al. (1990): (i) arid (9.0 % of the Spanish surface), with warm winters and summers and very low precipitation (< 400 mm); (ii) temperate Mediterranean (34.8 %), which is hot in winter and summer with low precipitation (400–600 mm); (iii) continental Mediterranean (30.9 %), which is cold in winter and hot in summer with low precipitation (400-600 mm); (iv) sub-Mediterranean (11.0 %); which is cold in winter and hot in summer with wetter precipitation (about 800 mm); (v) Atlantic (10.5 %), characterized by a mild winter and summer and high precipitation (> 1000 mm); and (vi) alpine (3.8 %), which is very cold in winter and cool in summer and characterized by high precipitation (> 1000 mm).

2.2 Forest inventory data

We used the permanent field plots from the Spanish national forest inventories to characterize the predominant species and stand structure variables across the study area. The Spanish national forest inventories (SFIs) are distributed with a systematic survey along the forested areas with a density of ∼ 1 plot km⁻². Specifically, we used the third inventory (SFI3) carried out in 1997–2008, which is the most recent and complete survey available to date. Spanish permanent field plots have four concentric circular subplots with radii of 5, 10, 15, and 25 m. Trees were identified and measured within variable circular size plots (5 m radius for trees with diameter at breast height (dbh) ≥ 7.5 cm, and trees with 2.5 ≤ dbh ≤ 7.5 cm are counted; 10 m radius for trees with dbh ≥ 12.5 cm; 15 m radius for trees with dbh ≥ 22.5 cm; and 25 m radius for trees with dbh ≥ 42.5 cm). Moreover, the shrub species or genus within the 10 m radius subplot was measured with its corresponding percent cover value and average height.

We focused on the field plots with > 70 % of plot basal area of the most predominant tree native species (Table S1) and an overall basal area of > 3 m² ha⁻¹ (to ensure that very sparse woodlands, which can hardly be considered a forest, were excluded), resulting in a total number of plots of 32 514. Plots were classified according to the dominant species in terms of basal area into five functional groups: temperate (e.g. Fagus sylvatica) and Mediterranean deciduous forests (e.g. Quercus pyrenaica), temperate (e.g. Pinus sylvestris) and Mediterranean coniferous forests (e.g. Pinus pinea), and sclerophylls (e.g. Quercus ilex) (see Fig. S3 for the SFI3 map of the plots selected and Table S1 for the classification of the tree species into these five groups).

2.3 Environmental variables

We downloaded the daily data of total precipitation and maximum, mean, and minimum temperature for the 5 years before and after the date of survey for every plot (that is, 10 years per plot) from E-OBS data at a 1 km horizontal resolution (Moreno and Hasenauer, 2016); ftp://palantir.boku.ac.at/Public/ClimateData/ (last access: 2 February 2023). Topographic variables of slope, aspect, and altitude were derived from a digital elevation model (DEM) at 25 m² of the Spanish National Orthophoto Program (PNOA) through the Instituto Geográfico Nacional (IGN) website (https://centrodedescargas.cnig.es/CentroDescargas/index.jsp#, last access: 14 February 2022). We used these topographic variables to obtain daily values of the relative humidity, solar radiation, and potential evapotranspiration (PET) per plot through the meteoland R package (De Cáceres et al., 2018). We characterized plot-level climatic conditions using three climatic indexes: climatic moisture (mean annual precipitation / mean annual PET), continentality (mean temperature of the hottest month with the mean temperature of the coldest month subtracted), and precipitation seasonality (standard deviation of monthly precipitation / mean of monthly precipitation) multiplied by 100).

A summary of the main stand structure and climate characteristic of the different functional groups and biomes is shown in Table S2.

2.4 MEDFATE model

The MEDFATE model (version 2.9.3) has been designed to simulate plant water balances and soil in structurally and compositionally heterogeneous forest stands (De Cáceres et al., 2015, 2023). MEDFATE uses daily weather as input and most processes are simulated at daily time steps. The aboveground stand structure is represented in terms of total crown ratio (CR), height (H), and leaf area index (LAI) of a set of woody plant cohorts. The soil is represented by vertical layers with different hydraulic properties, and each cohort may have a different root distribution specified using the depth corresponding to cumulative 50 % and 95 % of fine roots (Z₅₀ and Z₉₅, respectively).

Each day, the model updates the leaf area for (semi-) deciduous vegetation based on a simple phenological model that dictates and leaf drop. This model relies on the degree days required for leaf budburst (SGDD) parameter, representing the degree days needed for budburst (assuming evergreen plants maintain a consistent leaf area throughout the simulation). Subsequently, the model revises the light attenuation within the canopy, adhering to the Beer–Lambert model, as well as the canopy's water storage capacity, which signifies the minimum water quantity required to saturate the canopy. Following the canopy status update, the model addresses the input of water from precipitation. Prior to augmenting the soil layers' water content, the model initially deducts the water loss from rainfall due to interception and surface runoff. Rainfall interception loss is estimated using the simplified version of the Gash model (Gash et al., 1995), while runoff is calculated according to the United States Department of Agriculture Soil Conservation Service (USDA SCS) curve number method (Boughton, 1989). Processes related to lateral water transfer are omitted from consideration. Soil water storage capacity and water potential are derived from soil texture through pedotransfer functions (Saxton et al., 1986). When replenishing a specific soil layer, a portion of the water is assumed to directly infiltrate the layer below, as determined by microporosity (Granier et al., 1999). The water percolating from the deepest layer is presumed to be lost through deep drainage.

To assess plant transpiration, the model initially calculates a distinct estimation of the maximum transpiration for the entire plant community (including trees and shrubs) – that is, without considering the soil water deficiency. This calculation is done for each taxon and takes into account the atmospheric evaporative demand. It involves two taxon-specific parameters: $a_{T_{\max }}$ and $b_{T_{\max }}$ parameters. The estimation of maximum transpiration for the entire stand (E_max,stand(i)), excluding considerations for soil water deficit, relies on the daily Penman potential evapotranspiration (PET) and an empirical relationship established by Granier et al. (1999). However, MEDFATE modified the Granier equation with $a_{T_{\max }}$ and $b_{T_{\max }}$:

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{E_{\max ,\mathrm{stand}}\left(i\right)}{\mathrm{PET}}}=a_{T_{\max }}\cdot \mathrm{LAI}+b_{T_{\max }}\cdot {\mathrm{LAI}}^{\mathrm{2}},\end{array}$$

where $a_{T_{\max }}$ and $b_{T_{\max }}$ represent species-specific parameters (De Cáceres et al., 2023). Assuming reliable species-specific estimates are accessible for $a_{T_{\max }}$ and $b_{T_{\max }}$, the equation can be applied to calculate E_max,stand(i). Once E_max,stand(i) is determined for each species in the stand, the portion of shortwave radiation absorbed is employed to estimate its maximum transpiration (E_max(i)) from E_max,stand(i) (Korol et al., 1995).

Moreover, the cohort's transpiration is influenced by the vertical distribution of fine roots, the soil moisture profile, and two parameters specific to the taxon: Ψ extract and c extract. These parameters represent the soil water potential at which 50 % of the maximum transpiration occurs and the slope of a Weibull function that regulates the rate of transpiration decline, respectively. The plant's water status is expressed as plant water potential, denoted as Ψ plant, which is defined as the “average” of the soil water potential in the rhizosphere. The model monitors the effects of drought by assessing the proportion of hydraulic conductance lost due to stem cavitation, referred to as PLC. Increases in PLC occur whenever Ψ plant decreases, following a xylem vulnerability curve characterized by the parameters VC_stemc and VC_stemd. PLC sets limits on actual transpiration rates and does not decrease even when Ψ plant increases. It is worth noting that the effects of cavitation can only be reversed, i.e. PLC can be reduced, through the formation of new sapwood, as outlined by Choat et al. (2018).

2.5 Parameter estimation

Data from forest inventory plots included tree height (H) and tree diameter at breast height, which were used to obtain estimates of foliar biomass (henceforth leaf area, after multiplying by the specific leaf area, SLA) and crown ratio (CR) via species-specific allometries (see Tables S1–S3 of De Cáceres et al., 2023, for more details). In shrubs, foliar biomass is calculated from shrub height via species-specific allometries. In the model, SLA (specific leaf area; in mm² mg⁻¹), the ratio of leaf area to leaf dry mass, is constant for every species (De Cáceres et al., 2023). The leaf area index (LAI; in m² m⁻²) was calculated from the foliar biomass (in kg m⁻²) using a specific leaf area coefficient (SLA; in m² kg⁻¹) that is species-specific (LAI = foliar biomass ⋅ SLA). Taxon-specific parameter details are shown in the supplement of De Cáceres et al. (2023). Soil data of each forest inventory plot were extracted from the SoilGrids database (Hengl et al., 2017). For all plots, four soil layers down to a total depth of 4 m were initially considered, but the deepest layers were merged into a rocky layer (95 % of rocks) following the depth of the R horizon. A monotonous increase in rock fragment content across soil layers from the surface to the rocky layer was defined based on surface stoniness classes determined in SFI3 plot surveys.

2.6 MEDFATE simulations

MEDFATE was run on each selected SFI3 plot using daily weather data (temperature and precipitation, PET, and radiation and relative humidity) corresponding to a 10-year period centred on the year of the SFI3 sampling (1997–2008). We calculated the blue water as the sum of the runoff and deep water and the green water as the evapotranspiration (that included the sum of the transpiration, interception, and soil evaporation).

2.7 Model evaluation

MEDFATE predictions have already been evaluated at the forest stand scale in terms of soil moisture dynamics, plant transpiration, and water status in Mediterranean forests (De Cáceres et al., 2021, 2015). Given the focus of the present work, we evaluated regional-scale patterns of green and blue water predicted by MEDFATE against those produced by alternative methods. First, we did a comparison between the results for average blue water of MEDFATE with the average blue water from the Precipitation–Runoff Integrated Model (SIMPAL; the acronym is in Spanish) of the Spanish government. SIMPAL is calibrated with stations that measure streamflows across Spain, and it is interpolated at a 1 km² resolution for the whole country (Estrela et al., 2012). Second, we compared green water patterns predicted by MEDFATE with those of GLEAM, which is derived from satellite data and covers the world at a resolution of 0.25° (Martens et al., 2017). We used the GLEAM v3.8a that defines the evapotranspiration as the sum of transpiration (from short and tall vegetation), interception (from tall vegetation), soil and open water evaporation, and snow sublimation.

2.8 Statistical analysis

Our analyses focused on the relative blue water, the ratio between simulated annual mean blue water, and the total precipitation as response variable. The variation in relative blue water across biomes and functional groups was analysed using beta regression models (Ferrari and Cribari-Neto, 2004). We used beta regressions because the response variable was a proportion. A post hoc Tukey's test was applied for pairwise comparisons between groups (i.e. biomes and functional groups).

The effect and importance of climate, topographic, and stand structure variables in the partitioning of blue/green water of every functional group were analysed with XGBoost regression. XGBoost is a machine learning system that builds a sequential series of shallow regression trees with a gradient boosting technique (Chen and Guestrin, 2016). XGBoost regression was used because it allows us to analyse linear and non-linear responses and correlated variables. A regression tree is trained by splitting the input dataset into increasingly homogeneous subsets at each decision node and choosing the split that maximizes the distance between different terminal nodes. As predictors, 14 different climatic variables calculated as the mean annual values of the 10 years of climate data for every plot (mean, maximum, and minimum temperature; solar radiation; and moisture and continentality indexes) and the mean of seasonal variables (precipitation of spring, winter, summer, and autumn and precipitation seasonality index) were used. A topographic variable (slope) and two stand structure variables (basal area and LAI of the plot) were also included. Basal area serves as an indicator of canopy cover, as it is calculated by summing the trunk cross-sectional area at breast height of all trees per hectare, and higher basal area values can result in increased interception. LAI determines the transpiration of the trees and shrubs through the leaf surface. LAI of the plot was calculated from the sum of the LAI of trees and shrubs. One model was computed for each of the five functional groups. The models were tuned for finding the best hyperparameters and they were evaluated using a k-fold cross-validation procedure (10-fold with five repetitions), following a stratified random sample into k subsamples of the dataset. A forward feature selection method for selecting the predictor variables with a spatial cross-validation (Meyer et al., 2018) was used. The models were built for each predictor pair combination and the best was selected. The remaining predictors were tested by adding them to the best combination. This procedure continued until neither of the remaining predictors resulted in an improvement of the model. Variable importance in the retained variables was evaluated by computing the fractional contribution of each feature to the model based on the total gain of these variable splits (gain), the relative number of times a feature is used in trees (frequency), and the relative value of the feature observation (cover). We estimated R² to test the accuracy of the final models with 1 − (Sres/Stot), Sres being the sum of the squared differences between the observed values and the predicted values and Stot the total sum of squares, which is the sum of the squared differences between the observed values and the mean of the observed values. All statistical analyses were conducted with R v.4.1.1., with the packages xgboost (Chen et al., 2016) and caret (Kuhn, 2008).

3 Results

3.1 Model evaluation

The comparison indicated that the blue water of SIMPAL and the green water of GLEAM followed the same regional patterns as MEDFATE (Fig. 1). SIMPAL reaches higher values of blue water in the Atlantic biome and lower values in west temperate Mediterranean than MEDFATE. In opposition, GLEAM evapotranspiration is higher in the west temperate Mediterranean and lower in the Atlantic forest than MEDFATE. MEDFATE models are at a stand scale, whereas SIMPAL and GLEAM models are at a regional scale. At a stand scale, the evapotranspiration is high in the Atlantic forest, and then the blue water is lower. The west temperate Mediterranean is characterized by forests of low basal area and LAI (Table S2; Fig. S4). Therefore, the evapotranspiration is lower and blue water higher than surrounding forest. In other words, the differences observed seem to arise from the fact that MEDFATE takes into account the stand structural characteristics, which are difficult to represent in models based on interpolation (SIMPAL) or remote sensing (GLEAM).

Figure 1 Average evapotranspiration (ET) maps for SFI3 plots according to GLEAM (a) and MEDFATE (b) (Spearman correlation = 0.62; R²= 0.218) and average blue water maps for SFI3 plots according to SIMPAL (c) and MEDFATE (d) (Spearman correlation = 0.70; R²= 0.354).

3.2 Spatial patterns of blue water in the Spanish biomes

Forests with the highest absolute amount of blue water (> 500 mm) were mostly concentrated on the coast of the Atlantic Sea and in the Pyrenees (Fig. 2a). The south of the Spanish Peninsula in the temperate Mediterranean biome and some areas inside the continental Mediterranean biome had values between 300 and 500 mm. The rest of the country had lower blue water values (< 300 mm), mainly in the eastern and Mediterranean regions. Relative blue water showed a very similar spatial pattern, with the highest values located on the coast of the Atlantic Sea and the Pyrenees (Fig. 2b). In the southwest of Spain and in the Mediterranean mountain regions, the relative blue water was also high (> 50 %). Forests in the eastern Mediterranean region had the lowest values. The arid and continental Mediterranean biomes showed a very large percentage of plots (> 70 %) with low relative blue water (< 15 %), whereas the Atlantic and alpine biomes showed higher plots (> 50 %) with high blue water (> 30 %) (Table S3). Temperate Mediterranean and sub-Mediterranean biomes showed more intermediate values, although the temperate biome had the percentage of plots with the highest blue water (> 60 %) among all biomes. The results of the beta regression model and post hoc Tukey's test indicated that there were differences among biomes in the percentage of blue water (Table 1; Fig. 3; see Figs. S5–S6 for interception and soil evaporation differences): relative blue water was higher in the Atlantic and alpine biomes (mean over 25 %) followed by the temperate Mediterranean and sub-Mediterranean biomes, with the arid and continental Mediterranean biomes showing the lowest values (mean of 10 %–15 %). The absolute values showed a similar pattern to those of the relative blue (see Fig. S7).

Figure 2 Total (a) and relative (b) blue water in the Spanish Peninsula and Balearic Islands. In the grey area, no forest covers are shown. Maps were realized with a rasterization of the SFI3 plot results at 1 km².

Figure 3 Percentage of blue water in the different biomes. Different letters denote significant differences among biomes (p < 0.01) after Tukey's test. ARID: arid, TEME: temperate Mediterranean, COME: continental Mediterranean, SUME: sub-Mediterranean, ATLA: Atlantic, ALPI: alpine.

Table 1 Results of the beta regression models regarding relative blue water (%) for the different biomes, indicating the estimate, standard deviation (SD), statistic, and p value of each of them. A positive value indicates that an increment occurred in the biome with respect to the observed values in the intercept.

3.3 Patterns of blue water among forest functional groups

Beta regressions models and post hoc Tukey's test indicated that there were also differences among functional groups in the percentage of blue water (Table 2; Fig. 4; see Figs. S8–S9 for interception and soil evaporation differences). Relative blue water was high in temperate deciduous forests (35 %) followed by temperate coniferous forests, Mediterranean deciduous forests, and sclerophylls (20 %–25 %). Mediterranean coniferous forests had the lowest value (close to 10 %). The absolute values showed a similar pattern to that of the relative blue with the temperate deciduous with the highest values (see Fig. S10).

Figure 4 Percentage of blue water of the different functional groups of species. Different letters denote significant differences among forest types (p < 0.01) after Tukey's test. TEDE: temperate deciduous, MEDE: Mediterranean deciduous, TECO: temperate coniferous, MECO: Mediterranean coniferous, SCLE: sclerophylls.

Table 2 Results of the beta regression models regarding relative blue water (%) for the different forest functional groups, indicating the estimate, standard deviation (SD), statistic, and p value of each of them. A positive value indicates that an increment occurred in the functional group with respect to the observed values in the intercept.

The analysis of seasonality using monthly aggregated data showed that the amount of blue water in forests of the different functional groups strongly varied throughout the year. In the temperate deciduous forests, total blue water was higher than total evapotranspiration in autumn and winter (Fig. 5a), whereas in other functional groups, the values of the two flows were similar throughout the year, except in spring–summer, when the green water was higher. In the Mediterranean coniferous forests, evapotranspiration was higher than blue water all year long. Relative flows showed that blue water was at its minimum in the summer season in all groups and higher than evapotranspiration only in autumn–winter in temperate deciduous forests (Fig. 5b).

Figure 5 Total (a) and relative (b) green water and blue water and total precipitation by different functional groups throughout the year. Note that evapotranspiration can be higher than precipitation during summer months due to plant transpiration and bare-soil evaporation of water available in the soil. Relative evapotranspiration has been truncated to 100 % of precipitation with the aim to improve the visualization of water export. A line marks the seasons in the months.

3.4 Main drivers of the partitioning between blue and green water

XGBoost models showed high accuracy, with R² > 75 % in most functional groups (Table S4). These analyses revealed that LAI was the most important predictor of relative blue water in all functional groups (Fig. 6). The second and third most important variables were the climatic moisture index and winter precipitation, except for Mediterranean deciduous forests, for which it was autumn precipitation (see Figs. S11–S12 for spatial variation in climatic moisture and winter precipitation). The other variables retained in the models had low significance.

Figure 6 Average importance from gain, cover, and frequency values for every predictor variable retained in the XGBoost models for the five functional groups.

Figure 7 Partial dependence plots for the three main predictor variables(a, LAI; b, climatic moisture index; c, winter precipitation) and their effect on blue water percentage for the different functional groups.

Partial dependence plots showed that LAI had a strong and negative effect on relative blue water (Fig. 7). Nevertheless, this effect was different depending on the functional group. The temperate deciduous forest showed the highest percentage of blue water for the whole LAI gradient, being always at a value of over 20 %. The other forest types showed similar patterns among them, with a larger drop until LAI was < 3 and values of around 10 % afterwards. The only exception was the group of Mediterranean conifers, which showed lower percentages of blue water than the other groups for the LAI gradient. The climatic moisture index showed a positive effect in all groups. The three angiosperm groups had a similar slope and showed higher values than the two gymnosperm groups (Fig. 7b). Winter precipitation had a positive effect on all functional groups except on the Mediterranean deciduous forests (Fig. 7c). Temperate deciduous forests showed higher values through the whole winter precipitation gradient.

4 Discussion

4.1 Spatio-temporal patterns

Blue water spatial patterns showed that most Spanish forests have a proportion of the precipitation below 15 % as blue water (Fig. 3). This pattern was expected since most of the Spanish forests are in the Mediterranean climate, and in these dry environments, the evaporative demand is high and the majority of water evapotranspired or was intercepted by the canopy (Ungar et al., 2013). In contrast, in the Atlantic and alpine biomes, the proportion of rainfall as blue water can be higher than 25 % because these biomes are wetter and colder than the Mediterranean, and, thus, the evapotranspiration is lower and precipitation is higher (Kosugi and Katsuyama, 2007; Mastrotheodoros et al., 2020). Moreover, in some mountain areas in the continental Mediterranean biomes (e.g. Sierra Nevada and Sierra de la Demanda), the amount and proportion of blue water are also high (Fig. 2a, b). These regions are wetter than the rest of Mediterranean biomes and an increase in precipitation determines an increase in blue water (Helman et al., 2017). Interestingly, in other areas of southern Spain, in the temperate Mediterranean biome, the proportion of blue water can also be higher than 20 % (e.g. Sierra Morena; Table S3), not because they have higher rainfall but because this biome has very open forests, with LAI values lower than in the rest of biomes (Table S1). Transpiration is the main component of the green water in arid environments and decreases in wetter regions, where the interception increases their importance (see Fig. S13). This result is concordant with the global patterns of green water flux, where the transpiration is bigger in arid environments (Good et al., 2017).

The analysis of seasonal dynamics by forest type showed that in forests of all functional groups, there is more blue water during autumn and winter months. In contrast, the evapotranspiration was higher in spring–summer. In Mediterranean climates, there is a decoupling between rainfall and temperature (Baldocchi and Xu, 2007), with higher rainfall in spring and autumn and higher temperatures in summer. For this reason, plants have more photosynthetic activity and growth in spring (when rainfall is high) and more transpiration in summer (when temperature is the highest) (Baquedano and Castillo, 2006; Link et al., 2014), when trees use the groundwater stored in the soil (Jost et al., 2005). Water loss by runoff or deep drainage is higher in autumn–winter because water is not fully exploited by the forest in this lower-photosynthetic period, and, thus, green water is made lower than would be possible by falling rain (Helman et al., 2017; Williams et al., 2012). According to our results (Fig. 5), blue water is concentrated in this period of lower or dormant vegetation activity (i.e. late autumn and winter).

4.2 Differences in blue water between forest functional groups

The analysis of the percentage of blue water of the different functional groups of species showed that temperate deciduous forests had the highest value followed by Mediterranean deciduous forests (Fig. 4). The rest of forest types had lower values (Fig. 4), mainly in the case of Mediterranean coniferous forests, which have shown high green water of almost 90 %, as also observed in studies at the plot level (Ungar et al., 2013). The seasonal dynamics of these flows in the different forest types show that deciduous forests have more blue water in autumn and winter months than evergreen forests (Fig. 5). This indicates that the main factor that explains the differences in relative blue water between forest types is the shedding of leaves (or its lack of shedding) during the cold part of the year. Although deciduous trees are mainly water-spending since their transpiration is higher than that of evergreens (Baldocchi, 2020), they shed their leaves in autumn, and their interception and evapotranspiration are lower than those of evergreens in autumn/winter (Fig. 5). Moreover, due to their shape, coniferous needles intercept the rainfall more than the deciduous ones for the same area of leaves (Carlyle-Moses, 2004). Although green water is higher in deciduous than in evergreen forests in summer (Fig. 5), total annual blue water seems to be more determined by leaf phenology, which conditions the destination of the autumn–winter precipitation and determines the differences between deciduous and evergreen forests. Also, climate partial plots showed that at the same level of climatic moisture index and winter precipitation, deciduous forests have greater relative blue water than coniferous ones and sclerophylls (Fig. 7b–c).

However, once the differences between deciduous and evergreen forests are highlighted, the next aspect that stands out when comparing the different types of forests is that for each category of leaf phenology, the proportion of blue water depends mainly on climate. Thus, temperate deciduous forests have more relative blue water than Mediterranean deciduous ones (Fig. 4). Temperate deciduous forests, which are the ones exporting the most blue water in relation to their precipitation, are mainly concentrated in the Atlantic and alpine biomes, where precipitation is higher and transpiration is lower, while Mediterranean deciduous forests are concentrated in Mediterranean biomes, where precipitation is lower and transpiration is higher (Table S2). The same pattern happens in evergreen forests: temperate coniferous forests have a percentage of blue water that is higher than Mediterranean coniferous forests because precipitation in the biomes where they are more abundant is also higher than in the Mediterranean. The lowest value of blue water in Mediterranean coniferous forests is particularly relevant because this is the most abundant forest type in Spain, and it is predominant in Mediterranean biomes (Table S1). Finally, sclerophyll forests produce more blue water than Mediterranean coniferous forests although both are situated in Mediterranean biomes and are evergreen. Sclerophyll forests have higher transpiration than Mediterranean coniferous ones (del Campo et al., 2019; Sánchez-Costa et al., 2015), but the total LAI of their stands is lower (Table S2), which determines the different pattern obtained (as explained in the next section).

4.3 Differences in stand structure determining blue water

The analysis of the drivers affecting relative blue water showed that LAI is the most important predictor of blue water (Fig. 6). In our simulations, relative blue water was often very high (> 60 %) when LAI was < 1. Transpiration and interception are higher when the LAI of the stand is higher, and, therefore, the water is transpired and intercepted by leaves. However, the temperate deciduous forests showed higher blue water for the same LAI values. Deciduous trees present higher evapotranspiration than evergreens at the same level of LAI (Baldocchi et al., 2010), but, as described previously, the LAI effect is modulated by leaf phenology since blue water in these forests is concentrated in winter. When LAI is very low, evergreen forests can show higher relative blue water. For instance, the temperate Mediterranean biome showed a higher percentage of plots with a higher relative blue water than Atlantic or alpine biomes (Table S3) because the mean LAI of their forests is lower than in the rest of biomes (see Table S1). This points out the importance of stand structure and composition since in a drier environment, LAI values higher than 2–3 m² m⁻² can strongly reduce the relative blue water. This effect is higher if the forest is evergreen, which is the predominant type in the Mediterranean. Relative blue water of Mediterranean deciduous forests is higher than that of evergreens because a decrease in LAI in evergreen forests decreases tree transpiration but increases evaporation from the soil surface, which is high in the Mediterranean climate even in winter (Raz-Yaseef et al., 2010). The management of blue water in forests has been focused in modifying the basal area of the stand (del Campo et al., 2014). Nevertheless, in this study, basal area was a limited predictor of the blue water in comparison with LAI. LAI of the stand increases with basal area until a threshold, and then it remains stable (Fig. S14). Although the increase in basal area over time reduces blue water (Caldwell et al., 2016), a plateau is expected in the long term, since LAI will not increase more after a maximum is reached. For this reason, blue water did not reduce excessively when LAI was higher than 3 (Fig. 7a). The reduction in basal area in the forest through thinning is a recurrent way of increasing blue water (Ameztegui et al., 2017; del Campo et al., 2018; Simon and Ameztegui, 2023), but this non-linear relationship between basal area and LAI implies that only heavy management for reducing basal area would be effective for reducing the LAI of the stand and, thus, increasing blue water.

4.4 Limitations

Model simulations to estimate water flows are not the same as observing them and can lead to artificial results. In this sense, it is necessary to continue improving the ecohydrological model design parametrization to improve their accuracy, including the completion of traits of many species for the parametrization, a better soil database, or the inclusion of lateral water flows (in the case of the MEDFATE model). Model evaluation is another key aspect to increasing realism and usefulness of predictions. The model MEDFATE has been validated in Mediterranean forests on a local scale with field ground data, but it has not been validated on a country scale with field data like eddy flux or sap flow data, and it might have lower accuracy in the Atlantic or alpine biome or with some tree species. Nevertheless, the comparison between MEDFATE results and other land surface products showed similar regional patterns. Despite accepting these limitations, the lack of field data at bigger scales encourages the use of models to study water flows on a country or continental scale.

5 Conclusion

Ecohydrological models like MEDFATE allow us to model different species and stand structures, improving our insights into the role of forests in water flows. Thus, we can split the fluxes of water at a stand level and understand the role of species, climate, and stand structure variables in the different water flows: transpiration, interception, runoff, etc. Moreover, the incorporation of forest inventories into these models allows us to conduct analyses at large spatial scales across multiple climates, forest typologies, and species compositions. In the current context of global change, the knowledge of water fluxes in forests is essential to being able to define strategies to improve the water budget, especially in dry environments like the Mediterranean, where water is scarce and the previsions are in line with increasing drought periods. Thus, in the context of global change, we are suffering (Hoerling et al., 2012; Lionello et al., 2014), the low amount and proportion of blue water in a major part of Spain is an important problem in the water supply for the agriculture or human consumption in the region (Ellison et al., 2012). The consequences that this lack of blue water might have on the entire silvoagropastoral system in the whole region, aggravated by the harsher conditions expected in the future, should make us reflect on what kind of forest management could be more effective in the next few decades.