Introduction

Arid and semiarid regions, also referred to as drylands, are areas characterised by low rainfall and high temperatures, which translate to high evapotranspiration rates, low groundwater recharge and usually low water availability (Gee and Hillel 1988; Maliva and Missimer 2012; McEwan 2006; Missimer et al. 2012). These circumstances present ecosystems and the human population in these regions with water scarcity. Climate change is expected to exert further pressure on water resources in drylands and to exacerbate the existing issues. According to the UN’s 5th Intergovernmental Panel on Climate Change assessment report (IPCC 2014), mid-latitudes and dry regions will likely receive less precipitation, while the average worldwide temperature is expected to increase at least 1.5 °C by the end of the century (relative to the years 1850–1900). In this context, water management solutions are of paramount importance in drylands if they are to cope with the natural climatic conditions and adapt to climate change.

A good example of such drylands can be found on the Arabian Peninsula, namely in the Kingdom of Saudi Arabia (KSA). Rainfall in this country occurs chiefly as sporadic events of high intensity and low predictability. When these events take place, stormwater converges to the channels of ephemeral rivers (wadis) and often develops into flash floods (Gee and Hillel 1988; Maliva and Missimer 2012; Sen et al. 2011). The renewable groundwater resources in the KSA are mostly limited to unconfined to semiconfined aquifers which underlie wadis (Lopez et al. 2014; Maliva and Missimer 2012; Missimer et al. 2015, 2012). These are the so-called wadi aquifers. Recharge of surface water into wadi aquifers is not high and often ranges between 3 to 11% of the annual precipitation (Missimer et al. 2012).

The low recharge rates are often increased through the use of stormwater harvesting (SWH) systems and managed aquifer recharge (MAR). Stormwater harvesting is a term that refers to a set of solutions aimed to cope with water shortage by storing and utilising excess stormwater. MAR involves the intentional recharge of water through several techniques (Dillon et al. 2009). Representative examples of SWH and MAR in drylands can be found in the Kingdom of Saudi Arabia (Alataway and El Alfy 2019; Missimer et al. 2015), where the systems often rely on a combination of dams as the main storage and additional infiltration constructions as groundwater-recharge facilities (Chowdhury and Al-Zahrani 2015; Missimer et al. 2012; Sen et al. 2011). The dams are frequently located in the channels of wide-spread wadis to capture runoff from rainfall events for two primary purposes, namely flood prevention and groundwater recharge (Kalwa 2013; Lopez et al. 2015; Missimer et al. 2012; Sen et al. 2011). This sort of stormwater recharge system is known as a wadi dam (Missimer et al. 2015) or check dam (Dillon et al. 2019; Djuma et al. 2017) and is considered to be an in-channel modified MAR scheme (Standen et al. 2020).

When the infiltration rates in wadi reservoirs are low, water remains ponded for long periods, and the SWH system fails to satisfy its purposes. The higher the water level in the reservoir, the lower the flood buffering capacity; furthermore, the total groundwater recharge could be much lower than the stored volume due to evaporation. Losses to evaporation are as low as 5% of the water stored on an annual basis when the wadi reservoir is relatively new and as high as 80% in older structures (Lopez et al. 2014; Missimer et al. 2015). At a certain point, the excessive evaporation can also increase salinity (Missimer et al. 2015), and the prolongated surface storage of water offers potential for eutrophication and mosquito-breeding (DEC 2006; Philp et al. 2008).

The leading causes of low infiltration rates in wadi reservoirs are to a great extent the individual or combined effect of two phenomena: (1) the formation of a clogging layer (CL) on top of the existing soil profile due to the settlement of fine material in stagnant water (Alataway and El Alfy 2019; Kalwa 2013; Sen et al. 2011); (2) and the large thickness of the vadose zone (VZ), which can lead to low saturation and reduced hydraulic conductivity (Missimer et al. 2015; Walvoord et al. 2002).

There are several solutions available to address the low infiltration in wadi reservoirs. These solutions aim at different levels and components of the SWH system—for instance, some techniques deal with the CL, including removal of a few centimetres of topsoil, scratching of the surface, and ploughing (Al-Muttair et al. 1994; Maliva and Missimer 2012). Another set of solutions aims at infiltrating the stored water in areas with more favourable hydraulic conditions—for example, water can be released to a downstream portion of the wadi which is not capped by a clogging layer (Al-Muttair et al. 1994; Kalwa 2013; Lopez et al. 2014; Missimer et al. 2012) or it can be conducted in wells located downstream (Missimer et al. 2015). Furthermore, infiltration can also be augmented on-site by using structures such as vadose zone wells (also known as drywells; Alataway and El Alfy 2019; Kalwa et al. 2014) and trenches (Troeger and Wannous 2019). The decisive benefit of well and trench structures is the bridging of low permeability components such as the CL or an unsaturated zone. Depending on their construction, they present a low-cost and low-tech alternative that uses gravity-driven free drainage in comparison to well-based pumping solutions requiring continuous maintenance.

Several studies have assessed the use of trenches and wells to enhance surface (SF) infiltration, and the key findings are summarised in the Appendix. Sasidharan et al. (2018, 2019, 2020a) and Heilweil et al. (2015) focus on the numerical modelling of vadose zone wells and trenches, respectively. The research of Missimer et al. (2015) and Kalwa (2013) tries to bring a solution to the problem posed by evaporation in wadi dams. In the case of the former, they utilise groundwater modelling to explore the use of injection wells located downstream of the reservoir. Their models validate the use of the downstream scheme as it enhances infiltration. Kalwa (2013) uses mostly experimental data to compare three different solutions: infiltration through the bed of the downstream sediments, and infiltration directly within the reservoir with and without wells. He concludes that downstream infiltration and the wells are the best options to enhance wadi aquifer recharge than the sole reservoir. Händel et al. (2014, 2016) compare the infiltration from small-diameter wells with the infiltration from spreading basins and finds that wells yield more recharge than the spreading basins. Most recently, Sasidharan et al. (2020b) assess the effectiveness of infiltration basins and vadose zone wells in the context of managed aquifer recharge. In general, these articles point towards a better performance of vertical structures with high permeability (e.g. wells and trenches) in comparison with horizontal ones (e.g. spreading basins or direct surface infiltration). It is thought that the only work that addresses which of the infiltration solutions yields the highest infiltration in wadi reservoirs is that of Kalwa (2013). He states that his study leaves room for improvement due to the uncertainty in some parameter measurements and suggests that either wells or infiltration through the channel downstream of the reservoir bring about more groundwater recharge than surface infiltration within the reservoir.

Until now, and to the best knowledge of the authors, there is no systematic study that evaluates the performance of vadose zone wells and infiltration trenches to enhance on-site infiltration in wadi reservoirs, when compared with infiltration through the bed of the reservoir. Therefore, this study firstly evaluates the performance of wells and trenches in terms of recharge timing and quantity. Secondly, it explores the best infiltration solution and its optimal characteristics through a series of numerical exercises comprising the applicability of existing analytical equations and a sensitivity analysis primarily focused on different geometries of the wells and trenches. The results of this study are used to suggest potential water management solutions and optimisations of infiltration systems in drylands. In addition, this study provides a general framework and a basic understanding of vadose zone wells and infiltration trenches systems. The present study does not aim at characterising the optimal site conditions for the best performance of the infiltration solutions.

Methods

The well and the trench are modelled within a wadi reservoir, including the saturated and unsaturated zones. The models are in all cases conceptualised as 2D, corresponding to a cross-section or axisymmetric domain. A vertical cross-section is considered within the reservoir located far away from the borders of a channel or any structural barrier (e.g. a dam). Topography within the reservoir or in the underlying bedrock is not considered, and the cross-sections are oriented perpendicular to the hydraulic gradient of the wadi channel.

The scenarios are conceptualised using hydrogeological parameters representing arid and semiarid regions where a thick VZ or a CL pose a hindrance to infiltration. The selected values are representative for typical wadis in the KSA. As a significant variation, two different sets of models are evaluated:

  • Shallow water Table (WT) models. The thickness of the vadose zone is 10 m: this is a representative value for wadi aquifers in Saudi Arabia (Missimer et al. 2012).

  • Deep WT models. The thickness of the vadose zone is 50 m: this set of scenarios allows assessment of free drainage in comparison to the horizontal discharge of infiltrated water for shallow WT scenarios.

For both WT options, well and trench infiltration with simultaneous recharge through the surface of the reservoir (symbolised as “+ SF”) are investigated. These correspond to model setups 1 (shallow WT), 2 (shallow WT), 4 (deep WT), and 5 (deep WT) in Table 1. Additionally, these setups are compared to sole surface water (SF) infiltration in setups 3 (shallow WT) and 6 (deep WT). The scenarios in Table 1 will be referred to as the base scenarios, which are slightly modified in some sections.

Table 1 Base scenarios used in this study. Letter D is used as a prefix to denote a deep WT setting
Full size table

Two indicators are used to assess the effectiveness of the trench and the well to bypass the CL and the VZ: (1) the time at which recharge starts and, (2) the cumulative groundwater recharge over a year. The recharge rates are accounted for by a line located at the outer boundary of the domain which extends either vertically along the saturated thickness (shallow setups) or covers the lower horizontal domain boundary (deep setups). The quasi-steady-state recharge rates (QSSRR) of the scenarios are also estimated. They are used to carry out two additional analyses, namely the validation of the conceptual models and a sensitivity analysis.

The general structure of the different analyses carried out is the following—the scenarios in Table 1 are run, and the start of groundwater recharge is analysed (simulations of continuous infiltration). Subsequently, the complexity of the models is increased to transient conditions to allow for an assessment of the cumulative groundwater recharge (analysis of cumulative recharge). Next, a sensitivity analysis is carried out, exploring how the geometry of the well and the trench impact the QSSRR. This section is supplemented with a short analysis of the cost of infiltrating water in terms of the well and trench construction costs. Successively, simplified versions of the models presented in Table 1 are built, and their QSSRRs are compared against the infiltration rates obtained from analytical solutions. Note that, if quasi-steady-state conditions prevail, the recharge rate is equivalent to the infiltration rate. Finally, the analysed solutions are compared in terms of the cumulative recharge produced by the maximum number of wells and trenches fitted in a delimited area. Results and discussion are provided in the same section.

Conceptual background

The movement of water through the soil under unsaturated conditions is described by Richards’ equation (Richards 1931, Eq 1).

$$ \frac{\partial \theta (h)}{\partial t}=\frac{\partial }{\partial z}\left[K(h)\frac{\partial h}{\partial z}\right]+\frac{\partial K(h)}{\partial z}-S(h) $$
(1)

where h is the hydraulic head, K(h) the unsaturated hydraulic conductivity and S(h) a sink function that accounts for sinks as sources of water. This equation is a partial differential nonlinear equation which in this case represents one-dimensional (1D) vertical flow (z-direction) in a mixed form, i.e., there are two dependent variables, namely Ɵ and h. These two variables are linked through the water retention curve.

The water retention curve can be represented by mathematical expressions. One of the most frequently used ones was formulated by van Genuchten (1980; Eq. 2).

$$ \theta (h)=\frac{\theta_{\mathrm{s}}-{\theta}_{\mathrm{r}}}{{\left[1+{\left(-\alpha h\right)}^n\right]}^m}+{\theta}_{\mathrm{r}} $$
(2)

where α, n and m are fitting parameters (often referred to as the van Genuchten parameters), Ɵ(h) is the volumetric water content, Ɵs is the saturated volumetric water content, and Ɵr is the residual volumetric water content. It is commonly assumed that m = 1–1/n.

The parameter α is related to the inverse of the air entry value, and n is a measure of the pore size distribution. The residual volumetric water content is the water content measured at very negative pressure heads when the water retention curve approximates a flat region. Van Genuchten (1980) suggests that this value can be equitable to the wilting point.

There are several formulations to represent the unsaturated hydraulic conductivity, (K(h), as a function of the saturated hydraulic conductivity (Ks; e.g. Brooks and Corey 1964). In the numerical simulations conducted in this study, the equation by van Genuchten-Mualem is employed (Eq. 3).

$$ K(h)={K}_{\mathrm{s}}{S}_{\mathrm{e}}^I{\left[1-{\left(1-{S}_{\mathrm{e}}^{1/m}\right)}^m\right]}^2 $$
(3)

where K(h) is the unsaturated hydraulic conductivity, Ks the saturated hydraulic conductivity, Se the effective water saturation (Eq. 4), and I the pore-connectivity parameter. The latter parameter is on average 0.5 for diverse soil types (Mualem 1976).

$$ {S}_{\mathrm{e}}=\frac{\theta -{\theta}_{\mathrm{r}}}{\theta_{\mathrm{s}}-{\theta}_{\mathrm{r}}} $$
(4)

Conceptual models

The conceptual models employed in this study are identical to the ones used by Henao Casas (2019), who provides further details on the justification of the model parameters and characteristics.

Model dimensions

The trench is represented in a two-dimensional (2D) cross-section and the well in a 2D axisymmetric domain (pseudo-3D domain). Figure 1 shows the conceptual representation of the shallow and deep models. In the shallow scenarios, a WT is established from the bedrock and up to 10 m above the base of the models, while above, an unsaturated zone extends for 9 m. Finally, capping the VZ, there is a CL of 1 m of thickness. The dimensions of the parts of the system are based on the average values reported by Missimer et al. (2012) and the ranges provided by other authors for wadi deposit thickness (Al-Turbak et al. 1993; Masoud et al. 2019), depth to WT (Al-Shaibani 2008; Al-Turbak et al. 1993) and saturated thickness (Al-Turbak et al. 1993; Sorman et al. 1997).

Fig. 1
figure 1

Conceptual model of a the infiltration trench in a shallow water Table (WT) setting and b a vadose zone well in a deep WT setting. The region in colour represents the one simulated with HYDRUS 2D/3D. The black and white region is part of the assumption of symmetry concerning the axis of symmetry

Full size image

The trench has a depth of 3 m and a width of 1 m based on average values mentioned in the literature (Bouwer 2002; Schueler 1987). The vadose zone well is 5 m in depth and 1 m in diameter using the same criteria (Bouwer 2002; Liang et al. 2018; Sasidharan et al. 2019). In the shallow WT scenarios, the domain dimensions are set as 20 m in the vertical direction and 25 m in the horizontal one. In the deep WT scenarios, the domain corresponds to 49 m of VZ capped by 1 m of CL. In this case, the saturated thickness is excluded.

Boundary conditions

Boundary conditions (BCs) are shown in Fig. 1. In the shallow WT scenarios, a constant pressure head of 2.5 m is assigned to the horizontal upper infiltration BC, representing a constant water level within the reservoir. The infiltration from the trench and the well are represented using constant pressure heads in equilibrium with the upper BC. The vertical upper right boundary is set as the seepage face to allow the passage of excess flow, in accordance with Heilweil et al. (2015). The left BC is conceptually a symmetry boundary (and no flux technically), assuming in all cases axial symmetry. The horizontal lower BC is a no-flow boundary to represent the bedrock. The scenarios involving a deep WT are simulated using the same characteristics of the shallow WT models, except that the bottom of the former scenarios is limited by a water table (i.e. a constant pressure head of 0 m).

Hydraulic parameters

Material properties are summarised in Table 2. The material of the wadi aquifer is selected as ‘loamy sand’ from the soil catalogue of HYDRUS 2D/3D (Carsel and Parrish 1988; Šimunek and Sejna 2018a), whose hydraulic parameters are in the range of values found in the literature (Al-Shaibani 2008; Hussein et al. 1993; Kalwa 2013; Masoud et al. 2019; Missimer et al. 2015, Missimer et al. 2012; Rosas 2013, 2013; Sorman et al. 1997). The hydraulic conductivity is set as anisotropic and based on the ratio between the vertical and horizontal components reported by Missimer et al. (2012) and Sorman et al. (1997). The parameters of the ‘silty clay’ from the soil catalogue of HYDRUS 2D/3D are assigned to the CL in order to resemble the CL at the Al-Alb Dam (Kalwa 2013).

Table 2 Hydraulic parameters defined for the baseline scenario
Full size table

Numerical realisation

The software used to model the scenarios is HYDRUS 2D/3D version 3.01 (Šimunek and Sejna 2018a, b). The modelling area is discretised with an irregular mesh. The average element edge length is 0.7 m. Mesh refinements are applied along the well/trench contour (0.24 m), the WT (0.09 m), and the lower (0.13 m) and upper boundaries of the CL (0.37 m). The total number of nodes ranges between 5,960 and 6,950 for the set of models with a deep WT, and between 21,124 and 22,332 for those with a shallow WT (note: number of nodes is lower in deep WT scenarios as the saturated part of the simulation domain was left out).

Initial conditions

The bottom of the saturated thickness has a pressure head of 10 m. This pressure decreases linearly up to the WT, where the pressure head is 0 m. For the vadose zone, a pressure head of −2.3 m is assigned, and the scenarios are run until quasi-steady-state conditions are reached.

Simulations of continuous infiltration

The six scenarios presented in Table 1 are run and the time in which groundwater recharge starts is assessed (i.e. first arrival time of infiltration at the water table) and compared with the different means of infiltration. All scenarios are run until a quasi-steady state is attained. The respective recharge rates obtained (i.e. QSSRRs) are used as a reference in the sensitivity analysis.

Analysis of cumulative recharge

Events of periodic water infiltration into the wadi aquifer are simulated by considering the natural conditions in which wadi dams operate, and then translating them into the scenarios. The operation of the wadi dams implies input of water into the reservoir from rainfall and runoff, with simultaneous evaporation and infiltration into the wadi aquifer. The periodic simulations are configured with the same climatological conditions every year, and therefore they provide cumulative recharge representative of the long-term tendency of a year.

Yearly cumulative groundwater recharge from the well, the trench, and the SF are derived from the periodic simulations. Note that the recharge from these means of infiltration is in different domains, i.e., the well in volumes and the rest in areas. Therefore, the trench and the SF recharge are converted into volumes, and subsequently compared with the well recharge.

Two scenarios are considered to convert the trench recharge into volume: (1) the trench has the same storage volume as the well; (2) the trench has the same infiltration area (i.e. the area available for infiltrating water) as the well. Thus, the 2D cumulative recharge of the trench is multiplied by the lengths that assure the equivalence of volume (1.6 m) and area (2.7 m) to the well. In the case of the SF, the recharge is computed from a circular area which has the same radius of influence as the well (25 m).

Analytical solutions

Simplified versions of the models are set up to compare with analytical equations. These simplified models are equal to the ones described in the previous subsection, with a few exceptions. They do not have a clogging layer, nor anisotropy (the whole domain has a single hydraulic conductivity of 4 m/day). Besides, the entire length of the well and the trench up to the surface is set as a constant pressure head, and the hypothetical water level of the reservoir is established at the soil surface (pressure head of 0 m).

The simplified scenarios are run until quasi-steady-state conditions are attained. The QSSRRs of the well and the trench are compared with the infiltration rates calculated through the analytical equations.

Three analytical expressions are used to represent the well in a shallow WT setting. They correspond to the equation by Elrick et al. (1989; Eq. 5) and the method 1 of the USBR (1977; Eq. 6). The former equation is developed for the one-ponded height technique, which includes a term for the unsaturated water flow, while the latter is devised for a packer test in unconsolidated material at settings in which A ≥ 10r and A = H. The meaning and selected values of the variables in all the equations used in this section are presented in Table 3, and the schematic representation of the variables is depicted in Fig. 2.

$$ Q{\displaystyle \begin{array}{c}={K}_{\mathrm{s}}\left(\frac{2\ \uppi\ {H}^2}{C}+\uppi\ {r}^2+\frac{2\ \uppi\ H}{\upalpha\ C}\right)\\ {}C={H}^2\left(\frac{\frac{A}{H}{\sinh}^{-1}\left(\frac{A}{r}\right)}{A^2}-\frac{\sqrt{{\left(\frac{r}{H}\right)}^2+{\left(\frac{A}{H}\right)}^2}+\frac{r}{H}}{A^2}\right)\end{array}} $$
(5)
$$ {\displaystyle \begin{array}{c}Q=\frac{K_{\mathrm{s}}\ \left({C}_{\mathrm{s}}+4\right)\ r\ {T}_{\mathrm{u}}}{2}\\ {}{C}_{\mathrm{s}}=\frac{2\pi\ \left(\raisebox{1ex}{$A$}\!\left/ \!\raisebox{-1ex}{$r$}\right.\right)}{\ln\ \left(\raisebox{1ex}{$A$}\!\left/ \!\raisebox{-1ex}{$r$}\right.\right)}\end{array}} $$
(6)
Table 3 Variables of the analytical equations and the selected values
Full size table
Fig. 2
figure 2

Schematic showing the different variables used in the analytical equations and the sensitivity analysis

Full size image

The well in a deep WT setting is typified with the analytical equations of Glover (1953; Eq. 7) for a packer test and Bouwer’s equation (Bouwer 2002) for a vadose zone well (Eq. 8). Both equations are formulated for a deep WT setting in which A ≥ 10r. Glover’s equation is further restricted to configurations where H = A.

$$ {\displaystyle \begin{array}{c}Q={K}_{\mathrm{s}}\ r\ H\ {C}_{\mathrm{u}}\\ {}{C}_{\mathrm{u}}=\frac{2\uppi\ \left(\raisebox{1ex}{$H$}\!\left/ \!\raisebox{-1ex}{$r$}\right.\right)}{\sinh^{-1}\left(\raisebox{1ex}{$H$}\!\left/ \!\raisebox{-1ex}{$r$}\right.\right)-1}\end{array}} $$
(7)
$$ Q=\frac{2\ \uppi\ {K}_{\mathrm{s}}{H}^2}{\ln \left(\raisebox{1ex}{$2H$}\!\left/ \!\raisebox{-1ex}{$r$}\right.\right)-1} $$
(8)

For the trench, the equations used in all cases are the ones by Bouwer (2002; Eq. 9) and the empirical equation by Heilweil et al. (2015; Eq. 10). Note that trench infiltration according to Bouwer’s equation (Eq. 9) is 20% of the well infiltration calculated through Eq. 8.

$$ Q=\frac{0.4\ \uppi\ {K}_{\mathrm{s}}\ {H}^2}{\ln \left(4\ \frac{H}{W}\right)-1} $$
(9)
$$ q={K}_{\mathrm{s}}\ \left[2.79\ \ln (T)-5.47\right] $$
(10)

When choosing the values of the variables, the main factor to consider is the range of values in which the equations are valid and the resemblance to the base scenarios (Table 1).

Sensitivity analysis of well and trench dimensions and cost analysis

Sensitivity analysis

The dimensions of the trench and the well are modified with respect to the base scenarios. The modifications are performed in multiples of the original sizes; thus, the length and radius of the well are doubled and halved, and the trench depth and the width are modified in the same way. Subsequently, the modified scenarios are run until quasi-steady-state conditions are reached. The obtained QSSRRs are compared with those of the base scenarios.

Analysis of the construction costs

Costs of construction are allocated to the well and the trench base scenarios and all the variations considered in the sensitivity analysis (e.g. a well with twice the original diameter). The cost of construction is, in broad terms, a function of the well/trench volume. Based on the range of prices provided by Brown and Schueler (1997) for a vadose zone well and an infiltration trench (4–9$ per ft3 and a selected price of 7.5$ per ft3), volumes and construction costs are computed for each of the scenarios of the sensitivity analysis.

An indirect approach is used to estimate the cumulative recharge of the scenarios involving modifications in the original dimensions (note that the periodic simulations are strictly constructed for the base scenarios). The coefficient between the QSSRRs of the scenarios in the sensitivity analysis (numerator) and the base scenarios (denominator) is calculated. Subsequently, the obtained coefficient is multiplied by the cumulative recharge of the base scenarios and then divided by the respective construction cost to obtain the yearly cumulative recharge per dollar. Only the scenarios in a shallow WT setting are considered.

In addition to the scenarios of the sensitivity analysis, two more scenarios are evaluated. They are termed the “optimum” scenarios, and there is one for a well and one for a trench. They consist of the combination of vertical and horizontal dimensions, which ensure the highest QSSRRs per unit of volume. In the case of the well, this happens when the diameter is half of the original one, which is 0.5 m, and the length twice as long (10 m). The trench in the optimum scenario has both dimensions halved, so its depth is 1.5 m and the width 0.5 m.

Well and trench performance under areal constraints

The infiltration potential of wells and trenches on a fixed square area is assessed to mimic conditions with limited infiltration area. The maximum number of wells and trenches that can fit in a 1-km2 dam is calculated in accordance with the optimum distributions shown in Fig. 3. In these hypothetical distributions, the wells and trenches are positioned in such a way that their areas of influence do not overlap, enabling the arithmetic addition of the individual cumulative recharge. For the trench, the lateral influence extent (Li, Fig. 3) is 50 m for the shallow WT setting and 23 m for the deep WT configuration. In the case of the well, the influence radii (ri, Fig. 3) chosen were 25 and 17.5 m for shallow and deep WT settings, respectively.

Fig. 3
figure 3

Distribution of the a trenches and b wells in the hypothetical squared dam. (Ad: dam area; Li: trench lateral influence extent; ri: influence radius; W: trench width)

Full size image

The computed number of wells/trenches is multiplied by the yearly cumulative recharge of one well/trench. The resulting figures correspond to the annual cumulative recharge yielded by a battery of wells/trenches for both deep and shallow WT.

Results and discussion

Simulations of continuous infiltration

The evolution of the water content over time shows that the recharge from the well and the trench arrives considerably earlier than the recharge from the SF (Fig. 4). When the trench or the well are used, the recharge produced by either of these solutions takes place at around one day. On the other hand, infiltration from the SF starts after 81 days (8100% slower than with a well or a trench).

Fig. 4
figure 4

Evolution of water content over time for the three main shallow WT scenarios. The SF and the Trench + SF scenarios are in 2D domains, while the Well + SF scenario in a pseudo-3D domain. Note that the CL is fully saturated (Ɵs = 0.36)

Full size image

For the scenarios with a deep WT, the onset of recharge from either the trench or the well is observed around 20 days after the beginning of infiltration. When the water infiltrates solely through the SF, the arrival of the waterfront to the WT takes place around 450 days after the beginning of infiltration (2250% slower than with a well or a trench). The long arrival times obtained when relying solely on SF infiltration are corroborated by Sasidharan et al. (2020b) and Maples et al. (2019), who respectively point out that a CL and subsurface heterogeneity (which implies low conductivity facies) can delay recharge from 1 to more than 100 years.

In Fig. 4 (top row), note that the clogging layer has the effect of precluding full saturation in the VZ. When either the well or the trench is used, an additional portion of the vadose zone is seized to accommodate and infiltrate more water.

The times at which recharge starts point to the fact that the trench and the well effectively promote the bypassing of the CL and the VZ. This effect is reduced when the unsaturated zone increases but is still significant. The improvement of recharge would translate into more recharge of groundwater resources and fewer evaporation losses. However, it is important to bear in mind that the arrival times are not proportional to infiltration rates or cumulative recharge. As pointed out by several authors (Quichimbo et al. 2020; Brunner et al. 2009a, b) and shown later on in this article, deeper water tables, which imply longer recharge arrival times, are consistent with higher infiltration rates due to the increase of the hydraulic gradient in the vicinity of the water source.

Analysis of cumulative recharge

The comparison of the cumulative recharge from one trench, one well and a portion of wadi dam SF (Fig. 5) shows that the well is the best infiltration solution. It has the highest yearly cumulative recharge. Recharge from the well is around 1484% higher than the one from the SF, and between 336 and 463% higher than the yearly cumulative recharge from the trench. Nonetheless, the trench’s yearly cumulative recharge is still higher than the SF’s one, namely between 442 and 321%.

Fig. 5
figure 5

Visual comparison of the yearly cumulative recharge for the solutions evaluated in absolute values (left-hand-side axis) and as a percentage of the cumulative recharge computed in the SF scenario (right-hand-side axis). The scenarios are in a shallow WT setting. The numbers above the bars correspond to the yearly cumulative recharge in m3

Full size image

When the WT is deep, the well still exhibits the best performance, followed by the trench (Fig. 6). Nonetheless, the performance of the well is enhanced regarding the trench and the SF infiltration. In this case, the well yearly cumulative recharge is again one magnitude higher than the one from the SF (1642%), and six to eight times as high as the recharge from the trench (599–825%).

Fig. 6
figure 6

Visual comparison of the yearly cumulative recharge for the solutions evaluated in absolute values (left-hand-side axis) and as a percentage of the cumulative recharge computed in the SF scenario (right-hand-side axis). The scenarios are in a deep WT setting. The numbers above the bars correspond to the yearly cumulative recharge in m3

Full size image

The obtained results suggest that the best solution to increase groundwater recharge in wadi reservoirs would be the deployment of vadose zone wells, i.e., the combination of these wells with an infiltration basin. Other studies that analyse this combination support the present findings. Alrehaili and Tahir Hussein (2012) and Alataway and El Alfy (2019) find that the use of wells in wadi dams in the KSA reduces losses to evaporation. More specifically, Alataway and El Alfy (2019) conclude that wells can account for 44% more recharge and a decrease of 86% of the evaporation losses. In the same line, Sasidharan et al. (2020b) determine that the best solution to cope with a clogging layer in infiltration basins is the addition of drywells. Furthermore, the drywells offer the advantage of occupying little space, which makes them more suitable under area constraints when the number of solutions to be installed is relatively small.

Validation of the models with analytical equations

Overall, the recharge rates calculated with the simplified models are in the same order of magnitude as the infiltration rates from the analytical equations. In this way, the definitive models are, to a reasonable extent, validated.

The differences in recharge rates between the analytical equations and the simplified models of a well for a shallow WT show deviations of about 24% (Fig. 7). For the deep WT setting, these differences are even smaller, around 20% (Fig. 7). The discrepancies found could derive from the limitations of the analytical expressions. According to Stephens (1979), the infiltration rates obtained with Glover (1953) can differ up to 160% from the actual field values. Another reason that might account, to a minor extent, for the differences found is the fact that the values selected for H and r are at the limit of validity of the analytical equations. In this regard, Stephens (1979) pointed out that electric analogues conducted by Glover under the same conditions of the analytical equations yielded errors of 8% when H/r = 8, i.e., close to the limit of validity (H/r = 10).

Fig. 7
figure 7

Visual comparison of the recharge rates produced by analytical and numerical scenarios for the well in a a shallow and c deep WT setting, and the trench in a b shallow and d deep WT configuration. References include: Bouwer 2002; Elrick et al. 1989; Glover 1953; Heilweil et al. 2015; USBR 1977

Full size image

The trench in a shallow WT setting presents minor incongruences between analytical and numerical solutions (i.e. 5%; Fig. 7). Such small differences greatly validate the assumptions made for the conceptual models.

On the other hand, the highest difference during this analysis is found with the trench in a deep WT setting (Fig. 7). The infiltration rate obtained with the equation of Heilweil et al. (2015) is higher than the modelling results, with a difference of about 115% (Fig. 7). Such a difference might be a consequence of the simplifications implied in the equation by Heilweil et al. (2015), which does not consider, for instance, the depth of the trench, which is shown in this study (the sensitivity analysis) to affect the infiltration rates. Furthermore, the equation by Heilweil et al. (2015) could exaggerate the effect of a deeper water table given the fact that this formulation is empirical and based on a particular case study in which fractures play an important role and become open at greater depths (Heilweil et al. 2004). The equation by Bouwer (2002) is also limited since it does not account for the depth to the WT, which is proved in the sensitivity analysis of this study to be a vital parameter controlling trench infiltration rates.

Sensitivity analysis of well and trench dimensions

Sensitivity analysis

Figure 8 shows the sensitivity of the recharge rate concerning well depth and radius. It is evident that the well is more sensitive to changes in depth. Changes in the vertical dimension affect the QSSRRs more than changes in the horizontal direction (radius). For instance, doubling the length of the well in a deep WT setting increases the QSSRRs more than 230%, whereas doubling the radius, in contrast, only causes an increase of 135%.

Fig. 8
figure 8

Sensitivity analysis of the well. Recharge rates are compared with the recharge rate of the base scenarios (QSSRR quasi-steady-state recharge rate, L vertical length of the well; r radius of the well)

Full size image

The QSSRRs of the trench are likewise more sensitive to changes in its vertical dimensions (Fig. 9)—for example, when doubling the depth of the trench in a shallow WT the QSSRR increases to 151%. Contrastingly, changing the width of the trench has practically no effect on the QSSRRs, as attested by the horizontal lines in Fig. 9 for both shallow and deep WT settings. Heilweil et al. (2015) found similar results in their study. Based on trench sensitivity and the model calibration experiments, they concluded that the infiltration from the trench takes place predominantly through the trench walls and in a minor proportion through the floor.

Fig. 9
figure 9

Sensitivity analysis of the trench (QSSRR quasi-steady-state recharge rate, D depth of the trench, W width of the trench)

Full size image

Unlike the well, the QSSRR of the trench is dependent on the depth to the WT. Differences in the QSSRR between the deep WT and shallow WT scenarios range from 5 to 29%, and for most of the scenarios are above 17%. These differences coincide well with the observation from Heilweil et al. (2015), which conclude that the main factor controlling infiltration rates from the trench is the depth to the WT. Such a conclusion is stressed in their equation, which has the depth to the WT and the hydraulic conductivity as the only variables.

However, it has to be noted that the well was modelled in an axisymmetric domain, whereas the trench was modelled in a plane 2D-domain. This domain difference becomes evident when comparing the pressure heads in Fig. 4 (day 50). Here, the horizontal flow component in the well scenario is weakened with increasing the distance from the axis (due to the axisymmetric of the domain). On the other hand, the horizontal flow is not diminished in the trench scenario. This situation indicates a possible influence of the boundary condition set at the right margin of the domain. However, the vertical extent of the mounding at the left boundary is only 10 cm, a minimal value compared to the model dimensions. Furthermore, preliminary domain dimension tests showed that a longer domain yields slightly lower trench infiltration rates, which does not change the overall conclusions of this study.

Analysis of construction costs

For all the scenarios explored, the well always has the best yearly cumulative recharge per dollar (see Fig. 10). Furthermore, the optimum scenarios (i.e. the ones in which the well/trench has the highest QSSRRs per unit of volume) have the absolute best performance (~7 m3/$ for the trench and ~28 m3/$ for the well). Note that the well scenarios that involve a smaller radius (0.5r) yield the most recharge per dollar. In the case of the trench, the recharge is higher when vertical and horizontal dimensions are shortened.

Fig. 10
figure 10

Yearly cumulative recharge per dollar considering the scenarios of the sensitivity analysis (L vertical length of the well; r radius of the well; D depth of the trench, W width of the trench; optimum well: 2 L and 0.5r; optimum trench: 0.5 W and 0.5D)

Full size image

The cumulative recharge per dollar is the highest in the scenarios in which the infiltration solution has smaller dimensions (e.g. trench with half the original width) and the highest cumulative recharge per volume. A plausible reason for this pattern is that smaller dimensions translate into a steeper hydraulic gradient around the well/trench, increasing the specific recharge. Händel et al. (2014) arrived at an analogous conclusion. They observed that specific recharge rates (in their case per unit of screened area) diminished once the well approached the WT.

Well and trench performance under area constraints

In a shallow WT setting, 418 wells are fitted inside the 1-km2 square infiltration area. The corresponding yearly cumulative recharge is 2.82 × 106 m3. In the case of the trench, the combined length of trenches deployed within the hypothetical reservoir is 9,000 m, and the corresponding yearly cumulative recharge is 4.51 × 106 m3. The trench recharge is around 160% of the well recharge. The result indicates that a dam with a fixed area in a shallow WT setting can yield the highest yearly cumulative recharge when trenches are employed.

When the water table is deep, the resulting yearly cumulative recharge is 6.37 × 106 m3 and 6.23 × 106 m3 for the well and for the trench, respectively. In this case, the trench recharge is about 90% with respect to the well recharge. This result point to the fact that the trench can provide a higher recharge when the WT is close to the surface and that this ratio decreases as the depth to the WT increases. However, it must be kept in mind that the length of influence of the trench in a shallow WT setting was not fully addressed and that it might be longer, rendering the difference between trench and well recharge smaller.

Limitations

The numerical study conducted here simulates a limited part of the infiltration process taking place in the field. One of the most relevant operational obstacles for any infiltration solution, which is not examined in this study, is clogging. From the literature, it seems likely that this issue could be more severe in wells than in trenches or infiltration basins; furthermore, so far, there is no applicable method for regenerating drywells (Bouwer 2002). Trenches and infiltration basins, in contrast to drywells, can be redeveloped to a reasonable extent by mechanically removing the sediments accumulated on the bottom (Heilweil et al. 2015; Bouwer 2002).

Overall, the best way to deal with clogging is by conducting an adequate pretreatment of the infiltrating water. Edwards et al. (2016) and Sasidharan et al. (2018) elaborate on some pretreatment solutions for drywells. Pérez-Paricio (2000) summarises general recommendations to prevent different types of clogging in MAR systems. Following the findings from this author, physical clogging in superficial systems can be minimised to a reasonable level if the concentration of total suspended solids (TSS) in the infiltration water is lower than 10 mg/L. For wells, a TSS-content of 2 mg/L is recommended if the Ksat of the receptive medium is larger than 40 m/day. For smaller Ksat values (as in the present study), TSS should not be higher than 0.1 mg/L.

Another element which should be regarded in complementary research is a more detailed structure of cost given the simplicity of the assumptions considered. In reality, there are more factors involved in the cost of construction of a well or a trench than the extracted volume. For instance, the drilling costs can vary largely depending on the geological conditions and the planned well setup without considering the maintenance costs. Furthermore, the prices employed in the analysis are not updated. Still, the solution derived from this study presents a first estimation with respect to enhancements of surface infiltration.

In this study, the receptive medium is regarded as homogeneous, and the effects of changing the material properties are not assessed. In this sense, the available literature suggests that well, trench and wadi infiltration rates are positively correlated to higher hydraulic conductivity (Sasidharan et al. 2020a; Heilweil et al. 2015; Händel et al. 2014), while the sensitivity to other material parameters such as water content (θr), saturated water content (θs), n, and α is not significant (Händel et al. 2014; Heilweil et al. 2015; Quichimbo et al. 2020; Sasidharan et al. 2019).

Aquifer heterogeneity, on the other hand, has been proven to be an essential parameter when assessing recharge in the context of drywells and infiltration basins and is expected to vary largely in alluvial sediments (Brunner et al. 2009a, b). Sasidharan et al. (2020a) concluded that drywell recharge was higher in heterogeneous profiles than in homogenous ones, while Sasidharan et al. (2019, 2020b) proved that these results are explained by the interconnection of high-permeability lenses, which provide a larger surface area for infiltration. Maples et al. (2019) came to a similar conclusion studying the effect of heterogeneity in infiltration basins by means of numerical modelling. They found that interconnected coarse facies account for rapid infiltration, while fine-textured ones accommodate the majority of recharge in the long-term. Their numerical simulations also proved that the geological architecture of the subsurface plays a key role in groundwater recharge, showing variations of up to two orders of magnitude. Finally, Händel et al. (2014) found infiltration via small-diameter wells not to be affected as strongly by heterogeneity as infiltration basins, especially when the anisotropy of the sediment translates into lower vertical hydraulic conductivities (which is common in alluvial sediments (Song et al. 2009),

Even though the assumptions of this study render the models simple, they provide a fundamental understanding of the differences between the considered means of infiltration. Further research is needed to assess the role of the material properties on the infiltration solutions for wadi dams, and especially for trenches, which lack literature on the effect of subsurface heterogeneity. Finally, despite the fact that the present analysis is restricted to hydrogeological conditions typical of wadis in the KSA, this study constitutes a starting point for analogous research in other dryland areas around the world.

Conclusions

Implementing a vadose zone well or a trench can cause a much earlier visible recharge into the wadi aquifers than simple surface infiltration. Also, in terms of cumulative recharge, the wells and the trenches show significant advantages, with the wells being the best solution by far. In terms of geometry, the well and the trench are especially sensitive to changes in the vertical dimensions, and in the case of the trench, to the initial position of the water table. In general, the infiltration from the models presented here agrees well with analytical equations. Under reservoir area constraints, the trench provides a higher cumulative recharge per year, providing the WT is close to the surface. In a deep WT setting, the result is the opposite, and the wells can recharge slightly more than the trenches.

The trench and the well contribute to bypassing the VZ in terms of time until recharge starts, leading to higher volumes of recharge and eventually fewer evaporation losses. By using these structures, recharge starts between 2250 and 8100% faster than the recharge from the surface.

When the cumulative recharge of the well is compared with the cumulative recharge from the surface and a trench (with the same infiltration area), it has been shown that the well has the best performance. Its cumulative recharge over a year is between 336 and 825% of the cumulative recharge from the trench and between 1484 and 1642% of the simple SF-recharge. These values indicate that both the well and the trench help bypass the vadose zone in terms of recharge quantity. Furthermore, the well is especially effective in this regard.

In terms of cost, the well has the highest yearly cumulative recharge per dollar, recharging between 400 and 800% of the trench’s recharge, when equivalent modifications are compared.

The modelling approach explored in this study is validated using analytical equations. Infiltration rates obtained with the models are in the same order of magnitude as those obtained with the analytical equations. For the trench in a deep WT setting, a disagreement is found between both approaches; however, such discrepancy may very well be caused by the simplifications and assumptions involved in the definition of the analytical equations. Furthermore, the two equations considered for trench infiltration lack one of the two parameters for which a major sensitivity was found, namely depth to the WT and trench depth.

The well recharge rates have the highest sensitivity to changes in the vertical dimensions, while the trench recharge rates depend on both changes in the vertical dimension and the initial position of the water table. The highest recharge rates for the well (702 m3/day) and the trench (42 m2/day) are obtained when the vertical length is doubled, and the WT is deep.