Elucidating Past Climate, Land Cover and Groundwater Trends in the Breede Gouritz Water Management Area, as Well as Anticipated Outcomes for the Near Future

1. Introduction

Climate and land cover change (CLCC) have been influencing our hydrologic cycle since before its effects were discovered. Comprehensive studies on climate change (CC) in the context of water resources have been done, including groundwater resources, albeit to a lesser extent than surface waters. Land use and land cover (LULC) change is a relatively new area of research in the context of surface and groundwaters (Nkhonjera & Dinka, 2018; Hegerl et al., 2019) although its effects on water resources have been proven (Kundu et al., 2017; Wang et al., 2018). CC is influencing the hydrological cycle and, therefore, our water resources in the following ways: higher average worldwide temperatures, hot areas becoming even hotter and vice versa, wet areas becoming even wetter and vice versa, melting of glacial ice, sea-level rise and saltwater encroachment, changes in evapotranspiration rates, runoff rates, recharge and discharge rates, more intense storms, prolonged droughts and so forth (Jayakumar & Lee, 2017; Smyth et al., 2017). Ultimately, these perturbations influence our available freshwater supply. Dwindling freshwater supplies are under an even more significant threat due to the rising population and subsequent demand for water supplies (McCusker & Ramudzuli, 2007; Mogelgaard, 2012).

CC impact studies on water resources have been explored by many authors. For some time, global circulation models (GCMs) and regional circulation models (RCMs) have been downscaled and applied to more local scenarios (Kumar, 2012; Trzaska & Schnarr, 2014). However, it has been found that these results are not as accurate as when one uses field data such as historical climate records with an analysis on a regional-to-local scale. The climate is a regional phenomenon with local variations (McGregor, 2018; Wang et al., 2018). Historical data can aid in understanding precipitation, temperature and evapotranspiration trends, which can be applied to catchment scale scenarios. The latter is done to predict changes in streamflow, runoff, discharge, baseflow, recharge and ultimately water availability in future.

LULC impact studies have also increased (Correia et al., 2023). It is debated whether CC or LULC change significantly impacts our water resources (Andaryani et al., 2019; Olivares et al., 2019). Conclusions need to be drawn on a local scale. LULC change in the form of urbanization influences surface albedo, influencing evapotranspiration. Furthermore, urbanization affects runoff rates, consequently influencing recharge into groundwater stores. Deforestation and afforestation also influence water resources (Olivares et al., 2019; Berhail, 2019). Agricultural practices, whether croplands or pasturelands and the type of crops all have an impact on the circulation of water through the system and, ultimately, the freshwater supply (Fohrer et al., 2001; Hu et al., 2019; Touhidul-Mustafa et al., 2019). These variations are site-specific and must be accounted for appropriately in impact studies. Many LULC tools have been developed over the years, such as CLUE, DynaCLUE, and SWAT1 (Albhaisi et al., 2013; Shrestha et al., 2018; Adhikari et al., 2020). Moreover, calculations can be applied to LULC maps to determine the magnitude of the change and make future predictions.

The Breede Gouritz Water Management Area (BGWMA) is essential in South Africa as it is an area with vast agricultural activity, numerous export trades of agricultural stock, a robust tourism industry, property development and a supplier of fresh water to the City of Cape Town, which is strategically and economically important to South Africa. The City of Cape Town almost ran out of water entirely during a drought in South Africa between 2015 and 2017. Many research endeavours have been undertaken in the catchment, including a status quo report on all available surface and groundwater resources. Whilst it is critical to be proactive concerning the available groundwater resources and plan for the needs of future generations, no study has determined the future of groundwater resources in the area under current CLCC projections.

This study, therefore, presents a new analysis with advanced knowledge of CLCC trends over the past few decades. The first objective is to determine any underlying trends over the years via a change detection analysis. Secondly, statistical climate and groundwater depth projections will be made to determine the magnitude of anticipated trends in the near future. Lastly, a relationship between the climate and climate-related variables will be established by statistical means, and empirical relationships will be established between the former variables and land cover change.

2. Study Area Description

The selected study area is the Breede-Gouritz Water Management Area (BGWMA), governed by the Breede-Gouritz Catchment Management Agency (BGCMA). The boundaries are the Indian Ocean to the South, the Berg-Olifants Water Management Area (WMA) to the West, the Orange WMA to the North and the Mzimvubu-Tsitsikama WMA to the east (Figure 1). This WMA is mainly in the Western Cape, with overlaps in the Eastern and Northern Cape. It covers an area of about 72,000 km², and has multiple land uses (Figure 2). The protection, conservation, management and control of water resources in the catchment are mandated by the BGCMA. The BGWMA provides an exciting platform to assess CC and LC impacts.

Generally, river water quality and conditions are good, but there is a call to use water more sustainably. There are 717 dams recorded in the water management area, ranging from very small to large. However, only 27 dams hold at least 5 million m³ of water. Managing groundwater resources is challenging because there is not a lot of data available. There are boreholes in the area from which data can be extracted (for example, water level, quality and geology logs), which will help to assess each aquifer. These boreholes are, however, mostly privately owned, from which data are rarely reported to the Department of Water and Sanitation (DWS). Another reason groundwater management is problematic is the underlying geology’s complexity and resulting water flow paths. The surface geology is inconclusive regarding the underlying geology and aquifer type and needs further investigation.

Precipitation enormously affects groundwater recharge and groundwater storage. It is imperative to determine water levels to see if abstraction rates are responsible and sustainable (Department of Water and Sanitation, 2017). More local and site-specific information is necessary to make proper conclusions about the groundwater situation in the BGWMA.

3. Materials and Methods

3.1. Data Sources

3.1.1. Climate Data

Climate data, temperature (T) and precipitation (P), in particular, were obtained from two sources: the Agricultural Research Council (ARC) and the South African Weather Service (SAWS). T was measured in degrees Celsius (°C), and the monthly average of the maximum and minimum daily average was used. P was evaluated by the sum of the monthly precipitation in mm. The ARC provided information from 11 stations scattered across the BGWMA dating back to 1924. The SAWS provided information from 25 stations scattered across the BGWMA dating back to 1990. Both sources reported maximum and minimum T and P data observed on a daily time step. Linear interpolation was used to calculate missing T values. In contrast, in the case of P, cross-correlations with neighbouring stations between the ARC and SAWS datasets were calculated, whereafter, the regression equation was applied between the two datasets to fill in missing values. Only one station in each selected subregion was chosen for further analysis. Figure 3 displays the geographic distribution of the selected ARC and SAWS stations.

3.1.2. Actual Evapotranspiration

Actual evapotranspiration (ET) maps that were recorded on a monthly timescale (in mm) have been obtained from the WaPOR site. These maps date back from January 2009 to the present. The data are in raster format with a resolution of 250 m. Monthly totals were obtained by estimating ET in mm/day, multiplying it by the number of days in a dekad (ten days), and then using each month’s average. Relevant data were extracted per Köppen-Geiger climate subregion, 8 in total. The extracted pixel values were multiplied by a factor of 0.1 and used as a variable in a time-series analysis. No data pixel values were omitted. Moreover, individual shapefiles were created for the four selected subregions (BWh, BWk, Csb and Csa) to compare each LC class’s area change per subregion between 1990, 2013/14 and 2020.

3.1.3. Groundwater Data

Groundwater data from the Breede and Gouritz Catchments have been obtained from the National Groundwater Archives (NGA) and the Hydstra database from the DWS. All borehole recordings recorded static water levels. In some cases, an autographic recorder was used to record data in many boreholes, meaning multiple daily measurements could be made. However, for most boreholes, monthly observations were recorded. A dip meter was used for older recordings. After filtering the NGA data, only 35 boreholes from the Breede Catchment and 130 from the Gouritz Catchment had data recorded for at least ten years. From this dataset, reference boreholes (8 in total) were selected for each climate subregion where available, and a change detection analysis was run for each borehole. After filtering the Hydstra dataset, only 14 boreholes were selected for further analysis (Figure 3 and Table 1)

The selected boreholes were chosen based on specific criteria, which include proximity to a chosen weather station, preferably in the same quaternary catchment and on the same side of a watershed (at least the latter if the former is not possible), and the number of consistent data points available (at least ten years). For each location, boreholes were statistically correlated (section 4.1), and the resulting boreholes with the most significant correlations were chosen for time series analysis. Monthly averages were used for groundwater depth. Where there was missing information, linear interpolation was used to fill in missing gaps, assuming no significant change occurred between the previous and the subsequent measurement. Every effort has been made to choose appropriate boreholes with significant data and minimal missing values. However, only a few were deemed appropriate.

Groundwater depth (GWD) was analysed in this study (in metres); however, a few assumptions about the observed GWD were made during the analysis phase of this study. The first assumption is that GWD’s seasonal patterns indicated recharge patterns and were closely compared to P in each subregion. Secondly, it is assumed that the observed GWD is due to CC, mainly P since no abstraction data are available. Lastly, suppose GWD rises, in other words, the point at which water is first encountered is closer to the surface than expected. In that case, it is assumed that the respective aquifer has experienced recharge and a resultant increase in storage.

3.1.4. Land Cover Data

National land cover (NLC) maps determined the change in land cover (LC) between 1990 and 2020. These maps were generated using Multispectral Landsat 4/5 imagery at a 30 m resolution for the 1990 and 2013/14 maps. The 2020 map is an updated version of the 2018 map and was mapped with a 20 m resolution from Sentinel 2 imagery. The 1990 and 2013/14 maps had 72 classes, whereas the 2020 map had 73 LC classes. For this study, those classes were grouped into ten categories. These categories highlight significant LC classes while simplifying the comparison and analysing LC change (Figure 2).

3.2.1. Delineation of Climate and Land Cover Subregions

The BGWMA was divided into different climate subregions based on the Köppen-Geiger climate classification (Beck et al., 2020) (Figure 3) and into different LC subregions using the latest NLC map (Figure 2). Only four subregions were chosen for discussion: BWh, BWk, Csb and Csa. These subregions were chosen because they had significant available data for analysis from which firm conclusions could be drawn. Sufficient GWD data were also available in these subregions, which could be analysed individually and correlated with other variables. Lastly, ET patterns in these subregions and strong relationships between ET and variables other than T were also considered in this selection. The boreholes analysed for the BWh subregion are technically outside the quaternary catchment and climate subregion but are on the same side of the watershed, downstream of the Gamka River. It was the only available data.

3.2.2. Change Detection Analysis: Magnitude of the Impact of Temperature, Precipitation, Evapotranspiration, and Groundwater Depth

A change detection analysis was performed using Python software for T, P, ET and GWD from selected weather stations (variables T and P), boreholes (variable GWD) and climate subregions (variable ET). This change detection analysis determined existing data trends and was used to make near-future forecasts. The models were trained on 70% of the dataset, from which a validation forecast was made and compared to the actual values of the remaining 30% of the data using four evaluation metrics: mean squared error (MSE), root mean squared error (RMSE), mean absolute error (MAE) and mean absolute percentage error (MAPE). The model with the best evaluation metrics was chosen to forecast the near future. To verify the models, they were compared with decomposition plots that distinguish between a dataset’s trend, seasonality and residual noise. Furthermore, the Mann-Kendall (MK) test is a non-parametric test used to verify observed trends in the data. In summary, for each variable, the following processes and methods were employed in this order (Figure 4):

4. Results

Climate and climate-related results (T, P, ET, GWD variables) are reported per climate subregion. LC is discussed in section 4.2. T and P trends are related to the weather station in question but representative of the subregion. GWD trends are related to the closest weather station. ET trends are spatial.

4.1. Climate and Climate-Related Trends Per Subregion

4.1.1. BWh

The BWh subregion covers 8,16% of the study area and is characterised by an arid, desert environment and hot conditions (Figure 3). The Mean Annual Temperature (MAT) is 18.6 °C and the Mean Annual Precipitation (MAP) is less than 256 mm a year. The average annual P between 2014 and 2020 was 121 mm. The LC of most of this subregion is bare land, and much less so is some grassland areas. The boreholes are situated in a valley adjacent to the Swartberg Pass, at a higher elevation than the station, and towards the South of the weather station. The aquifer is fractured with a minimal yield.

T, P, ET and Hydstra GWD data are available for this subregion. T, ET and GW data are not stationary and indicate a trend. Boreholes J2N0580 and J2N0621 indicate a difference stationarity. According to the decomposition plots (Figure 5), a slightly positive trend for average monthly T between 18 °C and 20 °C is evident, and a simple, expected annual seasonality with hot summers and cold winters is observed. A decreasing P trend is observed between 20 mm and 15 mm monthly summative precipitation. P is higher at the beginning of the year, consistent with late summer P patterns over South Africa (Kruger, 2007), after which it declines in an oscillating fashion. The ET trend is typically higher at the beginning of the year. GWD increases early in the year (February/March), peaks midway through the year and declines sharply until the end of the year or early in the following year. For some boreholes, the recharge is faster than for others, such as J2N0620. The recharge observations indicate the reliance on P as the recharge patterns follow the P season, albeit at different rates in respective boreholes.

Correlations between variables are generally not very strong (Table 2); virtually no relationship exists between T and P; however, a significantly strong positive relationship exists between T and ET. ET and J2N0582 have a slight positive relationship with P. All the boreholes are correlated, some more pronounced than others, for example, J2N0580 and J2N0582 with borehole J2N0621. J2N0621 is also firmly positively correlated with ET. For the most part, Spearman’s rank indicates slightly stronger relationships, especially among the boreholes and T and ET. Pearson displays a more robust relationship (albeit still weak) between P and ET and ET and J2N0621. These positive correlations illustrate that one variable is likely to increase with the other one or vice versa (Pearson’s), albeit monotonic or not at a constant rate (Spearman’s).

The MK test confirms the decreasing trends observed in P and ET, respectively (Table 3). T displays an increasing trend. All of the boreholes’ GWDs show decreasing trends. An increase in T from 19.5 °C to 20.5°C between 2021 and 2028 is forecasted (Figure 6a). Moreover, a decrease in P from 6 mm to 0 mm monthly totals is expected between 2021 and 2028. ET indicates a decreasing forecast from 32 mm to 20 mm between 2022 and 2036. Borehole J2N0580 is expected to decrease from -18 m to -30 m, whereas borehole J2N0582’s future forecast is at a lower constant than the validation forecast, from -12.5 m to -13.75 meters. J2N0620 show a decrease from -11.25 m and -14.25 m. J2N0621 show a new constant of -14.25 m for the future phase, 4 meters below -10.25 m, the constant for the validation phase (Figure 6b). An increase in T is highly likely in this subregion, which means that if P increases, ET should increase. In the latter case, GWD in specific boreholes could, therefore, potentially also rise. However, a decrease in P is expected, which would negatively impact ET and GWD recharge.

4.1.2. BWk

The BWk subregion, as the second largest subregion, makes up 30.53% of this study area (Figure 3). Cold, aridic desert-like conditions characterise this subregion. The MAT is lower than 18 °C and the MAP is a function of the MAT depending on the precipitation season; in this summer precipitation area, the MAP is less than ten times the MAT plus 28. The MAP for Beaufort West between 1994 and 2020 was 258.6 mm. The LC consists primarily of bare land, shrubland, and a few areas that overlap with grassland areas and irrigated agriculture. The boreholes are all in the same quaternary catchment as the weather station and on the same side of the watershed. The Hydstra boreholes are in sporadic locations, including in the valley and cut out by the Gamka River to the North (J2N0001). Various fractured aquifers in this area have differing yields.

Good quality T and P data observed on a daily time-step for over 20 years with minimal missing values are available for this subregion. The same holds for ET data, but data have only been available since 2009. Good quality GWD data from the NGA and Hydstra datasets has been observed for multiple boreholes in the Beaufort West region. Beaufort West was analysed for all four variables and was used as the example modelling station. Beaufort West T, ET and GWD data for nine of the ten boreholes are not stationary. T, ET and four boreholes display trend stationarity.

According to the decomposition plots (Figure 7), Beaufort West indicates a gradual increase in T from 17 °C to 19.5 °C. An annual T seasonality, higher Ts in the summer and low Ts in the winter is observed. P varies between 20 mm and 55 mm over the entire period. It displays a robust seasonal component with relatively higher P in the summer months (December to February), followed by a sharp decrease and an oscillating climb from May/June until December again. ET is displaying an overall decreasing trend, with a similar seasonality observed before; a steady climb into the summer months, oscillating behaviour until March and then a decline until mid-year. NGA boreholes show a similar seasonal pattern. The year starts, and the GWD increases until February/March, after which it decreases until about August/September, although this pattern’s starting and ending times vary slightly in different boreholes.

Hydstra GWD datasets’ seasonality behaviours are similar, with an increase usually at the beginning of the year for a couple of months and then a decrease the rest of the year. J2N0001 has a simple seasonality with a sharp increase from December to mid-year, after which an immediate decrease for the rest of the season follows. J2N0043 display an increase from the end of the previous year, little variability for a couple of months and a sharp decrease for the rest of the season. The recharge lag time is relatively far behind the P peak, which could indicate a slow recharge. The GWD of 3222BC00179 is expected to rise and is situated 4.71 km Southwest and at a lower elevation than the station. H2N0041 and J2N0043 are 10.1 km N and 9.1 Northnorthwest from the station, respectively, in the valley at a higher elevation. This could indicate that the boreholes’ recharge lag time in the valley is longer.

A robust positive relationship between T and ET and a weak correlation between T and P (Table 4 and 5) is evident. The correlation between T and ET is the same as Pearson’s and Spearman’s, indicating a linear yet monotonic relationship. Pearson’s rank shows a more robust relationship between T and P and P and ET, which indicates a linear positive relationship (Table 5); P and ET are associated with warmer weather. All the boreholes display strong correlations except correlations with 3222BC00152 and 3222BC00179, although the latter correlation is still significant. A relatively significant negative relationship exists between T and 3222BC00152 and J2N0019, respectively. Spearman’s correlations are stronger among the Hydstra boreholes, whereas between the NGA boreholes, the correlations are more linear except for 3222BC00152. The MK test (Table 6) shows an increasing trend in T and borehole 3222BC00179 and a decrease in P, ET and all Hydstra boreholes.

The statistical models that proved the most suited include linear regression for T, P, and ET (Figure 8a). GWD models are a mix of all four models. T is expected to increase slightly from 18.75 °C in 2021 to 19.5 °C in 2040. P is expected to decrease from 15 mm to 8 mm monthly summative P for the same period as the T. Between 2021 and 2036, ET is expected to decrease from 40 mm to 32 mm per month. GWD from ten boreholes were modelled respectively, of which five display no significant trend (3222BC00170, 3222BC00152, 3222BC00171, J2N0019 and J2N0550). 3222BC00179 shows a rising trend according to the linear regression model from -13 m to -7 m between 2006 and 2028. GWD for boreholes J2N0001, J2N00041, J2N0043 and J2N0618 is expected to decrease from -23 m to -43 m, -35 m and -53 m, -30 m and -50 m and -32 m and -41 m respectively between 2020/22 and 2036/40 (Figure 6b,c–d).

Given the results, an increase in T is highly likely in the near future of this subregion, which could drive the process of ET at a more robust rate, provided there is sufficient P. However, P is likely to decrease, which could negatively impact ET. GWD recharge time seems to follow P seasons; however, it is at a much longer lag than previously observed.

4.1.3. Csb

Towards the West is the Csb subregion, which covers an area of 10.57% of the BGWMA (Figure 3). It also has a temperate climate but with dry and warm summers. The P in the driest month of summer is less than 39 mm, and a third of the P in the wettest month is still more significant than in the driest summer month. The LC consists of vast areas of dryland agriculture, spots of irrigated agriculture (minority), shrubland, grassland, bare land, indigenous forests and urban areas, especially towards the coast. Large water bodies are also observed. The weather station and borehole 3419AD00004 are in the same quaternary catchment and on the same side of the watershed. The borehole is situated on a golf course about 0.5 km from the ocean at its closest point; however, it is no longer used and was a fractured aquifer.

Good quality T and P data have been recorded daily for over 20 years. GWD data are also scarce in this subregion; however, some data were available for Hermanus, but only between 1956 and 1984. Nonetheless, the data were used for analysis, and Hermanus was chosen as the model station for this region. All data were stationary except for ET and GWD, exhibiting trend stationarity. A typical T seasonality is observed for Hermanus (Figure 9). P seasonality is almost opposite that of T; low P in summer and high P in winter. The same seasonality pattern is evident for ET as was previously described. GWD seasonality is opposite to T: low in summer and high in winter. GWD seasonality follows P: low in summer and high in winter. Recharge is relatively rapid, consistent with the rainy season, which illustrates the importance of recharge in the area for GWD recharge.

Relatively strong negative relationships have been confirmed between T and P for Hermanus (Table 7). This relationship and the negative correlation between P and ET are more robust with Spearman’s rank, indicating a negative monotonic relationship; P is associated with colder T, and ET is more prominent in warmer Ts. A solid positive and more linear relationship between T and ET confirms the latter observation. A negative relationship is evident between P and ET. The MK test (Table 8) confirmed an increasing trend for T, P and ET and no trend for GWD. The linear regression forecasting model proved to be the best model for the Hermanus area (Figure 10a,b). Between 2021 and 2040, T displays an insignificant trend from 17.4 °C to 17.65 °C, and P is expected to increase from 42.5 mm to 48 mm. ET shows an increase between 2022 and 2036 from 63.4 to 64.5 mm per month at a constant rate. GWD indicates no significant trend.

If significantly elevated T levels should occur in the future, ET is highly likely to increase, provided sufficient P. Although no significant GWD trend is observed in this dataset, P perturbations are highly likely to influence GWD since recharge follows P seasonality.

4.1.4. Csa

This temperate subregion with dry and hot summers (Figure 3), Csa, is similar to the climate of Csb. The MAT is 18.7 °C and the P in the wettest month (June) is more than three times that of the P in the driest month. This winter P area receives a MAP of 272.31 mm (between 1999 and 2020). Csa subregion covers an area of 1.18%. The LC consists of shrubland, mostly irrigated agriculture, a large waterbody, some urban areas (Worcester) and a few relatively small areas of dryland agriculture. Three of the boreholes are in the same quaternary catchment (30959/H1N0018, H1N0055), and all of them are on the same side of the watershed at more or less the same elevation, between 204 and 209 m. 30959/H1N0018 are adjacent to the Breede River, next to the riverbed, between the Breede River and Brandvlei dam. 3319CB00003 is between two tributaries on bare land. H1N0055 is the closest to the weather station, just over a kilometre on a farm. H2N0521 in the industrial area. 3319CB00015 is on a farm near Breede River and is an intergranular aquifer of high yield. All other aquifers are fractured aquifers.

All four variables were analysed in this subregion, and good quality T, P, ET and GWD data were available. T, P data were available from 1997 to the present, ET data from 2009 to the present, and NGA GWD data between 1982 and 1992. Three Hydstra GWD datasets were available for analysis, of which H1N0018 is a continuation of 30,959 (data gap of 6 years between 1993 and 1999). Worcester was used as the modelling weather station. The data are stationary according to the ADF and KPSS tests, except for a trend component for T and ET. Hydstra datasets indicate difference stationarity.

According to the decomposition plots (Figure 11), no clear trends are seen for T and P, but ET shows an increasing trend. The T seasonality is the typical sine-like curve. In contrast, P seasonality is almost opposite to T, with low P in the summer months and relatively higher P in winter in the middle of the year. The increase in P is rapid just after the new year until about April, after which it declines at a gradual pace. However, with fewer variations, a similar seasonality pattern for ET is observed for the other subregions. GWD seasonality is essentially the same for all three boreholes; a rising pattern is observed until April, then lowering to July/August, and then slowly rising again. Hydstra datasets indicate that groundwater starts the year with a decline, reaches its lowest point in February/March, rapidly increases, varies at its peak and then declines in the second half of the year until just after the new year. GWD recharge following the winter P season is relatively slow. The recharge lag time for H1N0018 and H1N0055 is about 6 to 7 months after the rainy season, whereas the other boreholes, between 5 and 10 km from the station, all show a faster recharge time.

This study’s most robust negative relationship between T and P is seen here (Table 9 and Table 10). A significantly strong positive linear relationship between T and ET and a negative relationship between ET and P is observed. The relationships between T and P and P and ET are more strongly negative, with Spearman’s rank indicating a non-linear monotonic relationship. T has a slightly negative relationship with H1N0055, which is more robust with Spearman’s rank, i.e., not linear. P has a slight positive relationship with H1N0018 and H1N0055. H1N0018 and H1N0055 have a strong positive relationship. The other Hydstra GWD relationships are positive but not as robust as between H1N0018 and H1N0055. The NGA boreholes are all significantly positively correlated. Since all NGA boreholes are positively correlated and 30,959 is the old H1N0018, the assumption is that most boreholes, except for H2N0521, are positively correlated and follow similar patterns. Therefore, a slight positive relationship between P and borehole H1N0055 could mean P is influential in GWD recharge.

The MK test (Table 11) results reveal a decreasing trend for T and borehole H2N0521 and an increasing trend for ET and boreholes 30959, 3319CB00015 and 3319CB00003. The models that best fit the variables in this region are a mix of linear regression, Prophet and Auto ARIMA (Figure 12a–c). Between 2021 and 2040, this model indicates a significant negative trend in T from 19 °C to 18 °C. No trend for P is observed. ET, however, increases from 67 mm to 71.2 mm between 2022 and 2036. NGA boreholes 30959, 3319CB00015 and 3319CB00003 all display rising trends; however, the projections could only be made up until 2002 (from -2.4 m to -2.25 m, -1.75 m to -0.75 m, and -1.75 m to -1.48 m respectively). Hydstra boreholes H1N0018 and H1N0055 display no significant trends, whereas H2N0521 lowers from -0.04 m to -0.09 m between 2022 and 2040. T is expected to decrease in this area, and P is expected to increase, which could result in lower levels of ET. Furthermore, recharge could be P-dependent; therefore, the higher the P, the higher the resulting potential recharge, especially in the vicinity of borehole 3319CB00015.

4.2. Land Cover Trends

LC change was determined from the NLC maps for 1990, 2013/2014, and 2020. From 1990 to 2020, the study area observed an increase in urban, dryland and irrigated agriculture, water surfaces, grassland and shrubland. The most significant change since 1990 has been an increase in shrublands by 8059.53 km². Commercial forests, natural forests, wetlands and bare lands have decreased. The most significant decrease is in bare lands with 2573.99 km² (Table 12).

The BWh and BWk subregions show anticipated decreases in ET trends, whereas the Csb and Csa subregions show anticipated increases in ET. Table 13 provides a percentage summary per LC for each 1990, 2013/14 and 2020 subregion. Between 1990 and 2013, the BWk subregion indicated a substantial decrease in indigenous forests, bare land, and wetlands, as well as a decrease in dryland agriculture and grass*. A decrease in bare land and an increase in shrublands, water surfaces**, urban areas, plantations, and irrigated agriculture is a positive sign of development in this subregion. The BWh subregion indicates losses in spatial areas of forests, plantations, dryland agriculture, water surfaces*, wetlands and grasslands*. Bare lands also decreased. Urban areas, shrublands, and irrigated agriculture show increases, indicating that development is occurring. Signs of desertification include grasslands being replaced by shrublands. In the Csb subregion, an increase in urban areas, indigenous forests, dryland and irrigated agriculture** is observed. In contrast, a decrease in plantations*, bare land, water surfaces*, wetlands, grass*- and shrublands is evident. Loss of bare land, urban increase and agricultural development suggest development in this subregion. In the Csa subregion, increases in urban areas, water surfaces, bare land, dryland and irrigated agriculture* and losses in forests, plantations*, wetlands, grass- and shrublands* have been observed. Development is also evident in this subregion. The loss of wetlands across the study and forests in BWk, BWh and Csa areas are highly likely due to urbanization and developments in irrigated agriculture.

* decrease since 1990, but increase since 2013/14, ** increase since 1990, but decrease since 2013/14.

5. Discussion

Except for one subregion, the Csa subregion, T is expected to increase by 2040 between 0.25 °C and 2.5 °C. The more significant increases are seen in the drier areas, such as the BWh and BWk subregions. These T trends are confirmed by previous studies that projected T increases by up to 2,5 °C by 2060 in Africa (Kahsay et al., 2019; Gemitzi et al., 2017; Karmalkar et al., 2010; McSweeney et al., 2010; Ragab & Prudhomme, 2001). There are only two instances where T and GWD are correlated; in both cases, they are negative. These instances include borehole J2N0019 in the BWk subregion and borehole H1N0055 in the Csa subregion. In both cases, the correlations are not significant. However, they are indicative enough of some underlying relationship between T and GWD, which essentially means that an increase in T could indirectly and negatively impact GWD (J2N0019) and vice versa (H1N0055), possibly due to elevated T enhancing the process of ET. This effect of elevated T on lowering the GWD due to enhanced ET has been found in other parts of the world (Salem et al., 2017; Gunawardhana & Kazama, 2012; Ranjan et al., 2006). The latter finding depends on an area’s degree of vegetation cover (Owuor et al., 2016) and the abundance or lack of P.

P patterns are highly variable; however, P is generally expected to decrease in most subregions between 5 and 10 mm monthly summative P. Certain areas show no anticipated near-future trends for P (for example, Csa). Places like the Csb subregion show an anticipated 5 mm increase in P. Coastal areas are less vulnerable to anticipated negative P trends. In contrast, drier and hotter areas are seeing declines in P by up to 10 mm monthly summative P (BWh – 5 mm, BWk – 7 mm). Since P is not vast in these drier areas, the P margins are already small, so a mm drop is alarmingly significant. Other studies have confirmed these P trends; increases are expected in some temperate regions and decrease in others, usually arid/semi-arid areas (Kahsay et al., 2019; Gemitzi et al., 2017, Burke and Stott, 2017; Nissen and Ulbrich, 2017).

No significant relationship between P and GWD could be statistically established; however, recharge seasonality seems to follow P seasonality in all areas with differing lag times. The lag times seem faster in more temperate areas like Csb and Csa than the drier arid/semi-arid areas like BWh and BWk. It has been found that there is a close relationship between P and groundwater storage potential and recharge potential is also highly dependent upon P, which will naturally influence GWD. Moreover, groundwater droughts have also become more common, especially in some temperate regions with an increase in hot months. However, the primary controller of groundwater droughts is decreased P. These findings are consistent with the findings from this study.

ET is expected to decrease between 2 mm and 12.5 mm in the drier areas on average monthly. In the more temperate climates (Csb, Csa), an increase of ET between 1.1 mm and 4.5 mm is expected. A close interworking relationship exists between T, ET, and P in the more temperate areas (Csb and Csa); T fuels the ET process, which does not directly determine P but plays a role in generating P. However, minimal changes in T and/or a reduction in P are expected in these areas, which will negatively impact ET. Moreover, T and ET display the same respective seasonal patterns: high in summer and low in winter, with some variations for ET at maximum and minimum peaks, confirming this correlation across all study areas. The BWh subregion is the only one to indicate significant positive relationships between ET and boreholes J2N0580 and J2N0621, the latter being more robust. It is inconclusive whether ET influences GWD or not; however, ET does seem to affect GWD in more arid/semi-arid environments like BWh. This observation is confirmed by previous studies that illustrate that higher ET levels could contribute to water loss, groundwater drought (Osei et al., 2019), and increased abstraction. However, ET levels could also be lowered by deforestation, which has been observed in subregions BWk, BWh and Csa (Berhail, 2019; Costa et al., 2003). The results of deforestation would be site-specific and need to be analysed on a local level.

ET has been confirmed to be influenced by LULC (Luo et al., 2018; Tamm et al., 2018; Zhang et al., 2018). The more arid study areas to the North/Northeast all indicate anticipated decreases in ET, which amount to a total area percentage of 38.7%. The more coastal, Southern and Western interior areas indicate anticipated increases in ET, which makes a total area percentage of 10.8%. However, coastal areas are usually more humid, which could bring ET levels down, especially with elevated levels of CO₂, and offset the effect of T (Guo et al., 2017; Snyder, 2017). The moderating effects of the coast on T have been verified in this study in the coastal areas of the Csb subregion, which is only projected to see a relatively minor increase in T but also a projected increase in P and a minor increase in ET.

GWD mostly show no significant increasing or decreasing trends. Of the 21 boreholes analysed, five show increasing trends between 2 and 4.5 m, eight show decreasing trends between -0,05 m and -22 m and eight show no trends. In the more temperate areas, only seven boreholes were analysed, of which three indicated ‘no trend’ and three indicated increasing trends (but these were old datasets and only forecasted until 2002). One indicated a decrease of 0.05 m. Six of these seven boreholes are situated in the Csa subregion. The recharge lag time after the rainy season is usually slower in the drier areas than in the more temperate areas. With elevated T levels in the more arid areas and ET levels, as a result, it is highly likely to influence recharge in these areas (refer to a previous paragraph). In most cases, the boreholes are all positively correlated, especially in the BWk subregion; however, this is not always the case. In the BWh and Csa subregions, the boreholes are located all over the subregion at different elevations and locations, which could also influence GWD.

In terms of LC, it is evident, in general, that there is large-scale development in the entire study area due to increases in urban areas, dryland and irrigated agriculture and the loss of bare land. Water surfaces and grass- and shrublands also seem to be increasing whilst wetlands and indigenous forests are decreasing. The latter could be due to the study area’s drier climate, lands cleared for development, or both. Increases in water surfaces could be artificial, such as the building of dams. Visual observations from Figure 2 confirm the loss of bare land and urban and agricultural development in the study area. Deforestation is also apparent due to anthropogenic activities, particularly increasing population and urbanization. Deforestation could cause changes in surface albedo and wind patterns and reduce the Leaf Area Index (LAI) and rooting depth, consequently leading to lower ET levels (Berhail, 2019; Costa et al., 2003). These trends could be the case for subregions BWk, BWh and Csa, which all show deforestation trends. The rate of deforestation is increasing, which is evident in the data (BWk); the decrease between 2013/14 and 2020 is much more significant than between 1990 and 2013/14. In some instances, forests increased between 1990 and 2013/14 but rapidly decreased since then (BWh and Csa).

Deforestation could potentially lead to increases in groundwater recharge in the area, provided there is sufficient P since it has been found that a reduction in the LAI can increase recharge (Owuor et al., 2016). This recharge increase is confirmed by Varet et al. (2009), who found that the rate at which groundwater is lost to ET is a function of rooting depth; thus, changes in water table levels due to deforestation (and the introduction of grasslands) should decrease ET rates. GWD usually correlate well with grass, crop and shrubland but not with forests (Bloomfield et al., 2019). On the contrary, places with increased vegetation density (for example, afforestation) can potentially cause a decrease in recharge, for example, in the Csb subregion. Moreover, forests play a role in the hydrological cycle by regulating atmospheric Ts, as seen in the relatively low anticipated T trend in Csb. The loss of wetlands for the entire area and all four subregions is due to large-scale development like urbanization and irrigated agriculture. Urbanization, conversely, can cause decreases in ET (Luo et al., 2018; Minnig et al., 2018) and a decrease in recharge; however, the influence on recharge can vary widely. The impact of LC on recharge will need to be assessed on a local level.

It has been found that grasslands being replaced by shrublands indicate desertification, which is caused by agricultural development such as dryland agriculture and drought. Most projections expect this trend to occur in which shrubland replaces grassland. In general, the study saw an increase in grasslands; however, in the respective subregions, grasslands are decreasing. The BWk, BWh and Csb subregions saw a decrease in grasslands between 1990 and 2013/14 and a rapid increase since then until 2020. The opposite is true for the Csa subregion. Shrublands are also decreasing in the Csa and Csb subregions but not in the BWk and BWh subregions. This trend, favouring shrub growth, is due to an increase in minimum Ts and a decrease in frost, which usually inhibits shrub growth. Furthermore, shrubs generally fare better in droughts due to deeper roots and mobilization of deeper soil water. Higher concentrations of CO₂ (also due to deforestation) accelerate the growth rate of shrubs. LULC has also been found to influence grass-to-shrub changes in semi-arid environments, especially in the case of high stocking rates and overgrazing (Acocks, 1953; Peters et al., 2006).

Model Performance and Evaluation

Ample data were available for this study area; however, sufficient correlation conclusions could not be made in each subregion due to the lack of groundwater data and too short or outdated datasets. However, the best use of the data was made to assess the groundwater situation of the BGWMA. Most of the datasets were non-stationary, often exhibiting an underlying trend. Furthermore, according to the decomposition plot, all the datasets show a strong seasonal pattern, which is expected from climate and climate-related data. For this study, the seasonal pattern was removed since the goal was not to model seasonal patterns but underlying trends.

Regarding the models’ performance and evaluation, an effort was made with each of the 33 models across all variables to apply the slightest adjustments (via rolling mean and/or seasonal differencing) to eliminate the seasonal pattern and model the underlying trend if present. This essentially meant that non-stationary datasets were used for modelling. The evaluation metrics were calculated in the validation phase. If no precise model fit was evident, the entire dataset was used to validate and choose an appropriate model or to confirm the model selection of another. T was the least complicated variable to model, followed by ET. P was highly variable and difficult to model. GWD was challenging to model in some areas like the BWk subregion and relatively easy in the Csa subregion. Lastly, 100% of the T, P and ET models and 85% of the GWD models had excellent MAPE values. Ninety-four percent of the models indicated a trend consistent with the MK test results, and 97% were consistent with the decomposition plots. These results indicate that the models are reliable in making near-future forecasts.

6. Conclusions and Recommendations

The results of this study prove that there are areas where anticipated trends in temperature, precipitation, evapotranspiration, groundwater depth and land cover should be considered in groundwater management. Moreover, these variables need to be analysed more locally to reach a satisfactory conclusion about water security, especially groundwater resources, in the BGWMA. The study also exposed not the lack of groundwater depth data but the lack of appropriate groundwater depth data necessary to draw firm conclusions. Although this study’s results are inconclusive, it is a first step to fully understanding the area’s underlying climate, groundwater and land cover change trends.

The more arid areas to the North and Northeast of the study area, BWh and BWk, are expected to become hotter and drier, which could be detrimental to communities dependent on the sole use of groundwater. Losses in wetlands, water surfaces, and indigenous forests, as well as increases in shrublands in both areas, indicate that the climate is getting more arid. Loss of vegetation, an increase in temperature, decreases in evapotranspiration and decreases in precipitation will likely result in decreased groundwater depth since groundwater depth recharge depends on precipitation, and recharge seasons follow precipitation seasons with differing lag times. The recommendation would be to monitor groundwater depth more intensely, especially in hotspots, and to limit groundwater use as much as possible. A particular focus on the development of irrigated agriculture will change the rooting depths of plants and could aid in the process of evapotranspiration. Development plans should also adopt water-efficient technologies, proper land management, avoiding overgrazing, and rainwater harvesting during wetter seasons (summer in these areas). Moreover, by planting crops that mature earlier and shortening growing cycles, water deficits and agricultural water use could be mitigated to some extent. Other measures include increasing irrigation efficiency, reducing evaporation losses through mulching or cover cropping, increasing soil water retention through reduced tillage and/or cover cropping, or enriching soil organic matter contents. These measures should be based on ‘transient numerical modelling’.

The Csb and Csa subregions are towards the West of the study, under a more temperate climate and show more favourable outcomes regarding near-future climate forecasts. In short, temperature is expected only to increase slightly, even decrease in other areas, and precipitation is expected to stay the same or increase. Robust relationships are observed between temperature, precipitation and evapotranspiration in these areas. groundwater depth recharge is also more rapid after the precipitation season, which illustrates the dependence of groundwater depth recharge on precipitation; therefore, groundwater depth could potentially rise if precipitation increases. The loss of wetlands, forests, grasslands and shrublands is due to large-scale development in the area. Signs of more favourable climatic conditions include increased water surfaces, dryland agriculture and indigenous forests. Even though the future in terms of precipitation and groundwater depth is slightly more optimistic in these temperate areas, more groundwater monitoring is recommended, especially in hotspots. In the Csa subregion, a particular focus should be on the further development of irrigated agriculture and focused recharge since no increases in precipitation are expected. A limit should be placed on drilling new boreholes and issuing licenses until there is a better understanding of groundwater storage in the area.

Future recommendations would be to examine the underlying correlations between the variables more intensely utilizing a multivariate analysis. The relationship between land cover, evapotranspiration and groundwater depth needs a particular focus. Each subregion needs to be studied intensely, and areas of interest in the respective subregions, such as where rapid urbanization or deforestation occurs, need to be identified and studied. Data, including water use or abstraction, are essential to separate climate change effects from abstraction.