Mapping the suitability of groundwater-dependent vegetation in a semi-arid Mediterranean area

Mapping the suitability of groundwater-dependent vegetation in semi-arid Mediterranean areas is fundamental for the sustainable management of groundwater resources and groundwater-dependent ecosystems (GDEs) under the risks of climate change scenarios. For the present study the distribution of deep-rooted woody species in southern Portugal was modeled using climatic, hydrological and topographic environmental variables. To do so, Quercus suber, Quercus ilex and Pinus pinea were used as proxy species to represent the groundwater-dependent vegetation (GDV). Model fitting was performed between the proxy species Kernel density and the selected environmental predictors using (1) a simple linear model and (2) a geographically weighted regression (GWR) to account for autocorrelation of the spatial data and residuals. When comparing the results of both models, the GWR modeling results showed improved goodness of fit as opposed to the simple linear model. Climatic indices were the main drivers of GDV density, followed by a much lower influence by groundwater depth, drainage density and slope. Groundwater depth did not appear to be as pertinent in the model as initially expected, accounting only for about 7 % of the total variation compared to 88 % for climate drivers.

The relative proportion of model predictor coefficients was used as weighting factors for multicriteria analysis to create a suitability map for the GDV in southern Portugal showing where the vegetation most likely relies on groundwater to cope with aridity. A validation of the resulting map was performed using independent data of the normalized difference water index (NDWI), a satellite-derived vegetation index. June, July and August of 2005 NDWI anomalies, for the years 1999–2009, were calculated to assess the response of active woody species in the region after an extreme drought. The results from the NDWI anomalies provided an overall good agreement with the suitability to host GDV. The model was considered to be reliable for predicting the distribution of the studied vegetation.

The methodology developed to map GDVs will allow for the prediction of the evolution of the distribution of GDV according to climate change and aid stakeholder decision-making concerning priority areas of water resource management.

1 Introduction

Mediterranean forests, woodlands and shrublands, mostly growing under restricted water availability, are one of the terrestrial biomes with a higher volume of groundwater being used by vegetation (Evaristo and McDonnell, 2017). Future predictions of decreased precipitation (Giorgi and Lionello, 2008; Nadezhdina et al., 2015), decreased runoff (Mourato et al., 2015) and aquifer recharge (Ertürk et al., 2014; Stigter et al., 2014) in the Mediterranean region threaten the sustainability of groundwater reservoirs and the corresponding dependent ecosystems. Therefore, sustainable management of groundwater resources and the groundwater-dependent ecosystems (GDEs) is of crucial importance.

A widely used classification of the GDEs was proposed by Eamus et al. (2006). This classification distinguishes three types: (1) aquifer and cave ecosystems, which include all subterranean waters, (2) ecosystems reliant on emerging groundwater (e.g., estuarine systems and wetlands; riverine systems), and (3) ecosystems reliant on resident groundwater (e.g., systems where plants remain physiologically active during extended drought periods without a visible water source). Mapping GDEs constitutes a first and fundamental step to their active management. Several approaches have been proposed, from local field surveys measuring plant transpiration of stable isotopes (Antunes et al., 2018) up to larger spatial scales involving remote-sensing techniques (e.g., normalized difference vegetation index – NDVI; Barron et al., 2014; Eamus et al., 2015; Howard and Merrifield, 2010), remote-sensing combined with ground-based observations (Lv et al., 2013), a geographic information system (GIS; Pérez Hoyos et al.,, 2016a) combining field surveys (Condesso de Melo et al., 2015) or even statistical approaches (Pérez Hoyos et al., 2016b).

Despite a wide-ranging body of literature reviewing the GDE's topics (Doody et al., 2017; Dresel et al., 2010; Münch and Conrad, 2007), most regional-scale studies do not include Mediterranean regions. Moreover, studies on ecosystems relying on resident groundwater frequently only focused on riparian environments (Lowry and Loheide, 2010; O'Grady et al., 2006), with few examples in Mediterranean areas (del Castillo et al., 2016; Fernandes, 2013; Hernández-Santana et al., 2008; Mendes et al., 2016). There is a clear knowledge gap on the identification of phreatophyte species reliant on resident groundwater and their associated vegetation (Robinson, 1958) in the Mediterranean region and the management actions that should be taken to decrease the adverse effects of climate change.

In the driest regions of the Mediterranean Basin, the persistent lack of water during the entire summer periods gave an adaptive advantage to the vegetation that could either avoid or escape drought by reaching deeper-stored water up to the point of entirely relying on groundwater (Chaves et al., 2003; Canadell et al., 1996; Miller et al., 2010). This drought-avoiding strategy is often associated with the development of a dimorphic root system in woody species (Dinis, 2014; David et al., 2013) or with hydraulic lift and/or hydraulic redistribution mechanisms (Orellana et al., 2012). Those mechanisms provide the ability to move water from deep soil layers, where water content is higher, to more shallow layers where water content is lower (Horton and Hart, 1998; Neumann and Cardon, 2012). Hydraulic lift and redistribution have been reported for several woody species of the Mediterranean Basin (David et al., 2007; Filella and Peñuelas, 2004) and noticeably for cork oak (Quercus suber L.; David et al., 2013; Kurz-Besson et al., 2006; Mendes et al., 2016). Mediterranean cork oak woodlands (montados) are agro-silvo-pastoral systems considered to be semi-natural ecosystems of the southwestern Mediterranean Basin (Joffre et al., 1999) that have already been referenced as a groundwater-dependent terrestrial ecosystem (Mendes et al., 2016). Montados must be continually maintained through human management by thinning and understory use through grazing, plowing and shrub clearing (Huntsinger and Bartolome, 1992), to maintain good productivity, biodiversity and ecosystems services (Bugalho et al., 2009). In the ecosystems of this geographical area, the dominant tree species are the cork oak (Quercus suber L.) and the Portuguese holm oak (Quercus ilex subs rotundifolia Lam. (Pinto-Correia et al., 2011). Additionally, stone pine (Pinus pinea L.) has become a commonly co-occurrent species in the last decades (Coelho and Campos, 2009). The use of groundwater has been frequently reported for both the Pinus (Antunes et al., 2018; Filella and Peñuelas, 2004; Grossiord et al., 2016; Peñuelas and Filella, 2003) and Quercus genre (Barbeta and Peñuelas, 2017; Kurz-Besson et al., 2006, 2014; Otieno et al., 2006). Furthermore, the contribution of groundwater to tree physiology has been shown to be of a greater magnitude for Quercus sp. compared with Pinus sp. (del Castillo et al., 2016; Evaristo and McDonnell, 2017).

Q. suber and Q. ilex have been associated with high resilience and adaptability to hydric and thermic stress and with recurrent droughts in the southern Mediterranean Basin (Barbero et al., 1992). In Italy and Portugal, during summer droughts, Q. ilex used a mixture of rainwater and groundwater and was able to take water from very dry soils (David et al., 2007; Valentini et al., 1992). An increasing contribution of groundwater in the summer has also been shown for this species (Barbeta et al., 2015). Similarly, Q. suber showed a seasonal shift in water sources, from shallow soil water in the spring to the beginning of the dry period followed by a progressive higher use of deeper water sources throughout the drought period (Otieno et al., 2006). In addition, the species roots are known to reach depths as deep as 13 m in southern Portugal (David et al., 2004). P. pinea has been recently included in the facultative phreatophyte species (Antunes et al., 2018). This species shows a very similar root system (Montero et al., 2004) compared to cork oak (David et al., 2013), with large sinker roots reaching 5 m depth (Canadell et al., 1996). Given the information available on water use strategies by the phreatophyte arboreous species of the cork oak woodlands, Q. ilex, Q. suber and P. pinea were considered to be proxies for arboreous vegetation that belongs to the GDE relying on resident groundwater (hereafter designed as groundwater-dependent vegetation – GDV).

GDV of the Mediterranean Basin is often neglected in research. Still, little is known about the GDV distribution, but research has already been done on the effects of climate change in specific species distribution, such as that for Q. suber, in the Mediterranean Basin (Duque-Lazo et al., 2018; Paulo et al., 2015). While the increase in atmospheric CO₂ and the rising temperature can boost tree growth (Barbeta and Peñuelas, 2017; Bussotti et al., 2013; Sardans and Peñuelas, 2004), water stress can have a counteracting effect on the growth of both Quercus ilex (López et al., 1997; Sabaté et al., 2002) and P. pinaster (Kurz-Besson et al., 2016). Therefore, it is of crucial importance to identify geographical areas where subsurface GDV is present and characterize the environmental conditions this vegetation type is thriving in. This would contribute to the understanding of how to manage these species under unfavorable future climatic conditions.

The aim of this study was to address the mentioned gaps by creating a suitability map of the arboreous phreatophyte species in southern Portugal, traducing their potential dependency on groundwater. We used an integrated multidisciplinary methodology combining a geospatial modeling approach based on the geographically weighted regression (GWR) and a GIS multicriteria analysis approach, both relying on forest inventory, edaphoclimatic conditions and topographic information. We expected this new integrated procedure to grant a more reliable estimation of the vegetation dependency on groundwater sources at the regional scale.

The mapping methodology was based on the occurrence of known subsurface phreatophyte species and well-known environmental conditions affecting water resource availability. Several environmental predictors were selected according to their expected impact on water use and storage and then used in GWR to model the density of Q. suber, Q. ilex and P. pinea occurrence in the Alentejo region (NUTSII) of southern Portugal. To our knowledge, very few applications of GWR have been used to model species distribution, and only recently its use has spread in ecological research (Hu et al., 2017; Li et al., 2016; Mazziotta et al., 2016). The coefficients obtained from the model equation for each predictor and expressed as the proportion of the total sum of absolute coefficients were used as weights to build the suitability map with GIS multifactor analysis, after reclassifying each relevant environmental driver. The resulting map was validated using the remote-sensed vegetation index, the normalized difference water index (NDWI).

Based on former knowledge gathered from field surveys conducted in the region (Antunes et al., 2018; Condesso de Melo et al., 2015; Kurz-Besson et al., 2006, 2014; Otieno et al., 2006; David et al., 2013; Pinto et al., 2013) on environmental conditions and the species ecophysiological needs, we hypothesized that (1) groundwater depth together with climatic conditions play one of the most important environmental roles in GDV's distribution and (2) groundwater depth between 1.5 and 15 m associated with xeric conditions should favor a higher density of GDV and thus a larger use of groundwater by the vegetation.

2 Material and methods

2.1 Study area

The administrative region of Alentejo (NUTSII; Fig. 1) covers an area of 31 604.9 km², between 37.22 and 39.39^∘ N in latitude and between 6.55 and 9.00^∘ W in longitude. This study area is characterized by a Mediterranean temperate mesothermic climate with hot and dry summers, defined as Csa in the Köppen classification (APA, 2016; ARH Alentejo, 2012a, b). It is characterized by a subhumid climate, which has recently quickly drifted to semi-arid conditions (Ministério da Agricultura do Mar do Ambiente e do Ordenamento do Território, 2013). A large proportion of the area (above 40 %) is covered by forestry systems (Autoridade Florestal Nacional and Ministério da Agricultura do Desenvolvimento Rural e das Pescas, 2010), providing a high economical value to the region and the country (de Sarmento and Dores, 2013).

Figure 1 Study area. On the left is the location of Alentejo in the Iberian Peninsula; on the right is the elevation characterization of the study area with the main river courses from Tagus, Sado and Guadiana basins (white lines). Names of the main rivers are indicated at to their location on the map.

2.2 Kernel density estimation of GDV

Presence datasets of Quercus suber, Quercus ilex and Pinus pinea of the last Portuguese forest inventory completed in 2010 (ICNF, 2013) were used to calculate Kernel density (commonly called a heat map) as a proxy for GDV suitability. The inventory registered the occurrence of each species on a 500 m mesh grid resolution, corresponding to a maximum occurrence of four counts per square kilometer. Only data points with one of the three proxy species selected as primary and secondary occupation were used. The resulting Kernel density was weighted according to tree cover percentage and was calculated using a quartic biweight distribution shape, a search radius of 10 km, and an output resolution of 0.018^∘, corresponding to a cell size of 1 km. This variable was computed using QGIS version 2.14.12 (QGIS Development Team, 2017).

2.3 Environmental variables

Species distribution is mostly affected by limiting factors controlling ecophysiological responses, disturbances and resources (Guisan and Thuiller, 2005). To characterize the study area in terms of GDV's suitability, environmental variables expected to affect GDV's density were selected according to their constraint on groundwater uptake and soil water storage. Within possible abiotic variables, landscape topography, geology, groundwater availability and regional climate were considered. The 12 selected variables for modeling purposes, retrieved from different data sources, are listed in Table 1. The software used in spatial analysis was ArcGIS^® software version 10.4.1 by Esri and R program software version 3.4.2 (R Development Core Team, 2016).

Table 1 Environmental variables for the characterization of the suitability of GDV in the study area.

2.3.1 Slope and soil characteristics

The NASA and METI ASTER GDEM product was retrieved from the online data pool, courtesy of the NASA Land Processes Distributed Active Archive Center (LP DAAC), USGS Earth Resources Observation and Science (EROS) Center, Sioux Falls, South Dakota (https://lpdaac.usgs.gov/data_access/data_pool, last access: 15 June 2017). The Spatial Analyst toolbox was used to calculate the slope from the digital elevation model. The slope was used as a proxy for the identification of shallow soil water interaction with vegetation.

The map of soil type was obtained from the Portuguese National Information System for the Environment – SNIAmb (© Agência Portuguesa do Ambiente, 2017, https://sniamb.apambiente.pt/?language=pt-pt, last access: August 2019) – and uniformized to the World Reference Base with the Harmonized World Soil Database v 1.2 (FAO et al., 2009). The vector map was converted to raster using the conversion toolbox. To reduce the analysis complexity involving the several soil types present in the map, soil types were regrouped into three classes according to their capacity to store or drain water (Table A1). The classification was based on the characteristics of each soil unit (available water storage capacity, drainage and topsoil texture) from the Harmonized World Soil Database v 1.2 (FAO et al., 2009). In the presence of dominant soil with low drainage capacity, a high clay fraction in the top soil and a high available water content, lower scores were given in association with decreased suitability for GDV by favoring non-GDV species. Otherwise, when soil characteristics suggested water storage at deeper soil depths and lower water content, drainage and sandy topsoil texture, higher scores were given.

Effective soil thickness (Table 1) was also considered for representing the maximum soil depth explored by the vegetation roots. It constrains the expansion and growth of the root system as well as the available amount of water that can be absorbed by roots.

2.3.2 Groundwater availability

Root access to water resources is one of the most limiting factors for GDV's growth and survival, especially during the dry season. The map of depth to the water table was interpolated from piezometric observations from the Portuguese National Information System on Water Resources (SNIRH) public database (http://snirh.apambiente.pt, last access: 31 March 2017) and the Study of Groundwater Resources of Alentejo (ERHSA; Chambel et al., 2007). Data points of large-diameter wells and piezometers were retrieved for the Alentejo region (Fig. 2) and sorted into undifferentiated, karst or porous geological types to model groundwater depth. In the studied area, piezometers are exclusively dedicated small-diameter boreholes for piezometric observations in areas with high abstraction volumes for public water supply. Large-diameter wells in this region usually have a low yield and are mainly devoted to private use and irrigation. Due to the large heterogeneity of geological media, groundwater depth was calculated separately for each sub-basin. A total of 3158 data points corresponding to large wells and piezometers were used, with uneven measurements between 1979 and 2017. For each piezometer an average depth was calculated from the available observations and used as a single value. In areas with an undifferentiated geological type, piezometric level and elevation were highly correlated (>0.9); thus a linear regression was applied to interpolate data. Ordinary kriging was preferred for the interpolation of karst and porous aquifers, combining large wells and piezometric data points. The ordinary kriging was calculated using a semi-variogram in which the sill, range and nugget were optimized to create the best fit of the model to the data. To build a surface layer of the depth to the water table, the interpolated surface of the groundwater level was subtracted from the digital elevation model. The Geostatistical Analyst toolbox was used for this task.

Figure 2 Large well and piezometer data points used for groundwater depth calculation. Squares represent piezometers data points, and triangles represent large well data points.

Drainage density is a measure of how well the basin is drained by stream channels. It is defined as the total length of channels per unit area. Drainage density was calculated for a 10 km grid size for the Alentejo region by the division of the 10 km square area (A) in square kilometers by the total stream length (L) in kilometers, as in Eq. (1):

$$\begin{array}{}\text{(1)}& D={\displaystyle \frac{L}{A}}.\end{array}$$

2.3.3 Regional climate

Temperature and precipitation datasets were obtained from the E-OBS (http://eca.knmi.nl/download/ensembles/ensembles.php, last access: 31 March 2017) public database (Haylock et al., 2008). The standardized precipitation evapotranspiration index (SPEI), aridity index (A_i) and ombrothermic indices were computed from long-term (1951–2010) monthly temperature and precipitation observations. The computation of potential evapotranspiration (PET) was performed according to Thornthwaite (1948) and was calculated using the SPEI package (Beguería and Vicente-Serrano, 2013) in the R program.

The SPEI, a multi-scalar drought index (Vicente-Serrano et al., 2010), was calculated over a 6-month interval to characterize drought severity in the area of study using the SPEI package (Beguería and Vicente-Serrano, 2013) for the R program. The SPEI is based on the normalization of the water balance calculated as the difference between cumulative precipitation and PET for a given period at monthly intervals. Normalized values of the SPEI typically range between −3 and 3. Drought events were considered to be severe when SPEI values were between −1.5 and −1.99 and extreme with values below −2 (Mckee et al., 1993). Severe and extreme SPEI predictors were computed as the number of months with severe or extreme drought, counted along the 60 years of the climate time series.

While the SPEI used in this study identifies geographical areas affected with more frequent extreme droughts, the aridity index distinguishes arid geographical areas prone to annual negative water balance (with low A_i value) from more mesic areas showing a positive annual water balance (with high A_i value). A_i gives information related to evapotranspiration processes and the rainfall deficit for potential vegetative growth. It was calculated following Eq. (2) according to Middleton et al. (1992), where PET is the average annual potential evapotranspiration and P is the average annual precipitation, both given in millimeters for the 60-year period of the climate time series. Drylands are defined by their degree of aridity in four classes: hyperarid (A_i<0.05), arid ($\mathrm{0.05}<A_{\mathrm{i}}<\mathrm{0.2}$), semi-arid ($\mathrm{0.2}<A_{\mathrm{i}}<\mathrm{0.5}$) and dry subhumid ($\mathrm{0.5}<A_{\mathrm{i}}<\mathrm{0.65}$; Middleton et al., 1992):

$$\begin{array}{}\text{(2)}& A_{\mathrm{i}}={\displaystyle \frac{P}{\mathrm{PET}}}.\end{array}$$

Ombrothermic indices were used to better characterize the bioclimatology of the study region (Rivas-Martínez et al., 2011) by evaluating soil water availability for plants during the driest months of the year. Four ombrothermic indices were calculated according to a specific section of the year stated in Table 1, and following Eq. (2), where P_p  s the positive annual precipitation (accumulated monthly precipitation when the average monthly mean temperature is higher than 0 ^∘C) and T_p is the positive annual temperature (total in tenths of degrees centigrade of the average monthly temperatures higher than 0^∘). The ombrothermic index presenting values below 2 for the analyzed months can be considered to be Mediterranean bioclimatically. For non-Mediterranean areas, there is no dry period in which, for at least 2 consecutive months, the precipitation is less than or equal to twice the temperature:

$$\begin{array}{}\text{(3)}& O={\displaystyle \frac{P_{\mathrm{p}}}{T_{\mathrm{p}}}}.\end{array}$$

2.4 Selection of model predictors

The full set of environmental variables was evaluated as potential predictors for the suitability of GDV (based on the Kernel density of the proxy species). A preliminary selection was carried out, first by computing Pearson's correlation coefficients between environmental variables and second by performing a principal components analysis (PCA) to detect multicollinearity. Covariates were discarded for modeling according to a sequential procedure. Whenever pairs of variables presented a correlation value above 0.4, the variable with the highest explained variance on the first axis of the PCA was selected. In addition, selected variables had to show the lowest possible correlation values between them. Variables showing low correlations and explaining a higher cumulative proportion of variability with the lowest number of the PCA axis were later selected as predictors for modeling. PCA was performed using the GeoDa software (Anselin et al., 2006), and Pearson's correlation coefficients were computed with the Spatial Analyst tool.

2.5 Model development

When fitting a linear regression model based on the selected variables, the normal distribution and stationarity of the response variable and residuals must be assured.

The Kernel density of the proxy GDV species, Q. suber, Q. ilex and P. pinea, showed a skewed normal distribution. Therefore, a square-root transformation of the data was applied on the response variable before model fitting. To be able to compare the resulting model coefficients and use them as weighting factors of the multicriteria analysis to build the suitability map, the predictor variables were normalized using the z-score function. This allows us to create standardized scores for each variable by subtracting the mean of all data points from each individual data point and then dividing those points by the standard deviation of all points so that the mean of each z predictor is zero and the deviation is 1.

Spatial autocorrelation and non-stationarity are common when using linear regression on spatial data. To overcome these issues, the GWR was used. This extension of the ordinary least-squares (OLS) linear regression considers the spatial non-stationarity in variable relationships and allows the use of spatially varying coefficients while minimizing spatial autocorrelation (Stewart Fotheringham et al., 1996). In this study, simple linear regression and GWR were both applied to the dataset, and their performances were compared. Models were fitted on a 5 % random subsample of the entire dataset (reaching a total of 6214 selected data points) due to computational restrictions and to a decrease in the spatial autocorrelation effect (Kühn, 2007). Bertrand et al. (2016) applied a similar methodology by using only 10 % of the available points, with a ∼10 km distance between points in the selected sub-sample. In spite of the subsampling, the minimum and maximum distance between two random data points were, respectively, 3.6 and 16.7 km, providing a good representation of local heterogeneity, as shown in Figs. 5 and 6. An additional analysis showing an excellent agreement between the two datasets is presented in Fig. A1.

Initially the model was constructed containing all selected predictors through the PCA and Pearson's correlation analysis. Afterwards, predictors were sequentially discarded to ascertain the model presenting lower second-order Akaike information criteria (AIC_c) and higher quasi-global R² values chosen to predict the suitability of GDV.

Adaptive Kernel bandwidths for the GWR model fitting were used due to the spatial irregularity of the random subsample. A local search radius was obtained by minimizing the cross-validation score (Bivand et al., 2008), thus minimizing the error of the local regressions. To analyze the performance of the GWR model alone, the local and global adjusted R² values were considered. To compare between the GWR model and the simple linear model, the distribution of the model residuals was used to identify clustered values as well as the AIC_c. The spatial autocorrelation of the models residuals was evaluated with the Moran's I test (Moran, 1950) calculated from the Spatial Statistics tool and also graphically. GWR model was fitted using the spgwr package from the R program (Bivand and Yu, 2017).

2.6 Suitability map building

To create the suitability map all predictor layers included in the GWR model were classified, similarly to Condesso de Melo et al. (2015) and Aksoy et al. (2017). The likelihood of an interaction between the vegetation and groundwater resources was scored from 1 to 3 for each predictor. Scores were assigned after bibliographic review and expert opinion. The higher the score, the higher the likelihood, with 1 corresponding to a weak likelihood and 3 indicating a very high likelihood.

Groundwater depth was divided in two classes according to the accessibility to shallow soil water above 1.5 m and the maximum rooting depth for Mediterranean woody species reaching 13 m, reported by Canadell et al. (1996). Throughout the paper water between 0 and 1.5 m depth was designated as shallow soil water, while water below 1.5 m depth was considered to be groundwater. The depth class between 0 and 1.5 m was based on the riparian vegetation in semi-arid Mediterranean areas, which is mainly composed of shrub communities (Salinas et al., 2000) and presents a mean rooting depth of 1.5 m (Silva and Rego, 2004). The most common tree-species rooting depth in riparian ecosystems is normally similar to the depth of fine sediment not reaching gravel substrates (Singer et al., 2012) and not reaching levels as deep as deep-rooted species. The minimum score was given to areas where groundwater depth was too shallow (below 1.5 m) to belong to emerging groundwater-dependent vegetation. Areas with a steep slope were considered to have superficial runoff and less recharge and influence tree density negatively (Costa et al., 2008). Those areas were treated as less suitable to GDV. Values of the ombrothermic index of the summer quarter and the immediately previous month (O₄) were split into three classes according to Jenks natural breaks, with higher suitability corresponding to higher aridity. The higher values of A_i corresponding to lower aridity had a score of 1 because a higher humid environment would decrease the necessity of the arboreous species to use deep water sources. Accordingly, an increase in aridity (lower values of A_i) has already been shown to increase tree decline (Waroux and Lambin, 2012), and thus lower A_i values corresponded to a score of 3, leaving the score 2 to intermediate values of A_i. Drainage density scoring was based on the drainage capability of the water through the hydrographical network of the river. A low drainage density (below 0.5) implies a high loss of water through runoff along the hydrographic network. This water lost for shallow soil horizons would be less available to the vegetation, thus favoring a higher use of water from deep groundwater reservoirs (Rodrigues et al., 2011).

A direct compilation of the predictor layers could have been performed for the multicriteria analysis. However, some predictors might have a stronger influence on GDV's distribution and density than others. Therefore, there was a need to define weighting factors for each layer of the final GIS multicriteria analysis. Yet, due to the intricate relations between all environmental predictors and their effects on the GDV, experts and stakeholders suggested very different scoring for a same layer. Instead the relative proportion of each predictor was used locally, according to the GWR model (Eq. 4), as weighting factors. The final GIS multicriteria analysis was performed using the Spatial Analyst tool by applying local model equations obtained for each of the 6214 coordinates of the Alentejo map (Eq. 4),

$$\begin{array}{}\text{(4)}& {\displaystyle }\begin{array}{rl}{S}_{\mathrm{GDV}}& =\mathrm{Intercept}+{\mathrm{coef}}_{p\mathrm{1}}\cdot \left[\mathrm{reclassified}\mathrm{value}{X}_{\mathrm{1}}\right]\\ & +{\mathrm{coef}}_{p\mathrm{2}}\cdot \left[\mathrm{reclassified}\mathrm{value}{X}_{\mathrm{2}}\right]\\ & +{\mathrm{coef}}_{p\mathrm{3}}\cdot \left[\mathrm{reclassified}\mathrm{value}{X}_{\mathrm{3}}\right]+\mathrm{\dots },\end{array}\end{array}$$

with S_GDV representing the suitability to groundwater-dependent vegetation, brackets representing the reclassified GIS X layer corresponding to the scoring, and coef_x indicating the relative proportion for the predictor x calculated as the ratio between the modulus of the local coefficient x and the sum of the modulus of all local coefficients.

According to this equation, lower values indicate a lower occurrence of groundwater use, representing a lower GDV suitability, while higher values correspond to a higher use of groundwater, representing higher GDV suitability. To allow for an easier interpretation, the data on suitability to GDV were subsequently classified based on their distribution value, according to Jenks natural breaks. This resulted in five suitability classes: “Very poor”, “Poor”, “Moderate”, “Good” and “Very Good”.

2.7 Map evaluation

Satellite-derived remote-sensing products have been widely used to follow the impact of drought on land cover and the vegetation dynamics (AghaKouchak et al., 2015). Vegetation indices offer excellent tools to assess and monitor plant changes and water stress (Asrar, 1989). The NDWI (Gao, 1996) is a satellite-derived index that aims to estimate fuel moisture content (Maki et al., 2004) and leaf water content at canopy level, being widely used for drought monitoring (Anderson et al., 2010; Gu et al., 2007; Ceccato et al., 2002a). This index was chosen to be more sensitive to canopy water content and a good proxy for water stress status in plants. Moreover, the NDWI has been shown to be best related to the greenness of cork oak woodland's canopy, expressed by the fraction of intercepted photosynthetically active radiation (Cerasoli et al., 2016).

In order to validate the GDV suitability map obtained in our study, we calculated anomalies of the NDWI (Gao, 1996) between an extremely dry year (2005) and the median value of the surrounding 10-year period (1999–2009). The NDWI is computed using the near-infrared (NIR) and the shortwave infrared (SWIR) reflectance, which makes it sensitive to changes in liquid water content and in vegetation canopies (Gao, 1996; Ceccato et al., 2002a, b).  The index computation (Eq. 5) was further adapted by Gond et al. (2004) to SPOT-VEGETATION instrument datasets, using NIR (0.84 µm) and mid-infrared (MIR) (1.64 µm) channels, as described by Hagolle et al. (2005):

$$\begin{array}{}\text{(5)}& \mathrm{NDWI}={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{NIR}}-{\mathit{\rho }}_{\mathrm{MIR}}}{{\mathit{\rho }}_{\mathrm{NIR}}+{\mathit{\rho }}_{\mathrm{MIR}}}}.\end{array}$$

Following Eq. (5), NDWI data were computed using B3 and MIR data acquired from the VEGETATION instrument aboard SPOT 4 and SPOT 5 satellites. Extraction and corrections procedures applied to optimize NDWI series are fully described in Gouveia et al. (2012).

The NDWI anomaly was computed as the difference between NDWI observed in June, July and August of 2005 and the median NDWI for the considered month for the period 1999 to 2009. June was selected to provide the best signal from a still fully active canopy of woody species while the herbaceous layer had usually already finished its annual cycle and dried out. The hydrological year of 2004–2005 was characterized by an extreme drought event over the Iberian Peninsula, where less than 40 % of the normal precipitation was registered in the southern area (Gouveia et al., 2009). Thus, in June 2005 the vegetation of the Alentejo region was already coping with an extreme long-term drought, which was well captured by the anomaly of the NDWI (negative values), as formerly shown by Gouveia et al. (2012).

2.8 Sensitivity analysis

Sensitivity analyses are conducted to identify model inputs that cause significant impact and/or uncertainty in the output. They can be used to identify key variables that should be the focus of attention to increase mode robustness in future research or to remove redundant inputs from the model equation because they do not have significant impact on the model output. Based on bootstrapping simulations (Tian et al., 2014), a sensitivity analysis was conducted on the GWR model by perturbing one input predictor at a time while keeping the rest of the equation unperturbed. To simulate perturbations, 10 000 values were randomly selected within the natural range of each input variable observed in the Alentejo region. Those random values were then used to run 10 000 simulations of the local equation of the GWR model for each of the 6214 coordinates of the geographical area. Local outputs corresponding to the predicted GDV density were then calculated for each perturbed input variable (A_i, O₄, W, D and s). The range of output values was calculated to reflect the sensibility of the model for the perturbed input variable. The overall sensibility of the model to all input variables was estimated as the absolute difference between the minimum output value and the sum of maximum output values of all predictors, thus representing the maximum possible output range observed after perturbing all predictors.

3 Results

3.1 Kernel density

Within the studied region of Portugal, the phreatophyte species Quercus suber, Quercus ilex and Pinus pinea were not distributed uniformly throughout the territory. Areas with higher Kernel density (or higher distribution likelihood) were mostly spread between the northern part of Alentejo region and the western part close to the coast, with values ranging between 900 and 1200 occurrences in a 10 km search radius (Fig. 3). Two clusters of high density also appeared below the Tagus river. The remaining study area presented mean density values, with very low densities in the area of the river Tagus and in the center–south.

Table 2 Effect of variable removal in the performance of GWR model linking the Kernel density of Quercus suber, Quercus ilex and Pinus pinea (S_GDV) to predictors of the aridity index (A_i), ombrothermic index of the summer quarter and the immediately previous month (O₄), slope (s), drainage density (D), groundwater depth (W), and soil type (S_t). The model with all predictors is highlighted in grey, and the final model used in this study is in bold.

Figure 3 Map of Kernel density weighted by cover percentage of Q. suber, Q. ilex and P. pinea. The scale unit represents the number of occurrences per 10 km search radius (∼314 km²).

3.2 Environmental conditions

The exploratory analysis of the variables performed through the PCA and the Pearson's correlation matrix confirmed the presence of multicollinearity. From the initial variables (Table 1), thickness (T), number of months with severe and extreme SPEIs (respectively, SPEI_s and SPEI_e), annual ombrothermic index (O), ombrothermic index of the hottest month of the summer quarter (O₁), and ombrothermic index of the summer quarter (O₃) were discarded, while the variable slope (s), drainage density (D), soil type (S_t), groundwater depth (W), A_i and O₄ were maintained for analysis (Fig. A2 and Table A2). A sequential removal of one predictor from the initial modeling including six variables was performed (Table 2), after which the model was reduced to five variables. Therefore, out of the initial 12 variables considered (Fig. 4) to explain the variation of the Kernel density of GDV in Alentejo, the following variables were endorsed: A_i, O₄, W, D and s.

Figure 4 Map of environmental layers used in model fitting. (a) Aridity index, (b) ombrothermic index of the summer quarter and the immediately previous month, (c) groundwater depth, (d) drainage density, and (e) slope.

In most parts of the Alentejo region, the slope was below 10 % (Fig. 4e) and coastal areas presented the lowest values and variability. The highest values of groundwater depth (Fig. 4c), reaching a maximum of 255 m, were found in the Atlantic margin of the study area, mainly in Tagus and Sado river basins. Several other small and confined areas in Alentejo also showed high values, corresponding to aquifers of porous or karst geological types. Most of the remaining study area showed groundwater depths ranging between 1.5 and 15 m. Figure 4a and b indicate the southeast of Alentejo as the driest area, given by minimum values of A_i (0.618) and much higher potential evapotranspiration than precipitation. Besides this, O₄ presented a maximum value (1.166) for this region (meaning that soil water availability was not compensated by the precipitation of the previous May–June–July–August months). This is also supported by the higher drainage density in the southeast, which indicates a lower prevalence of shallow soil water due to higher stream length by area.

Combining all variables, it was possible to distinguish two subregions with distinct conditions: the southeast of Alentejo and the Atlantic margin. The latter is mainly distinguished by its low-slope areas, shallower groundwater and more humid climatic conditions than the southeast of Alentejo.

3.3 Regression models

The best model to describe the GDV distribution was found through a sequential discard of each variable (Table 2) and corresponded to the model with a distinct lower AIC_c (18050.76) than the second-lowest AIC_c (27389.74) and showed an important increase in quasi-global R² values (from 0.926 for the second-best model to 0.992 for the best one). The best model fit was obtained with A_i, O₄, W, D and s. This final model was then applied to the GIS layers to map the suitability of GDV in Alentejo, according to Eq. (6):

$$\begin{array}{}\text{(6)}& {\displaystyle }\begin{array}{rl}{S}_{\mathrm{GDV}}& =\mathrm{Intercept}+A_{\mathrm{i}}{\mathrm{coef}}_p\cdot \left[\mathrm{reclassified}{A}_{\mathrm{i}}\mathrm{value}\right]\\ & +O_{\mathrm{4}}{\mathrm{coef}}_p\cdot \left[\mathrm{reclassified}{O}_{\mathrm{4}}\mathrm{value}\right]\\ & +W{\mathrm{coef}}_p\cdot \left[\mathrm{reclassified}W\mathrm{value}\right]\\ & +D{\mathrm{coef}}_p\cdot \left[\mathrm{reclassified}D\mathrm{value}\right]\\ & +s{\mathrm{coef}}_p\cdot \left[\mathrm{reclassified}s\mathrm{value}\right].\end{array}\end{array}$$

The local adjusted R² of the GWR model was highly variable throughout the study area, ranging from 0 to 0.99 (Fig. 5), however the local R² values below 0.5 corresponded to only 0.3 % of the data. The lower R² values were distributed throughout the Alentejo area, with no distinct pattern. The overall fit of the GWR model was high (Table 3). The adjusted regression coefficient indicated that 99 % of the variation in the data was explained by the GWR model, while only 2 % was explained by the simple linear model (Table 3). Accordingly, GWR had a substantially lower AIC_c when compared with the simple linear model, indicating a much better fit.

Figure 5 Spatial distribution of local R² from the fitting of the geographically weighted regression.

The spatial autocorrelation given by the Moran index (Griffith, 2009; Moran, 1950) retrieved from the geospatial distribution of residual values was significant for both the GWR and the linear models, indicating that observations are geospatially dependent on each other to a certain level. However, this dependence was substantially lower for the GWR model than for the linear model (z score of 50.24 and 147.56 respectively). In the GWR model (Fig. 6a) the positive and negative residual values were much more randomly scattered throughout the study region than in the linear model (Fig. 6b), highlighting a much better performance of the GWR, which minimized residual autocorrelation. Indeed, in the linear model (Fig. 6b), positive residuals were condensed in the right side of Tagus and Sado river basins, while negative values were mainly present on the left side of the Tagus river and in the center–south of Alentejo.

Figure 6 Spatial distribution of model residuals from the fitting of the geographically weighted regression (a) and simple linear model (b).

The spatial distribution of the coefficients of GWR predictors is presented in Fig. 7. They were later used for the computation of the GDV suitability score for each data point (Eq. 6). The coefficient variability was 3 times higher for the A_i when compared to O₄ (Fig. 8a), reaching 66 % and 22 %, respectively. For W, D and s, the coefficient variation was much lower, representing only about 6.2 %, 3.8 % and 1.2 % of the total variation observed in the coefficients, respectively. The remaining variables showed a median close to 0, and the O₄ was the second, with higher variability, followed by the W. The coefficient median values were, respectively, −3.40, 0.29, −0.015, −0.018 and 0.022 for A_i, O₄, W, D and s variables.

Table 3 Comparison of adjusted R² values and second-order Akaike information criterion (AIC_c) between the simple linear regression and the GWR model.

^* Quasi-global R².

Figure 7 Map of local model coefficients for each variable. (a) Aridity index, (b) ombrothermic index of the summer quarter and the immediately previous month, (c) groundwater depth, (d) drainage density, and (e) slope.

Figure 8 Box plot of GWR model coefficient values for each predictor (a), and box plot of the GWR model outputs, corresponding to GDV's density after each of the predictors was disturbed for the sensitivity analysis (b). A_i stands for aridity index, O₄ for the ombrothermic index of the hottest month of the summer quarter and the immediately previous month, W for the groundwater depth, D for the drainage density, and s for the slope. Error bars represent the 25th and 75th percentile, while crosses indicate the 95th percentile.

The distributions of negative coefficients were similar for A_i and the O₄ variables (Fig. 7a and b), with lower values in the southern coastal area and in the Tagus river watershed. The highest absolute values were mostly found for A_i in the southern area of the Alentejo region and in smaller patches in the northern region. In the central and eastern areas of Alentejo, a higher weight of the groundwater depth coefficient could be found (Fig. 7c), approximately matching a higher influence of slope (Fig. 7e). The groundwater depth seemed to have almost no influence on GDV density in the Tagus river watershed, being expressed by coefficients that were mostly zero around the riverbed (Fig. 7c). The coefficient distribution of D and O₄ shows some similarities, mostly in the center and southeast of Alentejo (Fig. 7d). Extreme values of O₄ coefficients were mostly concentrated in the eastern part of the Tagus watershed and in the southern coastal area included in the Sado watershed. Slope coefficient values showed the lowest amplitude throughout the study area (Fig. 7e), with prevailing high positive values gathered mainly in the center of the study area and in the Tagus river watershed (northwest of the study center).

3.4 GDV suitability map

The classification of the five endorsed environmental predictors is presented in Table 4, and their respective maps are presented in Fig. B1. The rivers Tagus and Sado had an overall large impact on GDV's suitability for each predictor, with the exception of W. This is due to a higher water availability reflected by the values of O₄ and D and lower slopes due to the alluvial plains of the Tagus river (Fig. B1b, d and e). Moreover, those regions presented higher humidity conditions (through analysis of the A_i in Fig. B1a) and groundwater depths outside the optimum range (Fig. B1c), therefore being less suitable for GDV. Optimal conditions for groundwater access were mainly gathered in the interior of the study region (Fig. B1c), with the exception of some confined aquifers in the northeast and southeast of the study region. Favorable slopes for GDV were mostly highlighted in the Tagus river basin area, where a good likelihood of interaction between GDV and groundwater could be identified (Fig. B1e).

Table 4 Classification scores for each predictor. A score of 3 was given to highly suitable areas, and a score of 1 was given to less suitable areas for GDV.

The final map illustrating the suitability to GDV is shown in Fig. 9. The largest classified area (8787 km²) presented a very poor suitability to GDV, corresponding to approximately a quarter of the total study area (29 %). This percentage was followed closely by the moderate suitability to GDV, which occupied 26 % (8000 km²). Overall, the two less suitable classes (very poor and poor) represented 47 % of the study area, whilst the two best ones and the moderate class (very good, good and moderate) represented 53 %. Consequently, most of the study area showed moderate to high suitability to GDV. The very good and good suitability classes cover an arch from the most southern and northeastern area of the Alentejo region, passing through the Sado and southern Guadiana river basins and close to the coastal line at 38^∘ N. Most of the center of the study area showed moderate to very good suitability to GDV, while the areas corresponding to the alluvial deposits of the Tagus river showed poor to very poor suitability.

Figure 9 Suitability map for groundwater-dependent vegetation.

The suitability to GDV in the Alentejo region was mainly driven by A_i, given that the highest coefficient variability was associated with the A_i predictor in the GWR model equation. Consequently a similar distribution pattern can be observed between the suitability map and the aridity index predictor (Figs. 4a and 9). Areas with good or very good suitability mostly matched areas of A_i with a score of 3, corresponding to aridity index values above 0.75 (Fig. B1a). On the other hand, the lowest suitability classes showed a good agreement with the lowest scores given to W (Fig. B1c), mostly in the coastal area and in the Tagus river basin.

3.5 Map evaluation

To evaluate the suitability map developed in the present study, the results were compared with the NDWI anomaly considering the month of June of the dry year of 2005 in the Alentejo area (Fig. 10). Both maps (Figs. 9 and 10) showed similar patterns, with higher presence of GDV satisfactorily matching areas with the lowest NDWI anomaly. From June to September in an extremely dry year, non-GDV plants can be expected to experience severe drought stress as in any regular summer period. Thus, those plants should show almost zero anomaly. In contrast, GDV plants coping well with usual summer drought can be expected to suffer an unusual stress under an extremely dry year, even having access to groundwater (Kurz-Besson et al., 2006, 2014; Otieno et al., 2006; David et al., 2013), with a negative impact of groundwater drawdown (Antunes et al., 2018). Therefore, GDV plants should show negative NDWI anomalies.

Figure 10 Spatial patterns of NDWI anomaly values considering the months of June, July and August of the extremely dry year of 2005, in reference to the same months of the period 1999–2009, in the Alentejo region. Dark brown colors (corresponding to extreme negative NDWI anomaly values) indicate the vegetation that experienced the highest loss of water in leaves in summer 2005 compared to the reference period 1999–2009, while light brown colors show NDWI anomaly values very close to the usual vegetation moisture condition of the considered month.

The NDWI anomaly was mostly negative over the Alentejo territory, indicating a lower leaf water content in June and July 2005 than usual. The loss of water attributed to the extreme drought mostly matched geographical areas with the highest GDV suitability (Fig. 9). Water loss was less pronounced in the central area of the Alentejo region between the Guadiana and Sado river basins, where the vegetation is less dense (Fig. 3). Areas with NDWI anomaly values of zero (indicating no NDWI change) were mostly distributed on the coastal area of the Atlantic ocean or close to riverbeds, namely in the Tagus and Sado floodplains, matching areas of very poor suitability for GDV in Fig. 9.

Despite an overall good agreement, the adequation between the density, suitability and NDWI maps was not perfect. Indeed, some patches showing a high vegetation occurrence and density and large NDWI anomalies also matched an area of very poor suitability for GDV.

3.6 Sensitivity analysis

The sensitivity of the model in response to the perturbation of each input variable (A_i, O₄, W, D and s) is presented in Fig. 11a to e. The overall sensitivity of the model is further presented in Fig. 11f. For any input variable, the model sensitivity (Fig. 11a to e) was higher where absolute values of local coefficients were also higher (Fig. 7a to e). The maximum impact on GDV's density, corresponding to the maximum output range observed after perturbation (Fig. 8b), was observed when perturbing the A_i, accounting for 66 % of the total variability. The second-highest impact was observed after perturbing the O₄, corresponding to 22 %. The variability in the model outputs observed after perturbing the remaining variables W, D and s accounted for 7 %, 4 % and 1 % of the total accumulated variability, respectively (Fig. 8b). The highest variability in the GWR model output was mostly observed in the central part of the southern half of the Alentejo region as well as close to the main channels of the Guadiana and Tagus rivers (Fig. 11f). Furthermore, areas with higher model sensitivity (Fig. 11f) significantly matched higher model performance expressed by R² (Fig. 5), assessed with a Kruskal–Wallis test (p<0.0001).

Figure 11 Sensitivity analysis performed on the GWR model by perturbing one of the predictors while keeping the rest of the model equation constant. Graphics present the output range of GDV's density when the aridity index (a), the ombrothermic index (b), the groundwater depth (c), the drainage density (d) or the slope variable (e) was perturbed. The maximum possible range combining all predictors (f). The 95th percentile was used for the maximum value of the color bar for a better statistical representation of the spatial variability.

4 Discussion

4.1 Modeling approach

The geographically weighted regression model has been used before in ecological studies (Li et al., 2016; Mazziotta et al., 2016) but never for the mapping of GDV, to our knowledge. This approach considerably improved the goodness of fit when compared with the linear model, with a coefficient of regression (R²) increasing from 0.02 to 0.99 at the global level and an obvious reduction of residual clustering. Despite those improvements, it has not been possible to completely eliminate the residual autocorrelation after fitting the GWR model.

Kernel density for the study area provided a strong indication of the presence and abundance of the tree species considered to be a GDV proxy for modeling. The Mediterranean cork woodlands dominate about 76 % of the Alentejo region (while only 7 % is covered by stone pine). In those systems, tree density is known to be a trade-off between climate drivers (Joffre et al., 1999; Gouveia and Freitas, 2008) and the need for space for pasture or cereal cultivation in the understory (Acácio and Holmgreen, 2014). In our study, the anthropologic management of agroforestry systems in the Alentejo region has not been taken into account. According to a recent study of Cabon et al. (2018) where thinning played an important role in Q. ilex density in a Mediterranean climate site, anthropologic management could, at least partially, explain the non-randomness of the residual distribution after GWR model fitting as well as the mismatches between the GDV and the NDWI evaluation maps.

Another explanation of the reminiscent autocorrelation after GWR fitting could be the lack of groundwater-dependent species in the model. For example, Pinus pinaster Aiton was excluded due to its more humid distribution in Portugal and due to conflicting conclusions taken from previous studies to pinpoint the species as a potential groundwater user (Bourke, 2004; Kurz-Besson et al., 2016). In addition, olive trees were also excluded, although the use of groundwater by an olive orchard has been recently proved (Ferreira et al., 2018), however with a weak contribution of groundwater to the daily root flow and thus with no significant impact of groundwater on the species physiological conditions.

Methods previously used by Doody et al. (2017) and Condesso de Melo et al. (2015) to map specific vegetation relied solely on expert opinion, e.g., the Delphi panel, to define weighting factors of environmental information for GIS multicriteria analysis. In our study, the GWR modeling approach was used to assess weighting factors for each environmental predictor in the study area to build a suitability map for the GDV in southern Portugal. This allowed an empirical determination of the local relevance of each environmental predictor in GDV distribution, thus avoiding the inevitable subjectivity of Delphi panels.

Also, by combining the GWR and GIS approaches, we believe that the final suitability map provides a more reliable indication of the higher likelihood for groundwater dependency and a safer appraisal of the relative contribution of groundwater by facultative deep-rooted phreatophytes species in the Alentejo region.

Modeling of the entire study area at a regional level did not provide satisfactory results. Therefore, we developed a general model varying locally according to local predictor coefficients. The local influence of each predictor was highly variable throughout the study area, especially for climatic predictors reflecting water availability and stress conditions. The application of the GWR model not only allowed for a localized approach, by decreasing the residual error and autocorrelation over the entire studied region, but also provided insights on how GDV's density can be explained by the main environmental drivers locally.

The GWR model appeared to be highly sensitive to coefficient fitting corresponding to a good model fit, as expected in a spatially varying model. As such, high coefficients are highly reliable in the GWR model in our study. Yet, the high spatial variability of local coefficients might reflect a weak physical meaning of the GWR model that challenges its direct application in other regions, even under similar climate conditions. Predictor coefficients showed a similar behavior in the spatial distribution of the coefficients. This was noticeable for the aridity index and the groundwater depth in the Tagus and Sado river basins. Groundwater depth had no influence on GDV's density in these areas, and similarly, the coefficient of the aridity index showed a negative effect of increased humidity on GDV's density. In addition, a cluster of low drainage density values matched these areas. Due to the lower variability and impact of the drainage density and slope on the GDV's density, these variables might not impact significantly this vegetation density in future climatic scenarios.

4.2 Suitability to groundwater-dependent vegetation

According to our results, more than half of the study area appeared to be suitable for GDV. However, one-fourth of the studied area showed lower suitability to GDV. The lower suitability to this vegetation in the more northern and western part of the studied area included the coastal area and the Tagus river basin. Those are the moist humid areas of the study area, where GDV is unlikely to rely on groundwater during the drought season because rainfall water stored in shallow soil horizons is mostly available.

The proxy species (cork oak, holm oak and stone pine) can grow perfectly under subhumid Mediterranean climate conditions without relying as much on groundwater to survive as in more xeric semi-arid areas (Abad Vinas et al., 2016). As facultative phreatophyte species, their presence and abundance is only an indication of a possible use of groundwater. The study provided by Pinto et al. (2013) showed that cork oak, for example, can perfectly thrive where very shallow groundwater is available while suffering drought stress where groundwater source is lower but still extracted by trees. Also, former studies have shown that in the extremely dry year of 2005, cork oak experienced a severe drought stress, close to the cavitation threshold, although its main water source was groundwater (David et al., 2013; Kurz-Besson et al., 2006, 2014). These findings can explain that part of the maximum density (Fig. 4) matches the area of very poor suitability for GDV (Fig. 9). Elsewhere, the better agreement between the two maps reflects the dominance of the aridity index on the vegetation's occurrence. Groundwater depth appeared to have a lower influence on GDV density than climate drivers, as reflected by the relative low magnitude of the W coefficient and outputs of our model outcomes. This surprisingly disagrees with our initial hypothesis because groundwater represents a notable proportion of the transpired water of deep-rooting phreatophytes, reaching up to 86 % of absorbed water during drought periods and representing about 30.5 % of the annual water absorbed by trees (David et al., 2013; Kurz-Besson et al., 2014). Nonetheless, this disagreement should be regarded cautiously due to the poor quality of piezometric data used and the complexity required for modeling the water table depths. Besides this, the linear relationship between water depth and topography applied to areas of undifferentiated geological type can be weakened by a complex nonlinear interaction between topography, aridity and subsurface conductivity (Condon and Maxell, 2015). Moreover, the high variability in geological media, topography and vegetation cover at the regional scale did not allow us to account for small changes in groundwater depth (<15 m), which has a huge impact on GDV suitability (Canadell et al., 1996; Stone and Kalisz, 1991). Indeed, a high spatial resolution of the hydrological database is essential in rigorously characterizing the spatial dynamics of groundwater depth between hydrographic basins (Lorenzo-Lacruz et al., 2017). Unfortunately, such resolution was not available for our study area.

The aridity and ombrothermic indices were the most important predictors of GDV density in the Alentejo region, according to our model outcomes. Our results agree with previous findings linking tree cover density and rooting depth to climate drivers such as aridity at a global scale (Zomer et al., 2009; Schenk and Jackson, 2002) and specifically for the Mediterranean oak woodland (Gouveia and Freitas, 2008; Joffre et al., 1999). Through previous studies showing the similarities in vegetation strategies to cope with water scarcity in the Mediterranean Basin (Vicente-Serrano et al., 2013) or the relationship between rooting depth and water table depth increasing with aridity at a global scale (Fan et al., 2017), we can admit that the most relevant climate drivers identified here are similarly important for mapping GDV in other semi-arid regions. In this study, the most important environmental variables that define GDV's density in a semi-arid region were identified, helping to fill the gap of knowledge for modeling this type of vegetation. However, the coefficients to be applied when modeling each variable need to be calculated locally due to their high spatial variability.

Temporal piezometric data would further help discriminate areas of optimal suitability to GDV, either during the wet or the dry seasons, because the seasonal trends in groundwater depth are essential under Mediterranean conditions. Investigative efforts should be invested to fill the gap either by improving the Portuguese piezometric monitoring network or by assimilating observations with remote-sensing products focused on soil moisture or groundwater monitoring. This has already been performed on a large regional scale such as that for GRACE satellite surveys, based on changes of Earth's gravitational field. So far, these technologies are not applicable to Portugal's scale, since the coarse spatial resolution of GRACE data only allows the monitoring of large reservoirs (Xiao et al., 2015).

4.3 Validation of the results

The understory of woodlands and the herbaceous layer of grasslands areas in southern Portugal usually end their annual life cycles in June (Paço et al., 2009), while the canopy of woody species is still fully active with maximum transpiration rates and photosynthetic activities (Kurz-Besson et al., 2014; David et al., 2007; Awada et al., 2003). This is an ideal period of the year to spot differential response of the canopy of woody species to extreme droughts events using satellite-derived vegetation indices (Gouveia et al., 2012).

The spatial patterns of the NDWI anomaly in June 2005 seem to indicate that the woody canopy showed a strong loss of canopy water in the areas where tree density and GDV suitability were higher (Figs. 3, 9 and 10). This occurred although trees minimized the loss of water in leaves with a strong stomatal limitation in response to drought (Kurz-Besson et al., 2014; Grant et al., 2010). In the most arid area of the region were holm oak is dominant but tree density is much lower, the NDWI anomaly was generally less negative, thus showing a lower water stress or higher canopy water content. Holm oak (Quercus ilex spp. rotundifolia) is well known as the most resilient species to dry and hot conditions in Portugal, due to its capacity to use groundwater, associated with a higher water use efficiency (David et al., 2007). Furthermore, the dynamics of the NDWI anomaly over the summer period (Figs. 10a–c) pointed out that the lower water stress status on the map is progressively spreading from the most arid areas to the milder ones from June to August 2005 despite the intensification of drought conditions. This endorses the idea that trees manage to cope with drought by relying on deeper water sources in response to drought, replenishing leaf water content despite the progression and intensification of drought conditions. Former studies support this statement by showing that groundwater uptake and hydraulic lift progressively took place after the onset of drought by promoting the formation of new roots reaching deeper soil layers and water sources, typically from July onwards, for cork oak in the Alentejo region (Kurz-Besson et al., 2006, 2014). Root elongation following a declining water table has also been reported in a review on the effect of groundwater fluctuations on phreatophyte vegetation (Naumburg et al., 2005).

Our results and the dynamics of NDWI over summer 2005 tend to corroborate the studies of Schenk and Jackson (2002) and Fan et al. (2017) by suggesting a larger and longer dependency of GDV on groundwater with higher aridity. Further investigation needs to be carried out across aridity gradients in Portugal and the Iberian Peninsula to fully validate this statement, though.

Overall, the map of suitability to GDV showed good agreement with the NDWI validation maps. The main areas showing good GDV suitability and highest NDWI anomalies are mostly matching in both maps. The good agreement between our GDV suitability maps and NDWI dynamic maps opens up the possibility of applying and extending the methodology to larger geographical areas such as the Iberian Peninsula and to the simulation of the impact of climate changes on the distribution of groundwater-dependent species in the Mediterranean Basin.

Simulations of future climate conditions based on RCP4.5 and RCP8.5 emission scenarios (Soares et al., 2015, 2017) predict a significant decrease in precipitation for the Guadiana basin and overall decrease for the southern region of Portugal by 2100. Agroforestry systems relying on groundwater resources, such as cork oak woodlands, may show a decrease in productivity and ecosystem services or even face sustainability failure. Many studies carried out on oak woodlands in Italy and Spain identified drought as the main driving factor of tree dieback and as the main climate warning threatening oak stand sustainability in the Mediterranean Basin (Gentilesca et al., 2017). An increase in aridity and drought frequency for the Mediterranean (Spinoni et al., 2017) will most probably induce a geographical shift of GDV vegetation toward milder and wetter climates (Lloret et al., 2004; Gonzalez, 2001).

4.4 Key limitations

The GWR modeling approach used to estimate weighting factors is mostly stochastic. Consequently, the large spatial variability and symmetrical fluctuations around zero (Fig. 8b) denote a weak physical meaning of the estimated coefficients, at least at the resolution chosen for the study. Also, the local nature of the regression coefficients makes the model difficult to directly apply in other regions, even with similar climate conditions, unless the methodology is properly fitted to local conditions and predictors.

With the methodology applied in this study, weighting factors can be easily evaluated solely from local and regional observations of the studied area. Nonetheless, the computation of model coefficients or expert opinion to assess weighting factors requires recurrent amendments, associated with updated environmental data, species distribution and revised expert knowledge (Doody et al., 2017).

The evolution of groundwater depth in response to climate change is difficult to model on a large scale based on piezometric observations because it requires excellent knowledge of the components and dynamics of water catchments. Therefore, a reliable estimation of the impact of climate change on GDV suitability in southern Portugal could only be performed in small-scale studies. However, the GWR model appeared to be much more sensitive to climate drivers than the other predictors, given that 88 % of the model output variability was covered by climate indices A_i and O₄. Nevertheless, changes in climate conditions only represent part of the water resources shortage issue in the future. Global-scale changes in human populations and economic progress also control water demand and supply, especially in arid and semi-arid regions (Vörösmarty et al., 2000). A decrease in useful water resources for human supply can induce an even higher pressure on groundwater resources (Döll, 2009), aggravating the water table drawdown caused by climate change (Ertürk et al., 2014). Therefore, additional updates of the model should include human consumption of groundwater resources, identifying areas of higher population density or intensive farming. Future model updates should also account for the interaction of deep-rooted species with the surrounding understory species. In particular, shrubs surviving the drought period, which can benefit from the redistribution of groundwater by deep-rooted species (Dawson, 1993; Zou et al., 2005), should be accounted for.

5 Conclusions

Our results show a highly dominant contribution of water scarcity of the last 30 years (Aridity and Ombrothermic indices) to the density and suitability of deep-rooted groundwater-dependent species in southern Portugal. Therefore, in geographical regions of the world with similar semi-arid climate conditions (Csa according to Köppen–Geigen classification, Peel et al., 2007) and similar physiological responses of the groundwater-dependent vegetation (Vicente-Serrano et al., 2013), the use of the aridity and ombrothermic indices could be used as the first approximation to model and map the deep-rooted phreatophyte species and the evolution of their distribution in response to climate changes. The contribution of groundwater depth was lower than initially expected; however, it might be underestimated due to the poor quality of the piezometric network, especially in the central area of the studied region.

The current pressure applied by human consumption of water sources has reinforced the concern on the future of economic activities dependent on groundwater resources. To address this issue, several countries have developed national strategies for the adaptation of water sources for agriculture and forests against climate change, including Portugal (FAO, 2007). In addition, local drought management as a long-term adaptation strategy has been one of the proposals by Iglesias et al. (2007) to reduce the climate change impact on groundwater resources in the Mediterranean. The preservation of Mediterranean agroforestry systems, such as cork oak woodlands and the recently associated P. pinea species, is of great importance due to their high socioeconomic value and their supply of valuable ecosystem services (Bugalho et al., 2011). Management policies in the long-term should account for groundwater resources monitoring, accompanied by defensive measures to ensure that agroforestry system sustainability and economical income from these Mediterranean ecosystems are not greatly and irreversibly threatened.

Our present study, and novel methodology, provides an important tool to help in delineating priority areas of action for species and groundwater management, at the regional level, to avoid the decline of productivity and cover density of the agroforestry systems of southern Portugal. This is important for guaranteeing the sustainability of the economical income for stakeholders linked to the agroforestry sector in that area. Furthermore, mapping vulnerable areas at a small scale (e.g., by hydrological basin), where reliable groundwater depth information is available, should provide further insights for stakeholders to promote local actions to mitigate climate change's impact on GDV.

Based on the methodology applied in this work, future predictions on GDV suitability, according to the RCP4.5 and RCP8.5 emission scenarios, will be shortly introduced, providing guidelines for future management of these ecosystems in the allocation of water resources.

Appendix A: Data exploration

Table A1 Classification scores for the soil type predictor.

Figure A1 Box plot of the main predictors for the final geographically weighted regression model fitting (a, b) and the response variable (c, d) for the total data (a, c) and for the 5 % subsample (b, d).

Figure A2 Correlation plot between all environmental variables expected to affect the presence of the groundwater-dependent vegetation. O₁, O₃ and O₄ are ombrothermic indices of, respectively, the hottest month of the summer quarter, the summer quarter, and the summer quarter and the immediately previous month; O is the annual ombrothermic index, SPEI_e and SPEI_s are, respectively, the number of months with extreme and severe standardized precipitation evapotranspiration index; A_i is aridity index; W is groundwater depth; D is the drainage density; T is thickness; and S_t refers to soil type.

Table A2 Correlations between predictor variables and principal component axis. The most important predictors for each axis (when squared correlation is above 0.3) are showed in bold. The cumulative proportion of variance explained by each principal component axis is shown at the bottom of the table. s is slope; A_i is aridity index; O, O₁, O₃ and O₄ are ombrothermic indices of, respectively, the year, the hottest month of the summer quarter, the summer quarter, and the summer quarter and the immediately previous month; SPEI_s and SPEI_e are, respectively, the number of months with severe and extreme standardized precipitation evapotranspiration index; W is groundwater depth; D is the drainage density; S_t refers to soil type; and T is thickness.

Appendix B: Predictor reclassification

Figure B1 Predictor maps after score reclassification. (a) Aridity index, (b) ombrothermic index of the summer quarter and the immediately previous month, (c) groundwater depth, (d) drainage density, and (e) slope.