A Novel Approach for Optimum Conjunctive Use Management of Groundwater and Surface Water Resources under Uncertainty

Gislar Edgar Kifanyi

Department of Civil Engineering, Mbeya University of Science and Technology, Mbeya, Tanzania.

DOI: 10.4236/ojmh.2023.131003                

Abstract

Uncertainty in determining optimum conjunctive water use lies not only on variability of hydrological cycle and climate but also on lack of adequate data and perfect knowledge about groundwater-surface water system interactions, errors in historic data and inherent variability of system parameters both in space and time. Simulation-optimization models are used for conjunctive water use management under uncertain conditions. However, direct application of such approach whereby all realizations are considered at every-iteration of the optimization process leads to a highly computational time-consuming optimization problem as the number of realizations increases. Hence, this study proposes a novel approach—a Retrospective Optimization Approximation (ROA) approach. In this approach, a simulation model was used to determine aquifer system responses (draw-downs) which were assembled as response matrices and incorporated in the optimization model (procedure) as coefficients in the constraints. The sample optimization sub-problems generated, were solved and analyzed through ROA-Active-Set procedure implemented under MATLAB code. The ROA-Active Set procedure solves and evaluates a sequence of sample path optimization sub-problems in an increasing number of realizations. The methodology was applied to a real-world conjunctive water use management problem found in Great Letaba River basin, South Africa. In the River basin, surface water source contributes 87% of the existing un-optimized total conjunctive water use withdrawal rate (6512.04 m³/day) and the remaining 13% is contributed by groundwater source. Through ROA approach, results indicate that the optimum percentages contribution of the surface and subsurface sources to the total water demand are 58% and 42% respectively. This implies that the existing percentage contribution can be increased or reduced by ±29% that is groundwater source can be increased by 29% while the surface water source contribution can be reduced by 29%. This reveals that the existing conjunctive water use practice is unsustainable wherein surface water system is overstressed while groundwater system is under-utilized. Through k-means sampling technique ROA-Active Set procedure was able to attain a converged maximum expected total optimum conjunctive water use withdrawal rate of 4.35 × 10⁴ m³/day within a relatively few numbers of iterations (6 to 8 iterations) in about 2.30 Hrs. In conclusion, results demonstrated that ROA approach is capable of managing real-world regional aquifers sustainable conjunctive water use practice under hydro-geological uncertainty conditions.

Keywords

A Novel Approach, Conjunctive Use Management, Retrospective Optimization Approximation (ROA), Uncertainty

1. Introduction

Water availability is a major determinant for any social-economic development. Climate change, rapid population growth, increasing water demand, mining and industrial development coupled with the dynamics in the hydrological cycle as well as the hydrogeological conditions of river systems have led to increased pressure on water resources in many developing countries [1]. One of the new methods currently accustomed for addressing water shortage problems in most developed countries is optimum conjunctive water use [2] [3]. However, optimum conjunctive use demands that the surface and subsurface reservoirs are thoroughly characterized if deterministic methods are to give accurate results. This is because in real world phenomena full characterization of surface and subsurface reservoirs is practically impossible. Surface and subsurface hydrology is highly uncertain [4] [5]. The parameters of uncertainty involved pose a great challenge for decision makers to have reliable water use policy and predictable returns from the water systems [6] [7] emphasized that water managers and decision makers must understand the surface water-groundwater system interactions, especially under uncertainty. This is because as in any system management, conjunctive water use system is also characterized by uncertainties. Hence, it is imperative that this uncertainty has to be explicitly considered when managing conjunctive water use [5] [8].

Different methods have been developed and used by various researchers to account for the uncertainty. Ndambuki [8] investigated Multi-Objective groundwater quantity management using stochastic methods for modeling cost of groundwater production and pumping rate decisions. The uncertain conditions of aquifer hydrogeology were analyzed for many hydrogeological model simulations in a Monte Carlo simulation and Second-Order-Cone Programming approaches. Schoups, et al. [9] used a multi-objective interannual optimization model to generate Pareto curves and evaluate trade-offs between sustainable agriculture and optimal reservoir operation. The uncertain conditions of hydrology were analyzed for many equally probable streamflow series in a Monte Carlo approach. Marques, et al. [10] presented modeling conjunctive water use operations and farm decisions with two-stage stochastic quadratic programming. Schoups, et al. [11] studied conjunctive water use management to alleviate drought for irrigated agriculture with a spatially distributed simulation-optimization model. Moreover, Kentel, et al. [12] developed a simulation-optimization model with constraints on drawdowns which were used to optimize additional groundwater withdrawal rates in a coastal aquifer system. Other studies include [13] [14] [15] also investigated optimum groundwater operations with a stochastic surface water system.

From the aforementioned literature, it is evidenced that simulation-optimization models are normally used for conjunctive water use management. However, note that direct application of the approach where all realizations are considered for every iteration of the optimization process leads to a highly computational time-consuming optimization as the number of realizations increases [16]. Moreover, Wang, et al. [16] opined that handling of uncertainty, mainly those arising from hydrogeological conditions, poses a major task in the application of optimization methods. Thus, they underscored the need for many geological realizations if high degree of uncertainty is to be resolved. Hence, the aim of this work is to propose an efficient method for managing real-world regional aquifers conjunctive water use withdrawal over multiple geological model realizations. This paper is based on the most recent emerged technique of Sample Average Approximation (SAA) method—the Retrospective Optimization Approximation (ROA) approach a stochastic simulation-optimization framework.

Traditionally, integration of simulation model with optimization-based management procedure (model) is a challenge and difficult task, takes substantial computational time to realize optimal solution [7]. Often embedding or response matrix approaches are used to combine simulation and optimization models [17]. In this work, response matrix approach was adopted. For more detail discussions on embedding and response matrix approaches one is referred to [2] [8] [17] [18] [19] and [21].

In short, the principal idea behind ROA approach is that it does not consider all realizations (samples) at every iteration of the optimization process. Alternatively, ROA procedure generates a sequence of approximate sample optimization sub-problems, which sequentially account for increasing number of realizations (sample size). ROA approach takes in the advantage of “warm start” method where an initial guess is updated by assigning the initial solution of the current sub-problem being the solution obtained from the previous sub-problem solution solved. ROA method is new approach in solving regional aquifers conjunctive water use management problems under uncertainty through simulation-optimization modelling framework. However, the author is aware that in current years, the approach has mostly been applied in other areas such as operation research [22] [23], petroleum engineering [24] [25] and well placement optimization problem [16] [26] [27] [28] [29]. Recently, the method has been used in water resources management through hypothetical example [4] [5] and groundwater resources management [30] but, their applications were mainly based on fixed and few decision variables. The ROA procedure developed and applied in this study is designed to solve a multiple variables regional aquifers conjunctive water use management problem under hydrogeological uncertain conditions.

2. Materials and Methods

2.1. Formulation of Conjunctive Use Optimization Problem

In the conjunctive use optimization problem, we seek to optimize the surface water and groundwater withdrawal rates of a number of spatially distributed pumping wells and abstraction points (surface water diversion points), considering aquifer hydraulic conductivity field (geological) uncertainty conditions. The uncertainty was addressed as assemblage of random aquifer system drawdowns (responses) due to unit withdrawal rate at every pumping well location. The assemblage set of the random responses (i.e., response matrix coefficients) represents simulation model in the simulation-optimization stochastic modeling framework.

Consider there exists a probability space $\left(\Omega ,F,P\right)$, where Ω, is a set of entirely possible outcomes;F is a set of collection of all subsets of Ω possible realizations; and P defines a probability measure. Thus, for a given total possible number of outcomes Ω, assume hydraulic conductivity field uncertainty realization be defined by $\varpi$, such that $\varpi \in \Omega$. Hence, for every realization $\varpi$ there will be different set of random response matrix denoted by $M_{\varpi }$ allied with random response matrix components represented as $a_{ij,\varpi }$. In general form, therefore, the conjunctive use optimization problem under the hydraulic conductivity field uncertainty realizations can be formulated as follows: Given a solution setY, such that $Y^{*}\in Y$, find a solution $Y^{\text{*}}$, that:

$\begin{array}{c}\text{Maximize}\end{array}F\left(Y\right)=\left[{\displaystyle {\sum }_{j=1}^{N_{gw}}{Y}_{gw,j}}+{\displaystyle {\sum }_{d=1}^{N_{sw}}{Y}_{sw,d}}\right],j=1,2,\cdots ,N_{gw};d=1,2,\cdots ,N_{sw}$ (1)

Subject to:

${\mathbb{E}}_{\Omega }\left[G\left(Y_{gw,j},M_{\varpi }\right)\right]\le b_i,i=1,2,\cdots ,N_c;\forall \varpi \in \Omega$ (2)

where ${\mathbb{E}}_{\Omega }$ represents the expectation function over the solution set of all possible outcomes Ω; and G is a numerical stochastic modeling process that computes the sample optimization sub-problems constraint functions $G\left(Y_{gw,j},M_{\varpi }\right)$ for a given $Y$ and the realization outcome $\varpi$. In which, $Y$ is a vector in space ${ℝ}^N$ such that $N=N_{gw}+N_{sw}$ defining total number of decision variables of the vector variables $Y_{gw}$ and $Y_{sw}$ such that $Y_{gw},Y_{sw}\in Y$. In which, $Y_{gw}$ and $Y_{sw}$ are groundwater withdrawal rates (a vector in space ${ℝ}^{N_{gw}}$ of the first elements of) and surface water withdrawal rates (a vector in space ${ℝ}^{N_{sw}}$ of the remaining elements of Y) respectively. $F\left(Y\right)$ is objective function evaluated over estimates of a random constrain function $G\left(Y_{gw,j},M_{\varpi }\right)$ by performing numerical simulations with the hydrogeological model defined by realization outcome $\varpi$ of uncertain parameter. Note that in the formulation (1) through (2), the vector variable $Y_{gw,j}$ denotes the groundwater pumping rate of a pumping well located at j; the vector variable $Y_{sw,d}$ is the surface water abstraction (withdrawal) rate at a point located at d; $b_i$ is a hydraulic head constraining value at control point i; $N_{gw}$ and $N_{sw}$ present total number of groundwater pumping wells and surface water diversion points, respectively; and $N_c$ is the total number of monitoring control points.

Recall the random response matrix denoted by $M_{\varpi }$ with allied random response matrix components $a_{ij,\varpi }$ or simply denoted by $\bar{\varpi }$. Assume randomness in hydraulic head drawdowns is only due to uncertainty conditions of aquifer system hydraulic conductivity field realizations and that $a_{ij,\bar{\varpi }}\in M_{\bar{\varpi }}$ such that $M_{\bar{\varpi }}\in \bar{\Omega }$. Moreover, assume that P is well defined with unknown distribution function but, what is known is expected mean values, and/or standard deviation of the random responses $\bar{\varpi }$. Consider a stochastic constraint function process as $G\left(Y_{gw},\bar{\varpi }\right)\in \bar{\Omega }$ to be defined as $G\left(Y_{gw},\bar{\varpi }\right)=M_{\bar{\varpi }}{Y}_{gw}$. Thus, the expected value of the function $G\left(Y_{gw},\bar{\varpi }\right)$ is defined as ( [2] [4] [5] ):

$\mathbb{E}\left[G\left(Y_{gw},\bar{\varpi }\right)\right]=\mathbb{E}\left[A_{\bar{\varpi }}{Y}_{gw}\right]$. (3)

Hence, the constraint inequality (2) can be formulated as $\mathbb{E}\left[G\left(Y_{gw},\bar{\varpi }\right)\right]\le b$ whereby the expectation $\left[G\left(Y_{gw},\bar{\varpi }\right)\right]={\displaystyle {\int }_{\bar{\Omega }}G\left(Y_{gw},\bar{\varpi }\right)\text{d}P\left(\bar{\varpi }\right)}$ is the corresponding expected value function. Consequently, inequality (2) can be approximated by applying Monte Carlo sampling based estimation methods (e.g., the ROA approach) by taking into account a sequence of finite set of generated samples of random response matrices for $\mathbb{N}$ realizations of $\left\{M_{{\varpi }_{{\mathbb{N}}_k}}\right\}=\left\{M_{{\bar{\varpi }}_{{\mathbb{N}}_1}},\cdots ,A_{{\bar{\varpi }}_{{\mathbb{N}}_k}}\right\}$. The approximate the expected constraint function $\mathbb{E}\left[G\left(Y_{gw},{\bar{\varpi }}_{{\mathbb{N}}_k}\right)\right]$ as ( [2] [4] [5] ):

$\mathbb{E}\left[G\left(Y_{gw},{\bar{\varpi }}_{{\mathbb{N}}_k}\right)\right]=\frac{1}{{\mathbb{N}}_k}{\displaystyle {\sum }_{z=1}^{{\mathbb{N}}_k}G}\left(Y_{gw},{\bar{\varpi }}_{{\mathbb{N}}_k}\right)$ (4)

Then progressively evaluate resulting sample optimization sub-problems of formulation (5) through (6) by using ROA method for $k=1,2,\cdots ,{\mathbb{N}}_{sp}$ (in which ${\mathbb{N}}_{sp}$ is total number of sample optimization sub-problems). Thus, by substituting the inequality constraint (2) with the approximations of expected constraint function, the corresponding Estimates Retrospective Conjunctive water use Sample Optimization Problems (ERCSOP) can be formulated as:

$\begin{array}{c}\text{Maximize}\end{array}{F}_{{\mathbb{N}}_k}\left(Y\right)=\left[{\displaystyle {\sum }_{j=1}^{N_{gw}}{Y}_{gw,j}}+{\displaystyle {\sum }_{d=1}^{N_{sw}}{Y}_{sw,d}}\right],j=1,2,\cdots ,N_{gw};d=1,2,\cdots ,N_{sw}$ (5)

Subject to:

$\frac{1}{{\mathbb{N}}_k}{\displaystyle {\sum }_{z=1}^{{\mathbb{N}}_k}G}\left(Y_{gw},{\bar{\varpi }}_{{\mathbb{N}}_k}\right)\le b,z=1,2,\cdots ,{\mathbb{N}}_k,k=1,2,\cdots ,{\mathbb{N}}_{sp}$ (6)

It should be noted that formulation (5) through (6) is deterministic optimization problem that can be solved by any appropriate standard core deterministic optimization solver. The conjunctive use optimization problems in the form of formulation (1) through (2) and formulation (5) through (6) are referred to as the true optimization problem and estimates optimization problem, respectively.

Note that in the optimization problem inequality (6), the term $\frac{1}{{\mathbb{N}}_k}$ defines the weight factor or probability associated with realizations $\varpi$. In this work, different number of realizations ${\mathbb{N}}_k$ (sample sizes) were considered for every sample optimization sub-problem. The performance objective function measured is the expected total optimal conjunctive use withdrawal rate. Decision variables measured are groundwater and surface water withdrawal rates (i.e., positive real values Y, such that $Y\in {ℝ}^N$, in which N is the total number of surface-subsurface abstraction points). The solution set $Y^{\text{*}}$ is assumed to be closed and bounded, hence the optimization sub-problems generated have finite set of feasible solutions.

2.2. ROA Approach

Retrospective Optimization Approximation (ROA) approach optimizes through sampling technique process. In which, instead of considering all samples (realizations) at each iteration of optimization, ROA procedure generates a sequence of sample sets of optimization sub-problems and subsequently solves and evaluates the optimization sub-problems generated in an increasing number of sample size (realizations) and updating initial guess solution through a “warm start” technique. In which, the initial guess is updated where the current sub-problem initial guess (solution) is just the solution from the previous sub-problem solved. The number of sample size (realizations) is increased sequentially from sub-problem to sub-problem resolved. Thus, the “warm start” technique applied ensures an increasing solution accuracy and a decreasing error of tolerance. Thus, ROA technique guarantees solution accuracy and speedy convergence to the true optimization problem solution. This is because in early iterations, ROA optimizations do not require much computational time as the number of realizations (sample sizes) is small. Similarly, at later iterations, the computational effort is inexpensive because the initial solutions are closer to the optimal solution of the true optimization problem. Hence, this is advantageous because the overall iterations required by the core optimizer is relatively fewer, and hence, computational savings can be attained compared to a direct optimization approach which considers all realizations at every iteration. The stochastic conjunctive use management optimization problem was solved and evaluated through ROA approach based on the following procedure illustrated through flow chart as shown in Figure 1.

2.3. Application of ROA Approach to a Real-World Regional Aquifer Water System

2.3.1. Study Area

The Retrospective Optimization Approximation (ROA) procedure developed

Figure 1. Flow chart for the ROA procedure.

was demonstrated the application to a real-world regional aquifer conjunctive water use system in the Great Letaba River (GLR) catchment (Figure 2(a)). The GLR catchment forms part of primary drainage region B of the Olifants River basin Water Management Area (WMA) which is located in the Mopani District of the Limpopo Province, South Africa. The GLR is one of the major tributaries of the Olifants River with drainage area of approximately 4952 km² [31].

2.3.2. Hydrology and Hydrogeological Settings

The study area is mainly drained by the Great Letaba River (also known as Groot Letaba River) and its tributaries. From the confluence of the Great Letaba and Klein Letaba Rivers, the Great Letaba River flows to the eastward through the Kruger National Park (KNP) until it joins with the Olifants River near the border to Mozambique. The study area regional water supply scheme uses water from the Great Letaba River and its tributaries to supply water to a number of towns including Polokwane (Pietersburg), Tzaneen, Duiwelskloof, Haenertsburg and to various villages. Also, extensive irrigation activities occurring within the catchment is supplied with water from this water supply system [32]. In the catchment, a total of 20 abstraction points is identified to be active. Figure 2(b) shows the Great Letaba River Basin quaternary catchments of the study area.

Groundwater pumping wells are also used to supply water to the study area and are spatially distributed throughout the catchment except in the quaternary catchment (QC) B81A where there is no pumping well (see Figure 2(b) and Figure 2(c)). Aquifers located in the catchment are predominantly secondary with exception of alluvial deposits along the main river system. Intergranular aquifers (ranging from unconsolidated to semi-consolidated materials with primary porosity occur in the study area. Figure 2(c) shows the pumping and observation wells distribution as well as the hydrogeological regions of the aquifer system.

In the study area, there exist 809 boreholes which are operational and another 294 boreholes are blocked or abandoned due to various reasons such as poor water quality. A total of eight (8) groundwater monitoring (control points) wells have been installed by the Government for groundwater level measurements (see Figure 2(c) blue dots). The existence of hot springs within the regional aquifer system suggests that the aquifer is confined except along the major rivers where

(a)(b)(c)

Figure 2. (a) Study area location map of Great Letaba River Basin; (b) Great Letaba River Basin Quaternary Catchments of the Study Area; (c) Great Letaba River Basin Pumping Wells, Observation Wells and Hydrogeological Regions of the Study Area Aquifer System.

localized alluvium aquifers occur including the eastern end of the catchment within the Kruger National Park (KNP) where unconfined aquifers occur. The aquifer water system is categorized by low magnitude of transmissivity values, ranging from a minimum of 7 m²/day to a maximum of 31 m²/day [33]. Recharge is mainly realized within the high elevation areas of the Great Escarpment (the Drakensburg Escarpment) mountains where rainfall is high above 1000 mm/annum and in the alluvium aquifers along major rivers. Total net recharge is about 126 mm/annum (DWAF, 2003) [32].

2.3.3. Conjunctive Water Use Conceptual and Simulation Models

Conceptual model was developed and used to build up conjunctive use simulation model. Figure 3 shows a schematized of the various components of the developed conjunctive water use conceptual model of the GLR catchment.

From Figure 3, DDA is Dams Direct Abstraction; RCDA is River Course Direct Abstraction; GWP is Groundwater Pumping; U/S and D/S are upstream and downstream, respectively. To delineate the study area regional aquifer water system from local aquifers, the whole study area was discretized into 150 × 100 cells, each in a model domain of a grid cell dimensional size of 1500 m × 1500 m. The area of interest was discretized through finite difference method. Figure 4 shows finite difference groundwater flow numerical simulation model.

Figure 3. Schematized conjunctive water use conceptual model.

Figure 4. Finite difference numerical simulation model of the study area.

It was realized that in the modelled domain, no-flow boundary appears at the northern and southern parts of the catchment, while constant flux boundary conditions exist for the two boundaries located at the north-west and south-west parts (these are areas bounded by the Great Escarpment (i.e., the Drakensburg Escarpment mountainous regions). Moreover, specific head boundary condition (also referred to as Dirichlet condition) for which head is considered to be independent of time was specified for the surface water bodies (i.e., rivers and Dams) and eastern boundary of the modelled domain. Other potential surface water bodies (such as wetlands and springs) were also included in the modelled area. In this study, steady-state simulations have been considered and hence, initial condition is not a major concern. This is because in any groundwater steady-state numerical simulations, the focus is to determine aquifer drawdown in response to externally imposed excitations/stresses such as groundwater pumping whereby relative heads as measured with respects to drawdown responses are of great importance rather than absolute head values. Due to limited availability of water demands data, the existing (current) water pumping (withdrawal) rates were considered as the minimum pumping (withdrawal) rate limits (i.e., lower bounds pumping (withdrawal) rates) which were required to satisfy the minimum water requirements of the competing water users. In this study, stochastic simulations were carried out to determine aquifer water system responses (drawdowns) due to unit withdrawal rate at each pumping well location for every set of realization of hydraulic conductivity field values realized. These aquifer water system responses (drawdowns) were assembled as response matrix and incorporated in the optimization model as hydraulic head constrain coefficients for the Retrospective Optimization Approximation (ROA) simulation-optimization problem analysis.

2.3.4. Model Data Input

Data on boreholes and surface water systems were found from GRIP Database [34] and the Department of Water and Sanitation of South Africa, respectively. Table 1 presents quaternary catchments names, named variables (combined pumping wells (CPW) withdrawal rates) and surface water diversion named variable (combined surface water diversion (CSWD) withdrawal rate) with their corresponding number of abstraction points, saturated mean aquifer thicknesses, and quaternary catchment polygon surface areas, which were used for the Retrospective Optimization Approximation (ROA) simulation-optimization problem analysis.

A total of 515 pumping wells with known capacity (average daily abstraction capacity ranges 62 - 70 m³/day) and 20 surface water abstraction points (i.e.,

Table 1. Summary characteristics of quaternary catchments (QC) and abstraction points used in study.

from river course and dams’ abstraction points with average daily abstraction capacity ranges 500 - 14,500 m³/day) were used. Hence, a total of 535 abstraction points (i.e., decision variables excluding slack variables) were used in simulation-optimization problem analysis. Table 2 presents model input parameter values of the river and aquifer systems properties.

2.3.5. Conjunctive Water Use Management under Hydrogeological Uncertainty

The overall objective was maximization of total conjunctive use withdrawal rates. Because of the water shortages of the study area for both irrigation and domestic use, surplus water from any water storage source can be used for irrigation and/or domestic. Hence, this objective seeks to maximize the amount of total conjunctive water production that can sustainably be withdrawn from both aquifer and surface water storage systems subjected to a number of constraints. In this work, the following objective and constraints were considered:

1) Conjunctive Use Objective Function

The objective function was formulated as follows:

$\text{Maximize}\left[F\left(Y\right)={\displaystyle {\sum }_{j=1}^{N_{pw}}{Y}_{gw,j}}+{\displaystyle {\sum }_d^{N_{sd}}{Y}_{sw,d}}\right],j=1,\cdots ,N_{pw};d=1,\cdots ,N_{sd}$ (7)

where $F\left(Y\right)$ is the objective function which defines the total of groundwater and surface water withdrawal rates; $Y_{gw}$ and $Y_{sw}$ are the spatially distributed surface and subsurface water withdrawal rates, respectively; $N_{pw}$ and $N_{sd}$ are the total number of pumping wells and surface water diversion points, respectively.

2) Conjunctive Use Constraints

The following objective function constraints were considered:

a) Drawdown Constraint

Drawdown constraint is meant for protection of the aquifer and ecosystem to circumventing excessive drawdowns. These constraints were formulated as follows:

${\displaystyle {\sum }_{j=1}^{N_{pw}}{a}_{i,j}{Y}_{gw,j}}\le b_i;i=1,\cdots ,N_c$ (8)

Table 2. Parameter values of river—aquifer systems properties.

in which $a_{i,j}$ is the response at control point i due to a unit pumping rate from pumping well located at j; $b_i$ is the allowable drawdown (response) at control point i; and $N_c$ is the total number of monitoring wells control points. Since the hydraulic conductivity field parameter values are sought to be uncertain, so the drawdowns $a_{i,j}$ are ought to be dependent on realizations of the hydraulic conductivity field, $\varpi$. In this case, the assemblage of these random response coefficients values $a_{i,j,\varpi }$ (which can simply be denoted by $\bar{\varpi }$ ) defines a finite set of independent identically distributed (i.i.d.) samples of random response matrices denoted by $M_{\varpi }$ of $\mathbb{N}$ realizations as $\left\{M_{{\varpi }_{{\mathbb{N}}_k}}\right\}=\left\{M_{{\bar{\varpi }}_{{\mathbb{N}}_1}},\cdots ,M_{{\bar{\varpi }}_{{\mathbb{N}}_k}}\right\}$, for $k=1,2,\cdots ,{\mathbb{N}}_{sp}$. Therefore, inequality (8) can be transformed into the following form [2] [4] [5] :

${\mathbb{E}}_{\text{Ω}}\left[G\left(Y_{gw,j},M_{{\varpi }_{{\mathbb{N}}_k}}\right)\right]\le b_i,i=1,2,\cdots ,N_c;k=1,2,\cdots ,{\mathbb{N}}_{sp};\forall \varpi \in \Omega$ (9)

Thus, the stochastic optimization inequality constraint (9) can be approximated as follow:

$\frac{1}{{\mathbb{N}}_k}{\displaystyle {\sum }_{z=1}^{{\mathbb{N}}_k}G}\left(Y_{gw},M_{{\bar{\varpi }}_{{\mathbb{N}}_k}}\right)\le b_i,i=1,2,\cdots ,N_c;k=1,2,\cdots ,{\mathbb{N}}_{sp};z=1,2,\cdots ,{\mathbb{N}}_k$ (10)

in which $M_{{\bar{\varpi }}_{{\mathbb{N}}_k}}$ is the k^th sample optimization sub-problem constraint response matrix, and all other parameters are as previously defined.

b) Total Recharge Constraint

To protect the aquifer from excessive exploitation, the total amount of water pumped from the aquifer was controlled so as not to surpass the total natural recharge. The constraint was considered as:

${\displaystyle {\sum }_{j=1}^{N_{pw}}{Y}_{gw,j}}\le \text{TWR}$ (11)

in which TWR is the total recharge of the aquifer system.

c) Base Flow to River Constraint

In this constraint aquifer heads at grid cells near to river course were restricted not to fall below riverbed bottom elevations. This constraint was set to maintain base flow water to river and ensures downstream water requirements. The constraint was formulated as follows:

$T{H}_{aq,j}\ge TRIVBO{T}_r$ (12)

where $T{H}_{aq,j}$ is the water head/water table level at aquifer grid cell location j near to river course; and $TRIVBO{T}_r$ is the riverbed bottom elevation at river reach cell r.

d) Total Water Demand Constraint

The aquifer and river (surface water) were considered as the only sources of water supply. Hence, this means that the premeditated optimal conjunctive withdrawal water policy must gratify at least the minimum total water demand without affecting negatively other water sources. This constraint was considered as follows:

${\displaystyle {\sum }_{j=1}^{N_{pw}}{Y}_{gw,j}}+{\displaystyle {\sum }_d^{N_{sd}}{Y}_{sw,d}}\ge \text{TWR}$ (13)

in which TWR is the total water requirement; and all other parameters are as previously defined.

e) Surface Water and Groundwater Withdrawal Constraints

To ensure positive real values, the water abstraction (withdrawal) rates at every pumping well and surface water abstraction point, were controlled to values between some values minimum and maximum withdrawal rates. These constraints were considered as follows:

$Y_{gw,j}^{\min }\le Y_{gw,j}\le Y_{gw,j}^{\max };j=1,2,\cdots ,N_{pw}$ (14)

$Y_{sw,d}^{\min }\le Y_{sw,d}\le Y_{sw,d}^{\max };d=1,2,\cdots ,N_{sd}$ (15)

where $Y_{gw}^{\min }$, $Y_{gw}^{\max }$ and $Y_{sw}^{\min }$, $Y_{sw}^{\max }$ are the minimum and maximum permissible groundwater and surface water withdrawal rates, respectively.

2.3.6. Statement of the Conjunctive Water Use Management Problem

The stochastic conjunctive water uses optimization problem solved, therefore, was formulated as follows:

$\text{Maximize}\left[F\left(Y\right)={\displaystyle {\sum }_{j=1}^{N_{pw}}{Y}_{gw,j}}+{\displaystyle {\sum }_d^{N_{sd}}{Y}_{sw,d}}\right],j=1,2,\cdots ,N_{pw};d=1,2,\cdots ,N_{sd}$ (16)

Subject to:

${\mathbb{E}}_{\Omega }\left[G\left(Y_{gw,j},M_{{\varpi }_{{\mathbb{N}}_k}}\right)\right]\le b_i,i=1,2,\cdots ,N_c;k=1,2,\cdots ,{\mathbb{N}}_{sp};\forall \varpi \in \Omega$ (17)

${\displaystyle {\sum }_{j=1}^{N_{pw}}{Y}_{gw,j}}\le \text{TWR}$ (18)

$T{H}_{aq,j}\ge TRIVBO{T}_r$ (19)

${\displaystyle {\sum }_{j=1}^{N_{pw}}{Y}_{gw,j}}+{\displaystyle {\sum }_d^{N_{sd}}{Y}_{sw,d}}\ge \text{TWR}$ (20)

$Y_{gw,j}^{\min }\le Y_{gw,j}\le Y_{gw,j}^{\max };j=1,2,\cdots ,N_{pw}$ (21)

$Y_{sw,d}^{\min }\le Y_{sw,d}\le Y_{sw,d}^{\max };d=1,2,\cdots ,N_{sd}$ (22)

The optimization problem (formulations (16) through (22)) is stochastic optimization problem. This is because the optimal solution realized depends on the outcomes of realizations of hydraulic conductivity field $\varpi$. The stochastic optimization problem solutions can be estimated through Retrospective Optimization Approximation (ROA) approach hence, the conjunctive water use sample optimization sub-problems solved, were considered as follows:

$\text{Maximize}\left[F\left(Y\right)={\displaystyle {\sum }_{j=1}^{N_{pw}}{Y}_{gw,j}}+{\displaystyle {\sum }_d^{N_{sd}}{Y}_{sw,d}}\right],j=1,\cdots ,N_{pw};d=1,\cdots ,N_{sd}$ (23)

Subject to:

$\frac{1}{{\mathbb{N}}_k}{\displaystyle {\sum }_{z=1}^{{\mathbb{N}}_k}G}\left(Y_{gw},M_{{\bar{\varpi }}_{{\mathbb{N}}_k}}\right)\le b_i,i=1,2,\cdots ,N_c;k=1,2,\cdots ,{\mathbb{N}}_{sp};z=1,2,\cdots ,{\mathbb{N}}_k$ (24)

${\displaystyle {\sum }_{j=1}^{N_{pw}}{Y}_{gw,j}}\le \text{TWR}$ (25)

$T{H}_{aq,j}\ge TRIVBO{T}_r$ (26)

${\displaystyle {\sum }_{j=1}^{N_{pw}}{Y}_{gw,j}}+{\displaystyle {\sum }_d^{N_{sd}}{Y}_{sw,d}}\ge \text{TWR}$ (27)

$Y_{gw,j}^{\min }\le Y_{gw,j}\le Y_{gw,j}^{\max };j=1,\cdots ,N_{pw}$ (28)

$Y_{sw,d}^{\min }\le Y_{sw,d}\le Y_{sw,d}^{\max };d=1,\cdots ,N_{sd}$ (29)

Note that the inequality (24) now becomes deterministic thus, formulation (23) through (29) is deterministic optimization problem which can be solved by any appropriate core deterministic optimization solver. This is the major advantage of ROA approach. The sample optimization sub-problems generated were solved and evaluated by using the ROA approach.

2.3.7. Formulation of Sample Path Optimization Sub-Problems

A correlation length of 100,000 m by 50,000 m in 2D-Dimensional, x, y-plane, was found to be sufficient to characterize input parameter uncertainties. In total, 500 realizations of uncertain aquifer hydraulic conductivity fields were produced. In the aquifer system, eight (8) groundwater control points were identified active. The water heads at these control points were restricted not to fall below certain prescribed maximum allowable limiting values. Twenty (20) abstractions points from river system network were identified active and the surface water levels (stages) at these abstraction points were constrained not to fall below riverbed bottom elevations. Collection of the aquifer system responses due to the 500 realizations resulted in 4020 constrain rows (i.e., observations rows). Hence, a response matrix of the size 4020 by 535 was produced and used to generate ten (10) sample optimization sub-problems of different number of sample rows (observation rows) for the optimization model. Sample sizes were designed heuristically. Table 3 presents the formulation of sample optimization sub-problems.

From Table 3, a sequence of 20, 30, 40, 60, 80, 100, 150, 200, 250, 500 realizations of hydraulic conductivity field realized, generated a sequence of 180, 260, 340, 500, 660, 820, 1220, 1620, 2020, 4020 of constraints, respectively (excluding total recharge, total demand and surface water base flow constrains). This sequence of constraints produced the corresponding ten (10) sample optimization sub-problems in a sequence of increasing number of observation rows (including total recharge, total demand, and surface water base flow constrains) of 183 × 535, 263 × 535, 343 × 535, 503 × 535, 663 × 535, 823 × 535, 1223 × 535, 1623 × 535, 2023 × 535, and 4023 × 535 (excluding lower and upper bounds, and nonnegative bound constraints). The last sample optimization sub-problem (i.e., SOSP10) is recognized to be the true conjunctive use optimization problem.

3. Results and Discussion

In total, 535 abstraction points (i.e., decision variables excluding slack variables) from groundwater and surface water sources were considered in simulation-optimization problem analysis. Currently, surface water source produces 87 percent of the total un-optimized conjunctive water use withdrawal rate and the remaining 13 percent is produced by groundwater source. Figure 5 shows overall percentages contributions of the existing (un-optimized) conjunctive use of surface water and groundwater sources in a pie chart view.

Table 3. Formulation of sample optimization sub-problems.

*Decision variables for conjunctive use withdrawal rates in which, the first 515 decision variables/columns represent groundwater withdrawal rates and the remaining 20 decision variables/columns represent surface water withdrawal rates.

Figure 5. Existing un-optimized overall percentage contribution of surface and subsurface water sources.

In this work, MODFLOW [35] simulation model and Active-Set core optimizer (Sequential Quadratic Programming (SQP) algorithm, implemented under MATLAB 2014a environment) codes were used for simulation-optimization problem analysis. The k-means clustering sampling technique was used for the realizations mapping. The objective of the optimization was to determine optimal withdrawal rates that maximize the expected total conjunctive use of surface and subsurface water withdrawal rates. To present the results in different graphical formats, the optimal conjunctive use solutions (i.e., the optimal groundwater withdrawal rates) obtained were sorted out according to their pumping wells location in the model domain and their associated quaternary catchments, and then for each quaternary catchment polygon a geometrical mean of the optimal groundwater withdrawal rates were computed. The optimal surface water mean withdrawal rates were also computed and set to a multiple of factor 0.1 so as to reduce the high differences between the surface and subsurface water withdrawal rates. Table 4 summarizes the Retrospective Optimization Approximation (ROA) conjunctive use sample path optimization sub-problems optimal mean withdrawal rates (in m³/day).

From Table 4, it can be seen that the expected total optimal conjunctive use objective function average (mean) values range from a minimum of 14835.72 m³/day to a maximum of 43453.36 m³/day. Moreover, the percentage error of tolerance decreases as the number of sample size (realizations) increases. The errors of tolerance range from 65.86 percent (i.e., at 20 realizations) to 0.00 percent (i.e., at 500 realizations). Note that the optimal conjunctive use average withdrawal volume rate solution values corresponding to the 500 hydraulic conductivity realizations varies from one another (see Table 4). This is because the optimal withdrawal volume rates designed based on Retrospective Optimization Approximation (ROA) approach depend on the outcomes of the uncertainty realizations. Figure 6 shows the conjunctive water use simulation-optimization sample optimization sub-problems optimal solutions.

Also, Figure 6 indicates that the sample optimization sub-problems solutions converge to the true conjunctive use optimization problem (i.e., SOSP10) as the number of realizations (sample size) increases. It can also be observed from

Table 4. ROA conjunctive use sample optimization sub-problems optimal mean withdrawal rates solutions (in m³/day).

Figure 6. Conjunctive water use simulation-optimization sample optimization sub-problems optimal solutions.

Figure 6 that after 30 realizations (sample size) the sample optimization sub-problems solutions are very close to true optimization problem solution. This is because in statistics, sample sizes below 30 are considered as small, hence, beyond which the convergence gains momentum and therefore, through Retrospective Optimization Approximation (ROA) approach it converges rapidly because the optimization process takes the advantage of “warm start” technique whereby initial guess solution of subsequent sub-problem is the solution of the current sub-problem. From Figure 6, it can also be seen that the ninth sample optimization sub-problem (SOSP9) optimal solution is almost equal to the original (true) optimization problem (SOSP10) solution. Thus, the tenth sample optimization problem (SOSP10) (in which all the 500 realizations generated were considered), converged with few iterations. This is because the initial guess solution of problem (SOSP10) is the solution of optimization sub-problem SOSP9 which is probably nearly equal to the true optimization problem (SOSP10) optimal solution. Figure 7(a) presents histogram diagram of the optimal conjunctive use withdrawal rates solution of the true optimization problem (SOSP10).

In Figure 6 and Figure 7(a), it can be observed that the highest groundwater withdrawal rates (>1000 m³/day) occurred within the quaternary catchments B81G, B81H, B81J, B81D with combined pumping wells CPW1, CPW2, CPW7, CPW8 (see Table 1) while the lowest volume rate (<250 m³/day) occurred at quaternary catchment B81C with combined pumping well CPW4. High differences in groundwater pumping rates are due to differences in aquifer hydraulic conductivity values and boundary conditions of the model domain. This is because quaternary catchments with combined pumping wells CPW1, CPW2, CPW7 and CPW8 are characterized by hydraulic conductivity values of relatively high magnitude values while the quaternary catchment with combined pumping well CPW4 falls within a relatively low magnitude of hydraulic conductivity value.

To evaluate the performance of the Retrospective Optimization Approximation (ROA) conjunctive use management model, the sample optimization sub-problems were resolved with changed initial solutions for three runs, which resulted

(a)(b)(c)

Figure 7. (a) Optimal conjunctive water use withdrawal rates of surface water and groundwater sources; (b) Performance of the ROA conjunctive use management model; (c) ROA optimal percentages of contribution of surface and subsurface water sources to the expected total optimal conjunctive use.

in different optimal solutions and, hence, different expected total conjunctive use objective function values. The expected total optimal conjunctive use withdrawal rates obtained were then evaluated over 500 realizations for three runs and averaged. Figure 7(b) shows the performance of the ROA conjunctive use management model with cluster sampling evaluated over 500 realizations.

From Figure 7(b), the Retrospective Optimization Approximation (ROA) expected total optimal conjunctive use withdrawal volume rate converged to its maximum mean withdrawal rate value of 4.35 × 104 m³/day within 7 to 8 iterations of which, surface water source contributed 36356.40 m³/day that is about 58 percent of the expected total optimal conjunctive use withdrawal rate and the remaining 7143.60 m³/day (i.e., 42 percent) was contributed by the groundwater source. Figure 7(c) presents ROA optimal percentages of contribution of surface and subsurface water sources to the expected total optimal conjunctive use withdrawal rate in a pie chart view.

Also, Figure 7(c) indicates that the overall percentages of contribution of surface and subsurface water sources to the total volume water withdrawal rate requirement were about 58 percent and 42 percent, respectively. This implies that if Retrospective Optimization Approximation (ROA) technique would be adopted for used in the study area, groundwater source can optimally (sustainably) be able to contribute up to 42 percent of the expected total optimal withdrawal rate. This is an increase of about 29 percent from the existing un-optimized groundwater source percentage (i.e., 13 percent) of contribution while the surface water source contribution would be reduced by 29 percent from the existing un-optimized surface water source percentage (i.e., 87 percent) of contribution.

4. Conclusion

Retrospective Optimization Approximation (ROA) methodology was applied to solve a real-world conjunctive water use stochastic simulation-optimization problem. ROA procedure solves and evaluates a sequence of optimization sub-problems in an increasing number of realizations. Results indicated that if ROA approach would be adopted for use in the study area, the existing groundwater source percentage of contribution to the total water requirement could increase up to 29 percent (i.e., from 13 to 42 percent) while the surface water source percentage of contribution could be reduced by 29 percent (i.e., from 87 to 58 percent). This implies that currently the river basin surface water storage system is water stressed and hence, reduction of 29 percent from the existing surface water source contribution could reduce stress on surface water storage system. This further demonstrates that ROA conjunctive water use management technique has potential to ensure sustainability of the limited water resources of river basins. Moreover, it has been revealed that through ROA method the expected total optimal conjunctive use objective function value converged to its maximum value within a relatively few iterations (6 to 8 iterations) in about 2.30 Hrs computational time. In conclusion, results demonstrated that the ROA-stochastic simulation-optimization approach is a promising technique for use in managing regional aquifers conjunctive water use under hydrogeological uncertainty conditions.

Acknowledgements

The financial assistance of the Tshwane University of Technology through the Faculty of Engineering and the Built Environment is gratefully acknowledged. The author also thanks the TUT SNOWS Project team for their support. The author further extends his gratitude to the Department of Water and Sanitation of South Africa for assisting with data.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Gyamfi, C., Ndambuki, J.M. and Salim, R.W. (2016) Hydrological Responses to Land Use/Cover Changes in the Olifants Basin, South Africa. Water, 8, 588. https://doi.org/10.3390/w8120588
[2]	Kifanyi, G.E. (2018) Optimizing the Conjunctive Use of Surface Water and Groundwater in Water Stressed River Basins: Case of Olifants River Basin, South Africa. PhD Thesis, Tshwane University of Technology, Tshwane.
[3]	Mahjoub, H., Mohammad, M. and Parsinejad, M. (2011) Conjunctive Use Modeling of Groundwater and Surface Water. Journal of Water Resource and Protection, 3, 726-734. https://doi.org/10.4236/jwarp.2011.310083
[4]	Kifanyi, G.E., Ndambuki, J.M. and Odai, S.N. (2016) A Quantitative Groundwater Resource Management under Uncertainty Using a Retrospective Optimization Framework. Sustainability, 9, 2.
[5]	Ndambuki, J.M., Kifanyi, G.E., Odai, S.N. and Gyamfi, C. (2018) Conjunctive Management of Surface and Groundwater Resources under Uncertainty: A Retrospective Optimization Approach. International Journal of Energy and Environmental Engineering, 12.
[6]	Regulwar, D.G. and Pradhan, V. S. (2013) Irrigation Planning with Conjunctive Use of Surface and Groundwater Using Fuzzy Resources. Journal of Water Resource and Protection, 2013, 816-822. https://doi.org/10.4236/jwarp.2013.58082
[7]	Safavi, H.R. and Bahreini, G.R. (2009) Conjunctive Simulation of Surface Water and Groundwater Resources under Uncertainty. Iranian Journal of Science and Technology, Transaction B: Engineering, 33, 79-94.
[8]	Ndambuki, J.M. (2001) Multi-Objective Groundwater Quantity Management: A Stochastic Approach. DUP Science, Delft University Press, Delft.
[9]	Schoups, G., Addams, C.L., Minjares, J.L. and Gorelick, S.M. (2006) Sustainable Conjunctive Water Management in Irrigated Agriculture: Model Formulation and Application to the Yaqui Valley. Water Resources Research, 42, W10417. https://doi.org/10.1029/2006WR004922
[10]	Marques, G. F., Lund, J. R. and Howitt, R. E. (2010) Modelling Conjunctive Use Operations and Farm Decisions with Two-Stage Stochastic Quadratic Programming. Journal of Water Resources Planning and Management, 136, 386-394. https://doi.org/10.1061/(ASCE)WR.1943-5452.0000045
[11]	Schoups, G., Addams, C.L., Minjares, J.L. and Gorelick, S.M. (2006) Reliable Conjunctive Use Rules for Sustainable Irrigated Agriculture and Reservoir Spill Control. Water Resources Research, 42, W12406. https://doi.org/10.1029/2006WR005007
[12]	Kentel, E. and Aral, M.M. (2007) Fuzzy Multi-Objective Decision Making Approach for Groundwater Resources Management. Journal of Hydraulic Engineering, 12, 206-217. https://doi.org/10.1061/(ASCE)1084-0699(2007)12:2(206)
[13]	Young, R.A. and Bredehoeft, J.D. (1972) Digital Computer Simulation for Solving Management Problems of Conjunctive Groundwater and Surface Water Systems. Water Resources Research, 8, 533-556. https://doi.org/10.1029/WR008i003p00533
[14]	Provencher, B. and Burt, O. (1994) Approximating the Optimal Groundwater Pumping Policy in a Multi-Aquifer Stochastic Conjunctive Use Setting. Water Resources Research, 30, 833-843. https://doi.org/10.1029/93WR02683
[15]	Gilling, D., McCarl, A. and Boadu, F. (2001) An Economic, Hydrologic and Environmental Assessment of Water Management Alternative Plans for South Central Texas Region. Journal of Agricultural and Applied Economics, 33, 59-68. https://doi.org/10.1017/S1074070800020782
[16]	Wang, H., Ciaurri, D.E., Durlofsky, L.J. and Cominelli, A. (2012) Optimal Well Placement under Uncertainty Using a Retrospective Optimization Framework. SPE Journal, 17, 112-121. https://doi.org/10.2118/141950-PA
[17]	Gorelick, S.M. (1983) A Review of Distributed Parameter Groundwater Management Modelling Methods. Water Resources Research, 19, 305-319. https://doi.org/10.1029/WR019i002p00305
[18]	Aguado, E. and Remson, I. (1974) Groundwater Hydraulics in Aquifer Management. Journal of the Hydraulics Division, 100, 103-118. https://doi.org/10.1061/JYCEAJ.0003848
[19]	Peralta R.C., Azarmnia H. and Takahashi, S. (1991) Embedding and Response Matrix Techniques for Maximizing Steady-State Groundwater Extraction: Computational Comparison. Groundwater, 29, 357-364. https://doi.org/10.1111/j.1745-6584.1991.tb00526.x
[20]	Gorelick, S.M.A. (1982) Model for Managing Sources of Groundwater Pollution. Water Resources Research, 18, 773-781. https://doi.org/10.1029/WR018i004p00773
[21]	Ndambuki, J.M., Otieno, F.A.O., Stroet, C.B.M. and Veling, E.J.M. (2000) Groundwater Management under Uncertainty: A Multi-Objective Approach. Water SA, 26, 35-42.
[22]	Wang, H.G. and Schmeiser, B.W. (2008) Discrete Stochastic Optimization Using Linear Interpolation. In: Proceedings Winter Simulation Conference, IEEE Press, Piscataway, 502-508. https://doi.org/10.1109/WSC.2008.4736106
[23]	Chen, H. and Schmeiser, B.W. (2001) Stochastic Root Finding via Retrospective Approximation. IIE Transactions, 33, 259-275. https://doi.org/10.1080/07408170108936827
[24]	Chen, Y., Oliver, D.S. and Zhang, D. (2009) Efficient Ensemble-Based Closed-Loop Production Optimization. SPE Journal, 14, 634-645. https://doi.org/10.2118/112873-PA
[25]	Wang, C., Li, G. and Reynolds, A.C. (2009) Production Optimization in Closed-Loop Reservoir Management. SPE Journal, 14, 506-523. https://doi.org/10.2118/109805-PA
[26]	Güyagüler, B. and Horne, R.N. (2004) Uncertainty Assessment of Well Placement Optimization. SPE Reservoir Evaluation & Engineering, 7, 24-32. https://doi.org/10.2118/87663-PA
[27]	Onwunalu, J.E. and Durlofsky, L.J. (2010) Application of a Particle Swarm Optimization Algorithm for Determining Optimum Well Location and Type. Computers & Geosciences, 14, 183-198. https://doi.org/10.1007/s10596-009-9142-1
[28]	Özdogan, U. and Horne, R.N. (2006) Optimization of Well Placement under Time-Dependent Uncertainty. SPE Reservoir Evaluation & Engineering, 9, 135-145. https://doi.org/10.2118/90091-PA
[29]	Yeten, B., Durlofsky, L.J. and Aziz, K. (2003) Optimization of Nonconventional Well Type, Location, and Trajectory. SPE Journal, 8, 200-210. https://doi.org/10.2118/86880-PA
[30]	Kifanyi, G.E., Ndambuki, J.M., Odai, S.N. and Gyamfi, C. (2019) Quantitative Management of Groundwater Resources in Regional Aquifers under Uncertainty: A Retrospective Optimization Approach. Groundwater for Sustainable Development, 8, 530-540. https://doi.org/10.1016/j.gsd.2019.02.005
[31]	DWAF (2006) Letaba Catchment Reserve Determination Study: Groundwater Report. Department of Water Affairs and Forestry, South Africa. Final Report February 2006, WSM Civil Engineers, Hydrogeologists and Project Managers (Pty) Ltd., Report No. RDM/B800/02/CON/ COMP/0504.
[32]	DWAF (2003) Luvuvhu/Letaba Water Management Area: Water Resources Situation Assessment. Department of Water Affairs and Forestry, South Africa. Report No. P02000/00/0301.
[33]	Holland, M. (2011) Hydrological Characterization of Crystalline Basement Aquifers within the Limpopo Province, South Africa. PhD Thesis, University of Pretoria, Pretoria.
[34]	GRIP-Limpopo (2013) GRIP Database for Borehole Point Data for the Limpopo Province. http://griplimpopo.co.za
[35]	Harbaugh, A.W., Banta, E.R., Hill, M.C. and McDonald, M.G. (2000) MODFLOW-2000, the U.S. Geological Survey Modular Ground-Water Model—User Guide to Modularization Concepts and the Ground-Water Flow Process. OFR, 2000-92. https://doi.org/10.3133/ofr200092