In general, there are detailed theoretical mathematical formulations of rainfall recharge and groundwater flow based on hydro-mechanical principles that require good software uses (Walton, 1970; Bear, 1979; Freeze and Cherry, 1979). This section presents practical formulations for groundwater recharge and unconfined aquifer groundwater flow upstream of subsurface dams. The basic form of the water balance equation, input, I, output, O, and the change in the storage can be written in terms of, ΔS, as follows.

## 8.1 Groundwater recharge

Groundwater recharge calculation methodologies are not repeated herein, because there are different old and recent methodologies for this purpose in the open literature. For instance, Subyani and Şen (1991) proposed a groundwater recharge study based on the recharge-outcrop relationship in the Arabian Peninsula. One of the methods is to record groundwater level fluctuations in observation wells and keep rainfall records near these wells, preferably at upstream meteorological stations. It is logically plausible that there is a direct relationship between groundwater level fluctuations and rainfall amounts, provided there are no groundwater withdrawals, which must also be taken into account in the water balance equation (Walton, 1970). Şen (2020) proposed a probabilistic methodology for estimating groundwater level fluctuations from rainfall records. It is suggested that the cumulative probability distribution functions (CDFs) of groundwater fluctuation and rainfall amounts are interrelated, so that the further the rainfall moves away from the mean, the greater the groundwater recharge. Lerner et al. (1990), Lerner (1997) and Rushton (1997) discussed the channel water budget method by considering surface water gains and transmission losses on the bases of streamflow data.

Another study concerns outcrop rocks, particularly in arid regions, and groundwater recharge calculations from occasional and rather scanty rainfall events. It is stated by Subyani and Şen (1991) that there is good agreement between predictions based on groundwater velocity and those based on isotope information.

## 8.2 Groundwater movement

The combination of the water balance (continuity) equation with Darcy groundwater velocity law provides the groundwater flow that governs the differential equation. The groundwater velocity, v, can be expressed in terms of the hydraulic conductivity, k, in the porous medium and hydraulic gradient, i, as follows.

v = ki (3)

There are different methodologies to find the value of k from ready tables or pumping test results depending on the geological formation or practically from sieve analysis as,

k = Cd2 (4)

where C is a constant. This method is based on a representative grain-size distribution effective diameter, d, (Krumbein, 1932) and the initial slope and intersection of the grain-size distribution curve (Al-Yamani and Şen, 1992).

Provided that R is the rate of groundwater recharge and A is the ratio of volume of water withdrawn per unit horizontal area, then net production can be expressed in terms of water density, \({\rho }\), and areal extensions \(\varDelta \text{x}\) and \(\varDelta \text{y}\) as,

$${\rho }\left(\text{R}-\text{A}\right)\varDelta \text{x}\varDelta \text{y}$$

5

In unconfined groundwater reservoirs (aquifers), the production of water depends on the porous medium including specific yield and retention; where the specific yield is equal to the groundwater aquifer storage coefficient, S. The amount of released water, \(\varDelta {\text{V}}_{\text{w}}\), or added to storage can be calculated as \(\varDelta {\text{V}}_{\text{w}}\)= S\(\varDelta \text{A}\varDelta \text{h}(\text{x},\text{y})\), and therefore the temporal ratio of water can be expressed as,

$${\rho }\frac{\varDelta {\text{V}}_{\text{w}}}{\varDelta \text{t}}= {\rho }\text{S}\varDelta \text{x}\varDelta \text{y}\frac{\varDelta \text{h}(\text{x},\text{y})}{\varDelta \text{t}}$$

6

where h is the groundwater level and t is the time. The substitution of this expression in continuity (water balance) equation leads to,

$$-{\rho }\text{d}\left(\frac{\partial {\text{V}}_{\text{x}}}{\partial \text{x}}+ \frac{\partial {\text{V}}_{\text{y}}}{\partial \text{y}}\right)\varDelta \text{x}\varDelta \text{y}+ {\rho }\text{U}\varDelta \text{x}\varDelta \text{y}-{\rho }{\text{V}}_{\text{z}2}\varDelta \text{x}\varDelta \text{y}+ {\rho }\left(\text{R}-\text{A}\right)\varDelta \text{x}\varDelta \text{y}= {\rho }\text{S}\varDelta \text{x}\varDelta \text{y}\frac{\varDelta \text{h}(\text{x},\text{y})}{\varDelta \text{t}}$$

7

where d is the aquifer thickness, U and D, are the upward and downward leakage (vertical) groundwater movement rates. It is possible to simplify the equation for two-dimensional groundwater flow motion as follows, assuming the water density as 1 gr/cm3 and as a single point.

\(-\text{b}\left(\frac{\partial {\text{V}}_{\text{x}}}{\partial \text{x}}+ \frac{\partial {\text{V}}_{\text{y}}}{\partial \text{y}}\right)+ {\text{V}}_{\text{z}1}-{\text{V}}_{\text{z}2}+ \left(\text{R}-\text{A}\right)=\text{S}\frac{\partial \text{h}}{\partial \text{t}}\) (8)

Darcy x- and y-direction groundwater flows and upward and downward groundwater table fluctuations can be expressed by the following equations.

$${\text{V}}_{\text{x}}= -{\text{K}}_{\text{x}\text{x}}\frac{\partial \text{h}(\text{x},\text{y})}{\partial \text{x}}$$

9

$${\text{V}}_{\text{y}}= -{\text{K}}_{\text{y}\text{y}}\frac{\partial \text{h}(\text{x},\text{y})}{\partial \text{y}}$$

10

$$\text{U}={\text{K}}_{1}\frac{\left[{h}_{1\left(x,y\right)}-h\left(x,y\right)\right]}{{\text{b}}_{1}}$$

11

and

$$\text{D}={\text{K}}_{2}\frac{\left[h\left(x,y\right)- {h}_{2}\left(x,y\right)\right]}{{\text{b}}_{2}}$$

12

Substitution of Darcy’s equation with these leakage terms in Eq. (8) yield,

\(\frac{\partial }{\partial \text{x}}\left[\text{b}{\text{K}}_{\text{x}\text{x}}\frac{\partial \text{h}\left(x,y\right)}{\partial \text{x}}\right]+ \frac{\partial }{\partial \text{y}}\left[\text{b}{\text{K}}_{\text{y}\text{y}}\frac{\partial \text{h}\left(x,y\right)}{\partial \text{y}}\right]+ {\text{K}}_{1}\frac{\left[{\text{h}}_{1\left(x,y\right)}-\text{h}\left(x,y\right)\right]}{{\text{b}}_{1}}+ {\text{K}}_{2}\frac{\left[\text{h}{\left(x,y\right)- \text{h}}_{2}\left(x,y\right)\right]}{{\text{b}}_{2}}+ \text{R}-\text{A}=\text{S}\frac{\partial \text{h}\left(x,y\right))}{\partial \text{t}}\) (13)

Finally the main equation for two-dimensional flow in an unconfined aquifer takes the following form.

$$\frac{\partial }{\partial \text{x}}\left[\left(\text{h}- {\eta }\right){\text{K}}_{\text{x}\text{x}}\frac{\partial \text{h}\left(x,y\right)}{\partial \text{x}}\right]+ \frac{\partial }{\partial \text{y}}\left[\left(\text{h}- {\eta }\right){\text{K}}_{\text{y}\text{y}}\frac{\partial \text{h}\left(x,y\right)}{\partial \text{y}}\right]+ {\text{K}}_{1}\frac{\left[{h}_{1}\left(x,y\right)-h\left(x,y\right)\right]}{{\text{b}}_{1}}+ {\text{K}}_{1}\frac{\left[h\left(x,y\right)-{h}_{2}\left(x,y\right)\right]}{{\text{b}}_{1}}+ \text{R}-\text{A}={\text{S}}_{\text{y}}\frac{\partial \text{h}(\text{x},\text{y})}{\partial \text{t}}$$

14

where \(\eta\)is the bottom elevation. The unknown h(x, y, t) in the above equations can be determined using suitable method and boundary and initial conditions.

Although each month generates a unique monthly average, the set of monthly averages remains the same for all designs. Therefore, the supply-side quantities for each case study and the generic model differs as a function of the size of the catchment area. The demand side of the equation has been developed specifically for each project, based on water demands for local use and landscape needs. The application of these equations was presented by Önder and Yılmaz (2005) for unconfined groundwater flow in subsurface dams