Next Article in Journal
The Role of Carbon to Nitrogen Ratio on the Performance of Denitrifying Biocathodes for Decentralized Wastewater Treatment
Next Article in Special Issue
Assessing the Impacts of Climate Change on Rainwater Harvesting: A Case Study for Eight Australian Capital Cities
Previous Article in Journal
A GIS-Based Model for Flood Shelter Locations and Pedestrian Evacuation Scenarios in a Rural Mountain Catchment in Romania
Previous Article in Special Issue
A Review and Analysis of Water Research, Development, and Management in Bangladesh
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Integrated Surface Water and Groundwater Modeling in Arid Environment, Al-Lusub Watershed, Saudi Arabia

Department of Hydrology and Water Resources Management, Faculty of Meteorology, Environment & Arid Land Agriculture, King Abdulaziz University, Jeddah 21589, Saudi Arabia
*
Authors to whom correspondence should be addressed.
Water 2022, 14(19), 3075; https://doi.org/10.3390/w14193075
Submission received: 25 August 2022 / Revised: 24 September 2022 / Accepted: 24 September 2022 / Published: 29 September 2022
(This article belongs to the Special Issue Sustainable Water Futures: Climate, Community and Circular Economy)

Abstract

:
Coupling surface and groundwater for a better understanding of the hydrologic response in arid basins is essential for managing water resources. The study aims to integrate surface water (SW) and groundwater (GW) modeling in an arid environment and to identify groundwater sources at Hadat Ash-Sham farm station in the Al-Lusub watershed, located in the western part of Saudi Arabia. The curve number (CN) method for SW modeling is primarily used to estimate GW potential recharge based on a developed correlation equation between event-based rainfall and potential recharge, which is processed in the Hydrologic Engineering Center-Hydrologic Modeling System (HEC-HMS) software. Monthly potential recharge is utilized to compute GW modeling in MODFLOW. Van Mullem’s equation is used to calculate hydraulic conductivity (K) from CN. The calibration of steady-state GW modeling reveals that the K (from Van Mullem’s equation) is within the range of measured K in the field. The transient groundwater modeling phase concludes that the groundwater system in the Al-Lusub watershed can be interpreted based on two different scenarios. The first is that an extra recharge from a nearby watershed flows through the geological structure into a downstream area. The second scenario involves potential recharge from surface water flowing on the Al-Lusub watershed’s mainstream bed. Validation reveals a strong relationship between predicted and observed water tables. The root mean square error (RMSE) for Scenarios 1 and 2 are 0.6 m and 0.7 m, respectively. Further investigation is needed to determine the region’s most common scenario.

1. Introduction

Nowadays, the world faces many challenges concerning the continuously increasing human population. One of the essential freshwater sources is groundwater. Although there are water desalination facilities, people cannot always rely on them, since it is too expensive. The dense human population in such an area impacts the land-use change due to human activities to fulfill their needs, e.g., agriculture, manufacturing industry, housing, etc. Under these circumstances, there is an increase in freshwater demand to supply for consuming sectors. Water resources management research may involve the application of many disciplines, such as hydrogeology, groundwater hydrology, and surface water hydrology [1].
Saudi Arabia is like most of the arid region countries with low rainfall and no permanent water river. This low amount of rainfall is the primary source of recharge for the groundwater system. This circumstance of an arid environment with high freshwater demand makes groundwater the most precious thing on earth. Fitts [2] stated that groundwater had become a vital source for the future.
Based on this information, it is crucial to study the correlation between surface water and groundwater on a large scale (i.e., watershed/catchment scale), especially in an arid environment. For the large scale (regional or watershed-scale), the study needs to integrate surface water and groundwater as a holistic study.
The integrated surface and groundwater studies in a watershed can be achieved through the modeling process. The simulation should describe, to some extent, the actual condition; thus, it is essential to include the validation process at the end of the modeling process to check model reliability. The model is then used for water resources management purposes.
Water resources research involving multiple disciplines, such as hydrogeology (groundwater hydrology) and surface water hydrology, is difficult to conduct today [1]. Both subjects are consistently identified as distinct disciplines. The relationship between groundwater hydrology and surface water hydrology is established by calculating the recharge from each rainfall event. The groundwater recharge is the amount of precipitation remaining after subtracting evaporation and surface runoff. The water then leaves the soil until it reaches the saturated zone [3,4].
For the large scale (watershed scale), integrated surface water and groundwater modeling typically employ loosely coupled models or iterative solution methods. Each model is interconnected by feeding the output of one model into another [5]. The modeling methods have been widely used for half a century, but the majority of them are applied in humid regions rather than the arid environment [6]; therefore, it is crucial to study the modeling application in the arid environment.
The objective of this study is developing an integration methodology between surface water and groundwater in the Al-Lusub watershed to interpret the existence of the water resources in the basin. The CN number method used in surface water modeling also determines the hydraulic conductivity value of the wadi alluvium and fractured rock. The coupling between surface and groundwater modeling is achieved through the development of a monthly potential recharge from rainfall, which is used as input for groundwater modeling. In the calibration phase, it is essential to define the realistically expected scenarios that provide a minimum model error and also provide a realistic interpretation of the measured groundwater response in the Al-Lusub watershed.

2. Study Area

Al-Lusub watershed lies in the western part of Saudi Arabia (Figure 1). It is the upstream part of the Usfan watershed, where the runoff water will continue to flow to the Red Sea. The geographical coordinate of this watershed is 39°39′18″ E–40°15′25″ E and 21°43′19″ N–22°10′7″ N. In the downstream area of the watershed, there is a farm that belongs to King Abdulaziz University, called the Hadat Ash-Sham Farm Experimental Station. This farm manages plantations (e.g., dates) and animals (e.g., cows, chickens, honey, etc.). For this purpose, the farm utilizes groundwater as the primary water source.
Rock outcrops cover the surface area of the watershed and are mostly impervious. The vegetation covers a small percentage of the alluvium area. The residential, agriculture and stream bed cover 1%, 14%, and 20% of the alluvium, respectively. Al-Lusub watershed soils have the potential for agriculture purposes [7].
Table 1 shows the morphometric information of the study area, which has an area of 901.16 km2, a perimeter of 270 km, a basin length of 72.02 km, a total stream length is 710.5 km, and a maximum stream length is 96.80 km. The sinuosity number can be derived by dividing the maximum stream length by basin length, and the result is 1.34. Based on Charlton [8], this value is classified as high sinuosity, which indicates the stream in the Al-Lusub watershed is a meandering type. Al-Lusub watershed lies 62–1765 m above sea level (Figure 2). The drain density is <2, which is classified as very coarse drainage [9]. The hypsometric curve number is 0.48, and it indicates the mainstream in the catchment is in the mature stage [10].

3. Materials

There are nine rainfall stations located surrounding Al-Lusub watershed. Only five stations have good quality and reliable data on these nine rainfall stations. The location of these five stations in Al-Lusub watershed can be seen in Figure 3. Using the Thiessen Polygon method, it can be concluded that only 2 rainfall stations influence the study area (i.e., J-104 and J-214). The Thiessen influence percentage for J-104 and J-214 are 25% and 75%, respectively. Limited data circumstance is typical for arid environments; therefore, it is important to use existing data in appropriate ways to solve water issues in the arid region.
The maximum rainfall data for each month and monthly rainfall, in Hadat Ash-Sham and Madrakah (i.e., J-104 and J-214, respectively) are shown in Figure 4. These data focus on a time window from November 2018–April 2021, since it will be used for the analysis and comply with the availability of other data (i.e., well data).
El-Hames [11] collected information from 60 wells across Al-Lusub watershed. The information includes water table elevation and the bottom of the shallow. The map location of these wells can be seen in Figure 5. Sharaf et al. [12] also determined that the average depth of the shallow aquifer in the Al-Lusub watershed downstream is 50 m. This conclusion is based on geophysical research conducted in Hadat Ash-Sham.
Hadat Ash-Sham is located in the Al-Lusub watershed’s downstream region. The average height of the alluvium in the region is 50 m [12]. Al-Amri [13] measured hydraulic conductivity utilizing two techniques: grain size analysis and infiltration test at the south of the farm, along the seasonal streams.
Alyamani and Şen [14] present a comprehensive grain-size analysis for calculating hydraulic conductivity. Al-Amri [13] collected 50 samples along the seasonal streams in the Hadat Ash-Sham region. These samples were also collected from different depths. The grain size analysis used the Hazen formula, which is an empirical method, and the equation can be seen in Equation (1).
K = C ( d 10 ) 2      
where K is the saturated hydraulic conductivity (cm/sec), d10 is the effective size (cm), and C is a constant considered equal to 100.
Al-Amri [13] also performed the constant head method, using the constant head permeability device in the lab to measure the hydraulic conductivity based on collected samples from the exact location. Then, the saturated K is calculated using Darcy’s (Equation (2))
K = Q L A H
where Q is the volumetric flow rate (cm3/min), L is the sample height (cm), A is π multiplied by squared sample radius, and H is the difference of heads (cm). This constant head method is suitable for measuring the saturated K in permeable soils (e.g., sand and gravel).
Al-Amri [13] experiments in Hadat Ash-Sham can be seen in Table 2. The result shows that for grain-size analysis, the K range is 0.95–35.17 m/day, while for the constant head method, the K range is 4.25–55.32 m/day. The average K for grain-size analysis and constant head method is 10.62 and 19.48 m/day. The median of K for grain-size analysis is 7.72, and the median of K for the constant head method is 18.49.
Many researchers also performed previous studies to investigate the K of alluvium. Al-Khatib [15] estimated the K alluvium to be equal to 103.68 m/day, Hussein and Bazuhair [16] stated that the K alluvium has ranged between 5.6–62.6 m/day, and Sharaf et al. [12] investigated the K alluvium equal to 95.04 m/day. El-Hames [11] performed the constant-head method, which is similar to Al-Amri [13], and the result was 86.4–216 m/day. The compilation of this information can be seen in Table 3. This information shows that hydraulic conductivity is generally highly variable [17], especially in arid environments. In the end, calibration is used to see the most representative K value of the alluvium.
In the downstream of Al-Lusub watershed, King Abdulaziz University installed an experimental farm for research purposes called Hadat Ash-Sham Farm Station. Inside the farm, 11 wells are actively being monitored (Figure 6). In the end, these data will be used for calibration and validation purposes during the transient groundwater modeling process. Figure 7 shows the graphs describing the water table changing in Hadat Ash-Sham Farm.
Budiman et al. [19] studied geostatistical analysis for this water table data from November 2018 to February 2020. The result shows that the groundwater flow directions have a typical pattern during the dry and wet seasons. During the dry season, the groundwater flows toward the outside of the farm in random directions. On the other hand, the groundwater flows toward the northwest direction during the wet season. This flow direction corresponds to the primary stream located south of the farm and recharges the farm area.

4. Methods

The conceptual model of this study can be seen in Figure 8. The watershed has two geological features on the surface: alluvium and fractured rock. The alluvium lies surrounding the stream channel, while the fractured rock is a rock that is close to the surface. The alluvium is the youngest sedimentary rock formed by the materials carried by the streams. The rock fracture is possibly caused by intensive weathering on the surface. Underneath both of these features, the rock is still in intact condition. Based on this information, the higher hydraulic conductivity occurs in alluvium and fractured rock, while the impermeable occurs in the intact rock. In this study, both the alluvium and fractured rock represent shallow aquifer.
If rainfall occurs on a dry stream channel in a watershed, the flow volume is reduced by infiltration in the bed, the banks, and the flood plain. These infiltration losses are called transmission losses [20]. The runoff will flow on the streams until reaching the outlet. The transmission loss will become potential recharge for the groundwater system until reaching the water table on the shallow aquifer.
The water table in the watershed will be high during the rainy season, and the rainfall intensity duration influences its raising, while during the dry season, the water table will drop to the lowest level.
Figure 9 shows the methodology flow chart implemented in this study. It is divided into three main stages: rainfall analysis, surface water modeling, and groundwater modeling (i.e., steady-state groundwater modeling and transient groundwater modeling). Rainfall analysis aims to derive the appropriate storm hyetograph patterns in the study area. Based on these hyetographs, the rainfall will be inputted for surface water modeling. The surface water modeling utilizes SCS curve number (CN) method in HEC-HMS software. The CN values used for surface water modeling will also be used in groundwater modeling to define hydraulic conductivity (K) using Van Mullem’s equation [21].
Surface water modeling aims to see the relation between rainfall and potential recharge. This relation is essential to estimate the potential recharge based on monthly rainfall data inputted for groundwater modeling. This information shows that the goodness of groundwater modeling depends on surface water modeling. The groundwater modeling comprises steady-state groundwater modeling and transient groundwater modeling. The purpose of steady-state groundwater modeling is to review the appropriate hydraulic conductivity (K). The purpose of transient groundwater modeling is to understand the watershed’s hydrology based on possible scenarios in the watershed.
This coupling of surface water modeling and groundwater modeling using HEC-HMS and MODFLOW 6 can be considered new in the arid environment, especially in Saudi Arabia. Moreover, it utilizes the same CN number as the parameter input for both models.

4.1. Rainfall Analysis

The most common storm duration in Saudi Arabia is 3 h [22]. Elfeki et al. [23] proposed the appropriate dimensionless cumulative hyetograph for some areas in Saudi Arabia. The adequate dimensionless cumulative hyetograph for Al-Lusub watershed can be seen in Figure 10, dedicated for Jeddah City on the first quartile.
Based on this dimensionless cumulative hyetograph, maximum daily rainfall data recorded by rainfall stations J-104 and J-214 (Table 2) are processed to develop hyetographs. Rainfall data with zero (0) rainfall are not required to be processed, and these increment hyetograph data will be inputted as rainfall during the surface water modeling stage.

4.2. Surface Water Modeling

Before performing surface water modeling, the first thing to do is delineate the watershed boundary using WMS. This watershed has a morphometric characteristic that WMS automatically exports this data for surface water modeling in HEC-HMS.
Then, it requires delineation of the alluvium inside the watershed. This step is done in the Google Earth application, where the actual surface condition can be seen; thus, alluvium delineation can be appropriately estimated. The result of this delineation is the separation between alluvium and rock boundaries (Figure 11). The total area of alluvium is 130.80 km2 (14.52%) and rock is 770.35 km2 (85.48%).
In performing integrated surface and groundwater modeling, it is essential to know the relationship between runoff and potential recharge when rainfall occurs in Al-Lusub watershed. This analysis will be computed using HEC-HMS software. In HEC-HMS, the loss method that will be used is SCS curve number (CN) method, and the transform method is SCS Unit Hydrograph.
Hawkins [24] introduced a well-known formula to estimate runoff using CN values (Equations (3)–(5)), where λ = 0.2.
Q = ( P I a ) 2 ( P I a ) + S
S = 25400 C N 254
I a = λ S = 0.2 S
The value of λ = 0.2 is insufficient to represent the initial abstraction ratio in an arid environment [25,26]. Farran et al. [27] studied that λ = 0.01 is appropriate to represent the basin in Saudi Arabia.
Farran et al. [28] explained that for arid environments, such as Saudi Arabia, there is a correlation between infiltration (F) and rainfall (P), as shown in Equation (6).
F = 77 %   P
This equation will be utilized as a reference during surface water modeling in Al-Lusub watershed. The ratio of infiltration and rainfall must be close to 77%.
Farran et al. [27] studied the most suitable CN number for alluvium and rock in arid regions, such as Saudi Arabia. They recommended CN alluvium is 40 and CN fractured rock is 72. These numbers will be used in surface water modeling to estimate recharge water, and Equation (7) will be utilized to calculate the CN average used for surface water modeling.
C N a v e r a g e = C N a l l u v i u m   A a l l u v i u m + C N f r a c t u r e d   r o c k   A f r a c t u r e d   r o c k A a l l u v i u m + A f r a c t u r e d   r o c k                                                      
In the arid region, Elfeki et al. [22] found that there is a correlation between CN and K. This correlation is presented using the empirical equation of Van Mullem’s [21] (Equations (8) and (9)).
C N = 6.2263 ln ( K 13 , 083 )    
K = 21 , 462 exp ( 0.16061   C N ) 0.6096
where K unit is m/day. By substituting CN alluvium and CN fractured rock (i.e., 40 and 72, respectively) into Equation 9 the hydraulic conductivity for alluvium and fractured rock in Al-Lusub watershed can be calculated. The K values derived from these steps will be input into the groundwater modeling step.
In HEC-HMS, it is required to input the lag time, which correlates with the time of concentration (Tc). Equation (10) shows how to calculate lag time from Tc [20]. There are many equations on how to calculate the time of concentration (Tc). Albishi et al. [29] introduced the formula to calculate Tc in an arid region with specifications in Saudi Arabia. The Tc formula can be seen in Equation (11). This equation is used to calculate appropriate Tc in Al-Lusub watershed.
Lag   Time = 0.6   Tc
Tc = Basin   Length 0.09 Basin   Slope 0.11
Surface water modeling in HEC-HMS can be run using maximum daily rainfall data for each month, average CN value, and lag time. The result of this stage is the correlation equation of the rainfall and potential recharge. This correlation equation will be used to calculate potential monthly recharge for transient groundwater modeling based on monthly rainfall.

4.3. Groundwater Modeling

In this study, groundwater modeling is focused on the shallow unconfined aquifer. The shallow aquifer comprised alluvium and fractured rock close to the surface, and it is believed that the rock close to the surface has fractures that give higher K than deeper rock.
Groundwater modeling is divided into two stages. The first stage is steady-state groundwater modeling, and the second stage is transient groundwater modeling.
Before running the MODFLOW 6 (with ModelMuse interface) model, producing the preparation or adjustment is required. The adjustment is divided into three stages (excluding boundary condition setup): unit and layer setup, grid construction, solver setup, and cell type setup for the unconfined aquifer.

4.3.1. Unit and Layer Setup

In this study, the time unit is in days, and the length unit is in meters (m). The total layer applied for this study is only one layer because the study will focus on the shallow unconfined aquifer.

4.3.2. Grid Construction

Grid construction during steady-state groundwater modeling focuses on the watershed area; therefore, the shapefile of the watershed boundary is imported into MODFLOW 6. The grid cell size is 200 m × 200 m (Figure 12).
The next step is to input the aquifer depth distribution, obtained from El-Hames [11]. This depth (from surface) distribution is divided into four groups, where the group delineation is according to the actual satellite image and considering the lithology of the geological map (Figure 13). Each group will provide average aquifer thickness from the surface. This method is chosen because if it uses the interpolation technique (e.g., IDW and Kriging), the contour sometimes crosses until above the surface, especially on the outer boundary.
The average depth for groups 2, 3, and 4 is 8 m, 13 m, and 19 m, respectively. The depth information for group 1 uses the information from another study. Sharaf et al. [12] investigated the depth of this shallow aquifer at the downstream of Al-Lusub watershed, using the geophysical technique, and the result is a 50 m depth average.
The aquifer depth information is used to calculate bottom aquifer elevation. In ArcGIS, the Raster Calculator tool is used to calculate aquifer elevation using Equation (12).
S h a l l o w   A q u i f e r B o t t o m   E l e v a t i o n = G r o u n d   S u r f a c e   E l e v a t i o n A q u i f e r   T h i c k n e s s  
The surface information is derived from Digital Elevation Model (DEM). After the elevation of the shallow aquifer bottom is generated, it is converted into a shapefile (i.e., points). This shapefile that contains aquifer bottom elevation is imported into MODFLOW 6 and adjusted to generate the bottom elevation layer.

4.3.3. Solver Setup

In this study, the MODFLOW 6 default solver setup is adjusted. Complexity is adjusted into complex instead of simple or moderate, and “Use Newton Formulation” option is checked. Linear acceleration is adjusted into BICGSTAB. The option of “Continue even if no convergence” can be checked if there is difficulty in reaching convergence during transient groundwater modeling.

4.4. Steady-State Groundwater Modeling

The purpose of steady-state groundwater modeling is to calibrate some physical properties (e.g., hydraulic conductivity, conductance, etc.) for each geological feature (i.e., alluvium and fractured rock). Many researchers have done several studies to investigate the value of K alluvium in the Al-Lusub watershed (Table 3 and Table 4).
The equation for the 2D groundwater model in steady-state conditions can be seen in Equation (13) [30,31].
x [ T x x h x ] + y [ T y y h y ] + R E = 0
where Tii is transmissivity for i = x, y directions, h is the hydraulic head, R is the recharge, and E is the evaporation.
Because the study will be applied to an unconfined aquifer, then T is given by Equation (14).
T = K ( h Z )
where Z is the bottom aquifer elevation, it is known that CN values for alluvium and rock are 40 and 72, respectively [27]. By using Equation 4.6 [21], the K value can be derived. K for alluvium is 21.2 m/day, and K fractured rock is 0.12 m/day (Table 4). K alluvium is still within the range of Al-Amri [13] and Hussein and Bazuhair [16]. Moreover, it is almost similar to K alluvium’s median, based on constant head analysis by Al-Amri [13] (Table 4).
In this study, K fractured rock is 0.12 m/day or 1.43 × 10−6 m/second. Based on literature book information from Freeze and Cherry [32], they explained that the K of permeable basalt range is 10−1–10−7 m/second, and the K of fractured crystalline rocks is 10−4–10−9 m/day; therefore, the K of fractured rock in this study is reliable.
In the 3D groundwater model, the direction of the K must be considered, and it is called the anisotropy of the K in different directions. The direction can be distinguished by being either radial (Kr) and vertical (Kz). Todd [33] studied that for the alluvium layer, Kz/Kr has a range between 0.1 to 0.5.

4.4.1. Steady-State Boundary Conditions

Boundary conditions are essential parameters to be inputted into MODFLOW 6. In this study, several boundary conditions need to be considered during steady-state groundwater modeling. These boundary conditions are as follows:
  • CHD: Time-Variant Specified-Head Package
This study applies the CHD package for the well on the upstream area and watershed outlet (Figure 14).
2.
DRN: Drain
In this study, the DRN package is applied to the stream network in the watershed (Figure 14). In the DRN package, one parameter can be calibrated, i.e., conductance.

4.4.2. Steady-State Validation

Validation is a stage that requires a comparison between the actual measurement and predicted result using statistical analysis. These well locations will record simulated water table changes during the modeling process. Another package in MODFLOW 6 will be utilized to perform this simulation called OBS. OBS is an observation utility package pane that allows users to observe the water table or hydraulic head value in the desired cell (i.e., observation well). Figure 15 shows the MODFLOW 6 visualization of the observation wells in this study. Most of the wells are located in the alluvium area.
One of the standard methods utilized to see the goodness of the model performance is using root mean squared error (RMSE) [34]. A good model is revealed by the RMSE value that must be close to zero [35]. The equation of RMSE can be seen in Equation (15).
R M S E = 1 n i = 1 n ( h ^ i h i ) 2
where h ^ i is the predicted data, h i is the actual observed data, and n is the sampling point numbers.
Another method is using R2 for the linear regression line, combining the data of observation and prediction. The R2 that is close to 1 indicated that there is a good relationship between predicted and actual observation.

4.5. Transient Groundwater Modeling

Transient groundwater modeling aims to simulate the groundwater system within the time sets, based on specific physical properties, after inputting several boundary conditions. Some result parameters will be monitored, such as water table or hydraulic head changing. The water table change must be similar to the actual change that occurred on the field. Some physical properties and boundary conditions can be calibrated based on the simulated water table changing to match the actual conditions.
The equation for transient groundwater modeling can be seen in Equation (16) by McDonald and Harbaugh [36]:
x [ T x x h x ] + y [ T y y h y ] + R E W = S h t
The hydraulic conductivity of the alluvium and the fractured rock is the same as the values applied for steady-state groundwater modeling (Table 5). The values are derived from specific CN values. There is another physical property that must be inputted, i.e., storativity. In an unconfined aquifer, the storativity includes the specific yield (Sy) or drainable porosity, given by following Equation (17).
S = S y + S s b
where S is storativity (dimensionless), Sy is specific yield (dimensionless), Ss is specific storage (L−1), and b is aquifer thickness (m). In an unconfined aquifer, Sy >> Ssb; therefore, it can be assumed that S ≈ Sy.
In Hadat Ash-Sham (the downstream area of the Al-Lusub watershed), Hussein and Bazuhair [16] studied that storativity for alluvium in this area ranges from 0.001 to 0.128. Lohman [37] examined that for unconfined aquifers, the storativity range is typically between 0.1 to 0.3.

4.5.1. Transient Time Setup

The time setup during transient groundwater modeling is made on a daily basis because there is no monthly basis option in MODFLOW 6. This daily basis must follow this study’s actual total days between 1 November 2018 and 30 April 2021. At the end of each month, the recharge will be stored, and there will be one day at the end of each month where recharge will occur. In this study, MODFLOW 6 software will run for 912 days, representing the daily time step from 1 November 2018 until 30 April 2021.

4.5.2. Transient Boundary Conditions

Several boundary conditions will be applied for transient groundwater modeling in this study. Each boundary package is used based on a selected scenario, and not all boundary packages will be used for each scenario.
  • Time-Variant Specified-Head Package
In the transient groundwater modeling of this study, the CHD boundary package will be applied for two locations, i.e., the upstream alluvium area and the outlet of the watershed. Due to there being no water table data in the upstream area from November 2018 until April 2021, it is assumed that the CHD at the upstream alluvium is equal to surface alluvium minus 1 m (1 m below the surface). CHD at the watershed outlet during transient groundwater modeling is identical to the CHD at the outlet during steady-state groundwater modeling. In the end, all CHD will be required to go through the calibration step until reaching a minor error.
2.
DRN: Drain
In transient groundwater modeling, the DRN boundary package will be applied to the stream network. The conductance value will be determined using calibration results from the steady-state groundwater modeling.
3.
ETS: Evapotranspiration Segments Package Pane
In this package, the user can adjust the evapotranspiration rate (m/day) and the influenced zone. Al-Amri [13] studied that the evaporation rate in Hadat Ash-Sham is almost equal to 10 mm/day or 0.01 m/day. ETS package is applied for the whole Al-Lusub watershed and all-time steps.
4.
RCH: Recharge Package Pane
In this study, the recharge is added using the potential monthly recharge based on monthly rainfall. The correlation equation of potential recharge and maximum daily rainfall on each month will be used to calculate potential monthly recharge. This correlation equation is the result of surface water modeling.
This study applies the potential recharge from monthly rainfall events for the whole watershed. The recharge might also be from other sources, influenced by certain geological conditions, such as geological structure and high permeability layer from neighboring watersheds. This recharge also will be applied for a specified time step, i.e., at the end of each month
The appearance in MODFLOW 6 that uses selected boundary conditions for transient groundwater modeling can be seen in Figure 16.
5.
Transient Validation
Validation is a stage to compare actual measurement and predicted results using statistical analysis. In this study, the transient groundwater modeling uses Hadat Ash-Sham Farm Station observation, where the actual data are recorded every week from November 2018 to April 2021. The average monthly water table data in Hadat Ash-Sham is shown in Figure 7B.
During the modeling process in MODFLOW 6, one observation well is placed at Hadat Ash-Sham Farm Station. This well will be recording water table changes in Hadat Ash-Sham. A package in MODFLOW 6 will be utilized to perform this simulation, called OBS. OBS is an observation utility package to observe the water table or hydraulic head value in the desired cell (i.e., observation well). Figure 17 shows the MODFLOW 6 visualization of the observation wells during transient groundwater modeling, located at Hadat Ash-Sham Farm Station.
The goodness of the model performance will be justified by root mean squared error (RMSE) in transient groundwater modeling. A good model is revealed by the RMSE value that must be close to zero [35].

5. Results and Discussions

5.1. Surface Water Modeling Results

It is known that the total area of Al-Lusub watershed is 901.15 km2, where the area of alluvium is 130.80 km2 (14.52%), and the area of rock is 770.35 km2 (85.48%). CN average in the Al-Lusub watershed is calculated by using the study of Farran et al. [27], where CN alluvium is equal to 40, and CN fractured rock is equal to 72. Initial abstraction (Ia) is calculated using λ = 0.01, where the result is 1.23 mm (Table 5). The result shows that the CN average of the Al-Lusub watershed is 67.36.
Tc and Tlag are calculated using Equations (10) and (11), and the Tc value is equal to 1.88 h, and the Tlag is equal to 67.69 min (Table 6). All this information is inputted into HEC-HMS to run surface water modeling.

5.1.1. Hyetograph and Hydrograph

The example results of surface water modeling are shown in Figure 18. The result shows the typical hydrograph in the arid region, characterized by a steep rising limb (a short time to reach the peak) and rapid recession to zero flow [22]. Figure 18 shows no runoff conditions occurred in November 2019. These conditions occurred because the rainfall was lower than the initial abstraction (1.23 mm).

5.1.2. Potential Recharge

Figure 19 shows the correlation between rainfall and recharge in Al-Lusub watershed for CN equal to 67.36. Equation (18) represents the correlation between rainfall (P) and potential recharge (R). Statistically, the R2 value (0.9935) is close to 1, where it is concluded that rainfall and recharge are highly correlated.
R = 0.7986 P
Equation (18) is almost similar to the study done by Farran et al. [28] where F = 77% P.
The next step is to calculate potential recharge using monthly rainfall data. The total rainfall is derived based on the Thiessen polygon, where the influence of J-104 is equal to 25%, and the influence of J-214 is equal to 75%. Potential monthly recharge is derived based on Equation (18), and the result can be seen in Table 7. This potential monthly recharge will be utilized in transient groundwater modeling and applied for the whole watershed.

5.2. Groundwater Modeling Results

5.2.1. Steady-State Groundwater Modeling

The result of steady-state groundwater modeling can be seen in Figure 20. The graph comprised actual water table measurement [11] and predicted water table from 60 observation well points.
This result follows K alluvium equal to 21.2 m/day and K fractured rock is equal to 0.12 m/day. These K numbers are related to CN alluvium and rock (i.e., 40 and 72, respectively), after being calculated using the equation of Van Mullem [21]. Furthermore, K alluvium is similar to the K median of constant head result from Al-Amri [13] (18.49 mm/day). The K for alluvium is distinguished into two directions, i.e., radial (Kr) and horizontal (Kz). The Kz equals 0.2 × Kr. This ratio agrees with the study by Todd [33]. On the other hand, in this study, Kr and Kz for fractured rock are the same (Kr = Kz).
The CHD value calibrated at the outlet is equal to 213 m asl, and the conductance is used for the DRN boundary package is equal to 100 m2/day.
All these results are validated using statistical analysis to see the goodness of the model. The RMSE is equal to 44.43 m, and R2 is equal to 0.9793. This RMSE value is considered minor because it is based on an extensive data range. The R2 close to one means a good relationship between the actual water table measurement and predicted water table data.

5.2.2. Transient Groundwater Modeling

The transient groundwater modeling using potential recharge from surface water modeling analysis is shown in Figure 21. The graph shows that the predicted line is far from the actual observation in Hadat Ash-Sham, and it is also below the lower percentile line.
The RMSE value of this modeling result is 3.09. This circumstance led to other assumptions that might be occurred in Al-Lusub watershed. This study applies two scenarios until the predicted water table line is similar to the actual average water table data and gets a smaller RMSE value. These scenarios are as follows:
  • Scenario 1:
Extra recharge from the neighboring watershed. This extra recharge flows through the geological structure (e.g., faults, fractures) that might have high permeability.
2.
Scenario 2:
Recharge from the ephemeral streams in Al-Lusub watershed. This scenario utilizes the RIV boundary package on MODFLOW 6.

Scenario 1: Extra Recharge from Neighboring Watershed

The geological analysis is performed to see the possible other water sources, which may be coming from the neighboring watershed. The geological fractures can be seen in Figure 22, based on the Saudi Arabia geological map for the Makkah region. The geological fractures are connecting the southern watershed with Al-Lusub watershed. This information shows that the additional recharge may be entering Al-Lusub watershed via geological structure (e.g., faults, fractures), which is highly permeable during storm events. The affected area is downstream of Al-Lusub watershed.
The mainstream (longest stream) is located near the downstream of Al-Lusub watershed. The direction of this stream is aimed to the north where the Al-Lusub watershed is located. This stream collects much water from the southern watershed; therefore, this stream might significantly influence additional recharge to Al-Lusub watershed during storm events.
Then, it is required to analyze the type of rainfall contour that will cause the southern watershed to deliver extra recharge into the Al-Lusub watershed. The hydrological analysis is performed to see the rainfall pattern in the region. Monthly rainfall data from stations in the region are utilized to develop the rainfall contour using the IDW method.
Examples in Figure 23 and Figure 24 show that there are two types of rainfall contour patterns that need to be considered:
  • Rainfall at upstream is higher than the rainfall at the surrounding area (Figure 23). If this circumstance occurred, then it is assumed that might be no extra recharge from southern watershed into Al-Lusub watershed. There are 8 months for this type of rainfall contour in this study (i.e., December 2019, February 2020, March 2020, July 2020, September 2020, December 2020, February 2021, and April 2021).
2.
Rainfall at upstream is lower than the rainfall at the surrounding area (Figure 24). This type of rainfall contour will possibly make extra recharge for the downstream Al-Lusub watershed. In this study, this type of contour happened in 13 months (i.e., November 2018, December 2018, January 2019, February 2019, May 2019, June 2019, August 2019, September 2019, October 2019, November 2019, December 2019, November 2020, and January 2021).
After proceeding with the calibration step, the extra recharge added downstream of Al-Lusub watershed is equal to five times the potential monthly recharge from Al-Lusub watershed. Figure 25 shows an additional RCH boundary package for downstream of the Al-Lusub watershed.
The result of Scenario 1 can be seen in Figure 26. It shows that the water table prediction line is close to the actual average water table and located between the upper and lower percentile lines. The RMSE value is 0.6 or lies between 0 and 1, indicating good model quality.
The yearly water balance for the entire system in Al-Lusub watershed, based on Scenario 1, is shown in Table 8. The cumulative water input and output for year 2019 are 105,402,119 m3 and 105,378,735 m3, while for 2020, they are 127,549,634 m3 and 127,525,717 m3, respectively. The percent discrepancy for both years is 0.02% or close to zero, which means that the model has good conservation of mass, and the error is minimum.

Scenario 2: Recharge from River, Using RIV Boundary Package on MODFLOW 6

This scenario is applied based on the assumption that the recharge downstream of Al-Lusub watershed is mainly from the streams. It is based on the alluvium appearance upstream of the Al-Lusub watershed, where the alluvium width is narrow. This narrow alluvium lies in the rocky mountainous areas; therefore, the role of upstream alluvium to increase the water table in the downstream area in Hadat Ash-Sham Farm is not significant. There is also a study by Budiman et al. [19] in Hadat Ash-Sham in which they concluded that the primary source for the downstream of Al-Lusub watershed is coming from the streams.
Then, it is required to calculate the stream water depth that lies in the alluvium area. The rainfall excess can be derived based on monthly rainfall and potential recharge data. Rainfall excess is multiplied by the watershed area (901.16 km2) to get the excess volume, and then divided by the alluvium area (130.80 km2). The result is the water depth for the alluvium area (Table 9).
Then the stream water depth is inputted into MODFLOW 6 using the RIV boundary package. The conductance value is 100 m2/day, the same as the conductance from steady-state groundwater modeling. The starting time is the same as the rainfall occurrences. There will be two parameters that will be selected after going through the calibration stage. These parameters are stream bottom and ending time. It is assumed that the value of stream bottom is equal to surface −1 m, then after the calibration process, the ending time for all periods is equal to 18 days from the starting time.
In this scenario, the DRN boundary package is not included. The RIV boundary package already includes the DRN function, i.e., to drain the groundwater system into the stream if the water table is above the water head on the stream. The MODFLOW 6 appearance for Scenario 2 can be seen in Figure 27, which involves several parameters.
The result of transient groundwater modeling for Scenario 2, involving the RIV boundary package on MODFLOW 6 can be seen in Figure 28. The graph shows that the modeled water table line is close to the average water table line in Hadat Ash-Sham Farm Station. The modeled water table is also located between the upper and lower percentile lines. The RMSE for this scenario is 0.70, which is positioned between 0 and 1; therefore, the modeling quality is good.
The yearly water balance for the entire system in the Al-Lusub watershed, based on Scenario 2, is shown in Table 10. The cumulative water input and output for year 2019 are 105,402,119 m3 and 105,378,735 m3, while for 2020 are 127,549,634 m3 and 127,525,717 m3, respectively. The percent discrepancy for both years is 0.02% or close to zero, which means that the model has good conservation of mass, and the error is minimum.

6. Conclusions and Recommendations

6.1. Conclusions

The CN values for surface water modeling of the basin have shown low values (for the alluvium and the fractured rock which are 40 and 72, respectively). These values interpret the significant losses from rainfall that transform to infiltration in arid environments.
The result of surface water modeling leads to a correlation between rainfall and potential recharge (potential recharge = 0.7986 rainfall). This result is confirmed by a recent study (about 77% of rainfall) [28].
The average K alluvium and K fractured rock values are 21.2 m/day and 0.12 m/day, respectively. The K alluvium is similar to the K obtained from the constant-head measurements conducted downstream of the Al-Lusub watershed. The value of K fractured rock is within the magnitude of the ranges available in the literature [32].
In steady-state groundwater modeling, the comparison of predicted and observed water table showed that the R2 = 0.98, which is classified as good; therefore, these K values are reliable.
In transient groundwater modeling using monthly potential recharge, it is concluded that two scenarios need to be addressed. The first scenario involves extra recharge from neighboring watersheds that potentially flow through the geological structure (i.e., faults, fractures). Based on the geological fractures, the extra recharge only influences the downstream area. The primary source of this extra recharge is from the mainstream of the southern watershed. The extra recharge from the neighboring watershed is only applied when the rainfall upstream of the Al-Lusub watershed is lower than the rainfall in the surrounding area. The water balance in the model is achieved with a relative error of 0.02%.
The second scenario assumes that the downstream area’s main recharge is from the ephemeral stream and the flow is running in the channel for 18 days. The mass balance error is also 0.02%, which is acceptable.
The calibration result for transient groundwater modeling shows that the RMSE for both scenarios is pretty good. Based on the water balance analysis, the percent discrepancy for both scenarios of transient groundwater is close to zero, which means that the model has good conservation of mass.

6.2. Recommendations

The first scenario of transient groundwater modeling shows the possibility of extra recharge that flows from neighboring watersheds into the downstream area of the Al-Lusub watershed. The primary source is from the mainstream of the southern watershed. To prove this assumption, it is recommended to perform a tracer test study or install new observation wells downstream of the Al-Lusub watershed and the southern watershed. The water flow movement can be identified from the water table contour of these wells.
The mainstream of Al-Lusub watershed is the primary source of groundwater recharge for the downstream area, where Hadat Ash-Sham Farm Station is located, according to the second scenario of transient groundwater modeling. Based on this information, it is essential to preserve natural stream conditions and develop artificial groundwater recharge schemes, so that the groundwater in the Hadat Ash-Sham region can be preserved.

Author Contributions

Conceptualization, N.A.-A., J.B. and A.E.; methodology, N.A.-A., J.B. and A.E.; software, J.B.; validation, J.B.; formal analysis, N.A.-A., J.B. and A.E.; investigation, N.A.-A., J.B. and A.E.; resources, N.A.-A., J.B. and A.E.; data curation, N.A.-A., J.B. and A.E.; writing—original draft preparation, J.B.; writing— review and editing, N.A.-A., J.B. and A.E.; visualization, J.B.; supervision, N.A.-A. and A.E.; project administration, N.A.-A. and A.E.; funding acquisition, N.A-A. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Data Availability Statement

Not applicable.

Acknowledgments

The authors would like to express gratitude to the staff of the Department of Hydrology and Water Resources Management, Faculty of Meteorology, Environment, and Arid Land Agriculture, King Abdulaziz University, who measured and collected data for this research.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Barthel, R. A call for more fundamental science in regional hydrogeology. Hydrogeol. J. 2014, 22, 507–510. [Google Scholar] [CrossRef]
  2. Fitts, C.R. Groundwater Science; Elsevier: Amsterdam, The Netherlands, 2002. [Google Scholar]
  3. Barthel, R. Common problematic aspects of coupling hydrological models with groundwater flow models on the river catchment scale. Adv. Geosci. 2006, 9, 63–71. [Google Scholar] [CrossRef]
  4. Scanlon, B.R.; Keese, K.E.; Flint, A.L.; Flint, L.E.; Gaye, C.B.; Edmunds, W.M.; Simmers, I. Global synthesis of groundwater recharge in semiarid and arid regions. Hydrol. Process. Int. J. 2006, 20, 3335–3370. [Google Scholar] [CrossRef]
  5. Barthel, R. HESS Opinions “Integration of groundwater and surface water research: An interdisciplinary problem?”. Hydrol. Earth Syst. Sci. 2014, 18, 2615–2628. [Google Scholar] [CrossRef]
  6. Wheater, H.S. Hydrological processes in arid and semi-arid areas. Hydrol. Wadi Syst. 2002, 5–22. [Google Scholar] [CrossRef]
  7. Ewea, H. Towards an Integrated Groundwater Management in Hadat Al Sham Area, Makkah Al-Mukarramah Region. J. King Abdulaziz Univ. 2011, 22, 99–116. [Google Scholar] [CrossRef]
  8. Charlton, R. Fundamentals of Fluvial Geomorphology; Routledge: London, UK, 2007. [Google Scholar]
  9. Horton, R.E. Drainage-basin characteristics. Eos Trans. Am. Geophys. Union 1932, 13, 350–361. [Google Scholar] [CrossRef]
  10. Strahler, A.N. Quantitative analysis of watershed geomorphology. Eos. Trans. Am. Geophys. Union 1957, 38, 913–920. [Google Scholar] [CrossRef]
  11. El-Hames, A. Determination of groundwater availability in shallow arid region aquifers utilizing GIS technology: A case study in Hada Al-Sham, Western Saudi Arabia. Hydrogeol. J. 2005, 13, 640–648. [Google Scholar] [CrossRef]
  12. Sharaf, M.; Al-Bassam, A.; Bayumi, T.; Qari, M. Hydrogeological and Hydrochemical Investigation of the Cretaceous-Quaternary Sedimentary Sequence East of Jeddah City; King Abdulaziz City for Science and Technology: Riyadh, Saudi Arabia, 2002. [Google Scholar]
  13. Al-Amri, N.S. Hydrogeological Assesment of the Haddat Ash Sham Region; University of Birmingham: Jeddah, Saudi Arabia, 2005. [Google Scholar]
  14. Alyamani, M.S.; Şen, Z. Determination of hydraulic conductivity from complete grain-size distribution curves. Groundwater 1993, 31, 551–555. [Google Scholar] [CrossRef]
  15. Al-Khatib, E. Hydrogeology of Usfan District. A Thesis of Master of Science. Master’s Thesis, Institute of Applied Geology, King Abdulaziz University, Jeddah, Saudi Arabia, 1977. [Google Scholar]
  16. Hussein, M.; Bazuhair, A. Groundwater in Haddat Al Sham-Al Bayada area, western Saudi Arabia. Arab Gulf J. Sci. Res. 1992, 10, 23–43. [Google Scholar]
  17. Gelhar, L.; LW, G. Effects of hydraulic conductivity variations on groundwater flows. In Proceedings of the Second International Symposium on Stochastic Hydraulics, Lund, Sweden, 2—4 August 1976; pp. 409–431. [Google Scholar]
  18. Al-Juhani, A. Evaluation of Groundwater Resources at Hadat Ash Sham Area: Wadi Ghulah (Western Saudi Arabia). Master’s Thesis, Faculty of Earth Sciences, King Abdulaziz University, Jeddah, Saudi Arabia, 2002. [Google Scholar]
  19. Budiman, J.S.; Al-Amri, N.S.; Chaabani, A.; Elfeki, A.M. Geostatistical based framework for spatial modeling of groundwater level during dry and wet seasons in an arid region: A case study at Hadat Ash-Sham experimental station, Saudi Arabia. Stoch. Environ. Res. Risk Assess. 2021, 36, 2085–2099. [Google Scholar] [CrossRef]
  20. Kent, K.M. Travel time, time of concentration and lag. In National Engineering Handbook; NRCS: Washington, DC, USA, 1972; Volume 4. [Google Scholar]
  21. Van Mullem, J.A. Applications of the Green-Ampt Infiltration Model to Watersheds in Montana and Wyoming; College of Engineering, Montana State University-Bozeman: Bozeman, MT, USA, 1989. [Google Scholar]
  22. Elfeki, A.; Masoud, M.; Basahi, J.; Zaidi, S. A unified approach for hydrological modeling of arid catchments for flood hazards assessment: Case study of wadi Itwad, southwest of Saudi Arabia. Arab. J. Geosci. 2020, 13, 490. [Google Scholar] [CrossRef]
  23. Elfeki, A.M.; Ewea, H.A.; Al-Amri, N.S. Development of storm hyetographs for flood forecasting in the Kingdom of Saudi Arabia. Arab. J. Geosci. 2014, 7, 4387–4398. [Google Scholar] [CrossRef]
  24. Hawkins, R.H. The importance of accurate curve numbers in the estimation of storm runoff 1. JAWRA J. Am. Water Resour. Assoc. 1975, 11, 887–890. [Google Scholar] [CrossRef]
  25. Farran, M.M.; Elfeki, A.M. Statistical analysis of NRCS curve number (NRCS-CN) in arid basins based on historical data. Arab. J. Geosci. 2020, 13, 31. [Google Scholar] [CrossRef]
  26. Yuan, Y.; Nie, W.; McCutcheon, S.C.; Taguas, E.V. Initial abstraction and curve numbers for semiarid watersheds in Southeastern Arizona. Hydrol. Process. 2014, 28, 774–783. [Google Scholar] [CrossRef]
  27. Farran, M.M.; Elfeki, A.; Elhag, M.; Chaabani, A. A comparative study of the estimation methods for NRCS curve number of natural arid basins and the impact on flash flood predications. Arab. J. Geosci. 2021, 14, 1–23. [Google Scholar] [CrossRef]
  28. Farran, M.M.; Al-Amri, N.; Elfeki, A.M. Aquifer recharge from flash floods in the arid environment: A mass balance approach at the catchment scale. Hydrol. Process. 2021, 35, e14318. [Google Scholar] [CrossRef]
  29. Albishi, M.; Bahrawi, J.; Elfeki, A. Empirical equations for flood analysis in arid zones: The Ari-Zo model. Arab. J. Geosci. 2017, 10, 51. [Google Scholar] [CrossRef]
  30. Batelaan, O.; Smedt, F.D.; Becker, P.D.; Huybrechts, W. Characterization of a Regional Groundwater Discharge Area by Combined Analysis of Hydrochemistry, Remote Sensing and Groundwater Modelling; Balkema: Rotterdam, The Netherlands, 1998. [Google Scholar]
  31. Harbaugh, A.W.; Banta, E.R.; Hill, M.C.; McDonald, M.G. Modflow-2000, the U.S. Geological Survey Modular Ground-Water Model-User Guide to Modularization Concepts and the Ground-Water Flow Process; Open-File Report; U.S. Geological Survey: Reston, VA, USA, 2000; p. 134. [Google Scholar]
  32. Freeze, R.A.; Cherry, J.A. Groundwater; Prentice-Hall Inc.: Englewood Cliffs, NJ, USA, 1979; Volume 7632, p. 604. [Google Scholar]
  33. Todd, D. Groundwater Hydrology; Jon Wiley Sons Inc.: New York, NY, USA, 1980; p. 535. [Google Scholar]
  34. Chai, T.; Draxler, R.R. Root mean square error (RMSE) or mean absolute error (MAE)?—Arguments against avoiding RMSE in the literature. Geosci. Model Dev. 2014, 7, 1247–1250. [Google Scholar] [CrossRef] [Green Version]
  35. Johnston, K.; Verhoef, J.; Krivoruchko, K.; Lucas, N. ArcGIS Geostatistical Analyst Tutorial. USA ESRI 2003, 256–258. Available online: https://dusk.geo.orst.edu/gis/geostat_analyst.pdf (accessed on 2 September 2022).
  36. McDonald, M.G.; Harbaugh, A.W. The history of MODFLOW. Groundwater 2003, 41, 280–283. [Google Scholar] [CrossRef] [PubMed]
  37. Lohman, S.W. Ground-Water Hydraulics; US Government Printing Office: Washington, DC, USA, 1972; Volume 708. [Google Scholar]
Figure 1. Schemes follow the same formatting.
Figure 1. Schemes follow the same formatting.
Water 14 03075 g001
Figure 2. Al-Lusub watershed drainage configuration according to digital elevation model.
Figure 2. Al-Lusub watershed drainage configuration according to digital elevation model.
Water 14 03075 g002
Figure 3. Rainfall stations in Al-Lusub watershed overlayed with Thiessen Polygon.
Figure 3. Rainfall stations in Al-Lusub watershed overlayed with Thiessen Polygon.
Water 14 03075 g003
Figure 4. Rainfall data for station J-104 and J-214: (A) maximum daily rainfall for each month; (B) monthly rainfall.
Figure 4. Rainfall data for station J-104 and J-214: (A) maximum daily rainfall for each month; (B) monthly rainfall.
Water 14 03075 g004
Figure 5. Well location based on El-Hames (2005) study. The data contain water table elevation and bottom of shallow aquifer information.
Figure 5. Well location based on El-Hames (2005) study. The data contain water table elevation and bottom of shallow aquifer information.
Water 14 03075 g005
Figure 6. Location of Hadat Ash-Sham Farm Station.
Figure 6. Location of Hadat Ash-Sham Farm Station.
Water 14 03075 g006
Figure 7. (A) Water table data from 11 wells in Hadat Ash-Sham Farm Station; (B) The monthly average, upper percentile (95%), and lower percentile (5%) of water table data from 11 wells, Hadat Ash-Sham Farm Station.
Figure 7. (A) Water table data from 11 wells in Hadat Ash-Sham Farm Station; (B) The monthly average, upper percentile (95%), and lower percentile (5%) of water table data from 11 wells, Hadat Ash-Sham Farm Station.
Water 14 03075 g007
Figure 8. Conceptual model of the integrated surface and groundwater modeling: (A) The watershed contains alluvium and fractured rock on the surface; (B) Cross-section from a–b, where the water table elevation is getting higher to the upstream; (C) Cross-section from c–d, where during rainy season water table is high; (D) Cross-section from c–d, where during rainy season water table is low.
Figure 8. Conceptual model of the integrated surface and groundwater modeling: (A) The watershed contains alluvium and fractured rock on the surface; (B) Cross-section from a–b, where the water table elevation is getting higher to the upstream; (C) Cross-section from c–d, where during rainy season water table is high; (D) Cross-section from c–d, where during rainy season water table is low.
Water 14 03075 g008
Figure 9. The methodology flow chart in this study.
Figure 9. The methodology flow chart in this study.
Water 14 03075 g009
Figure 10. Dimensionless cumulative hyetograph, first quartile for Jeddah City (Elfeki et al., 2014).
Figure 10. Dimensionless cumulative hyetograph, first quartile for Jeddah City (Elfeki et al., 2014).
Water 14 03075 g010
Figure 11. Alluvium and rock delineation in Al-Lusub watershed.
Figure 11. Alluvium and rock delineation in Al-Lusub watershed.
Water 14 03075 g011
Figure 12. Grid construction for Al-Lusub watershed. The cell size is 200 m × 200 m.
Figure 12. Grid construction for Al-Lusub watershed. The cell size is 200 m × 200 m.
Water 14 03075 g012
Figure 13. Shallow aquifer depth distribution overlayed with satellite image.
Figure 13. Shallow aquifer depth distribution overlayed with satellite image.
Water 14 03075 g013
Figure 14. Steady-state groundwater modeling appearances on MODFLOW 6, involve several boundary conditions packages (i.e., CHD and DRN).
Figure 14. Steady-state groundwater modeling appearances on MODFLOW 6, involve several boundary conditions packages (i.e., CHD and DRN).
Water 14 03075 g014
Figure 15. OBS (observation utility package) is located on existed observation well.
Figure 15. OBS (observation utility package) is located on existed observation well.
Water 14 03075 g015
Figure 16. Transient groundwater modeling appearances on MODFLOW 6, involve several boundary conditions packages (i.e., CHD, DRN, EVT, and RCH). The EVT and RCH are applied to the whole watershed.
Figure 16. Transient groundwater modeling appearances on MODFLOW 6, involve several boundary conditions packages (i.e., CHD, DRN, EVT, and RCH). The EVT and RCH are applied to the whole watershed.
Water 14 03075 g016
Figure 17. OBS (observation utility package) during transient groundwater modeling stage, that located at Hadat Ash-Sham Farm Station.
Figure 17. OBS (observation utility package) during transient groundwater modeling stage, that located at Hadat Ash-Sham Farm Station.
Water 14 03075 g017
Figure 18. Hyetograph and hydrograph for maximum daily rainfall in August, September, October, and November 2019.
Figure 18. Hyetograph and hydrograph for maximum daily rainfall in August, September, October, and November 2019.
Water 14 03075 g018
Figure 19. Graph of rainfall vs. potential recharge using CN 67.36 based on maximum daily rainfall on each month.
Figure 19. Graph of rainfall vs. potential recharge using CN 67.36 based on maximum daily rainfall on each month.
Water 14 03075 g019
Figure 20. Scatter plot of actual water table measurement and predicted water table.
Figure 20. Scatter plot of actual water table measurement and predicted water table.
Water 14 03075 g020
Figure 21. Monthly rainfall and water table time series as the result of transient groundwater modeling using potential recharge from monthly rainfall.
Figure 21. Monthly rainfall and water table time series as the result of transient groundwater modeling using potential recharge from monthly rainfall.
Water 14 03075 g021
Figure 22. Geological fractures in the study area overlayed with the streams Al-Lusub watershed and southern watershed.
Figure 22. Geological fractures in the study area overlayed with the streams Al-Lusub watershed and southern watershed.
Water 14 03075 g022
Figure 23. The rainfall contours projected over the fractures map to show the case of no extra recharge from the southern watershed into Al-Lusub watershed.
Figure 23. The rainfall contours projected over the fractures map to show the case of no extra recharge from the southern watershed into Al-Lusub watershed.
Water 14 03075 g023
Figure 24. The rainfall contours projected over the fractures map to show the potential extra recharge from the southern watershed into Al-Lusub watershed in the downstream area.
Figure 24. The rainfall contours projected over the fractures map to show the potential extra recharge from the southern watershed into Al-Lusub watershed in the downstream area.
Water 14 03075 g024
Figure 25. Transient groundwater modeling appearances on MODFLOW 6 for scenario 1, involving several boundary conditions packages (i.e., CHD, DRN, EVT, and RCH). The EVT and RCH are applied to the whole watershed, and the RCH for extra recharge is applied for the downstream of the Al-Lusub watershed.
Figure 25. Transient groundwater modeling appearances on MODFLOW 6 for scenario 1, involving several boundary conditions packages (i.e., CHD, DRN, EVT, and RCH). The EVT and RCH are applied to the whole watershed, and the RCH for extra recharge is applied for the downstream of the Al-Lusub watershed.
Water 14 03075 g025
Figure 26. Monthly rainfall and water table time series as the result of transient groundwater modeling for Scenario 1 (extra recharge from neighboring watershed).
Figure 26. Monthly rainfall and water table time series as the result of transient groundwater modeling for Scenario 1 (extra recharge from neighboring watershed).
Water 14 03075 g026
Figure 27. Transient groundwater modeling appearances on MODFLOW 6 in Scenario 2, involve several boundary conditions packages (i.e., CHD, EVT, RCH, and RIV). The EVT and RCH are applied to the whole watershed.
Figure 27. Transient groundwater modeling appearances on MODFLOW 6 in Scenario 2, involve several boundary conditions packages (i.e., CHD, EVT, RCH, and RIV). The EVT and RCH are applied to the whole watershed.
Water 14 03075 g027
Figure 28. Monthly rainfall and water table time-series results of transient groundwater modeling for Scenario 2 (discharge from the ephemeral stream in Al-Lusub watershed).
Figure 28. Monthly rainfall and water table time-series results of transient groundwater modeling for Scenario 2 (discharge from the ephemeral stream in Al-Lusub watershed).
Water 14 03075 g028
Table 1. Morphometric information of Al-Lusub watershed.
Table 1. Morphometric information of Al-Lusub watershed.
NoParametersValue
1Area901.16 km2
2Perimeter270 km
3Basin Length72.02 km
4Basin Slope0.1063 m/m
5Total Stream Length710.5 km
6Maximum Stream Length96.80 km
7Maximum Stream Slope0.0109 m/m
8Sinuosity1.34
9Drainage Density0.0008
10Hypsometry Integral0.48
Table 2. The hydraulic conductivity measurement by Al-Amri [13] using grain-size analysis and constant-head method.
Table 2. The hydraulic conductivity measurement by Al-Amri [13] using grain-size analysis and constant-head method.
K Grain Size Analysis (m/day)K Constant-Head (m/day)
Min0.954.25
Max35.1755.32
Mean10.6219.48
Median7.7218.49
Table 3. Hydraulic conductivity data for the alluvium in Al-Lusub watershed.
Table 3. Hydraulic conductivity data for the alluvium in Al-Lusub watershed.
Hydraulic Conductivity (m/day)Source
4.25–55.32Al-Amri [13]
86.4–216El-Hames [11]
43.2–69.12Al-Juhani [18]
95.04Sharaf et al. [12]
5.6-62.6Hussein and Bazuhair [16]
103.68Al-Khatib [15]
Table 4. Hydraulic conductivity of the alluvium and fractured rock based on the CN value.
Table 4. Hydraulic conductivity of the alluvium and fractured rock based on the CN value.
CNK (m/Day)
Alluvium4021.2145
Fractured Rock720.1243
Table 5. Predicted CN value in Al-Lusub watershed.
Table 5. Predicted CN value in Al-Lusub watershed.
CN Alluvium CN Fractured Rock CN AverageInitial Abstraction (mm)
407267.361.23
Table 6. Tc and Tlag value of Al-Lusub watershed.
Table 6. Tc and Tlag value of Al-Lusub watershed.
Basin Length (km)Basin Slope (m/m)Tc (h)T lag, 0.6Tc (h)
72.020.10631.881.12
Table 7. Potential recharge based on monthly rainfall data.
Table 7. Potential recharge based on monthly rainfall data.
MonthMonthly Rainfall Data (mm)Total Monthly Rainfall (mm)Potential Monthly Recharge (mm)
Station J-104Station J-214
Nov-1831.4-7.856.87
Dec-180.2-0.050.05
Jan-1912.2-3.053.03
Feb-190.8-0.20.2
Mar-19----
Apr-19----
May-1937.253.249.239.89
Jun-191.4-0.350.35
Jul-19----
Aug-1928.43.49.658.30
Sep-199.815.41411.78
Oct-1918.414.215.2512.77
Nov-19-1.41.051.05
Dec-19610.69.458.14
Jan-20----
Feb-200.62.62.12.1
Mar-200.61.81.51.5
Apr-20----
May-20----
Jul-19----
Jul-200.41.41.151.15
Aug-20----
Sep-202285.669.756.26
Oct-2025.245.640.532.94
Nov-202.6-0.650.65
Dec-2012.637.631.3525.63
Jan-211.6-0.40.4
Feb-214.42116.8514.05
Mar-21----
Apr-21-75.254.79
Table 8. Yearly water balance for transient groundwater modeling Scenario 1.
Table 8. Yearly water balance for transient groundwater modeling Scenario 1.
YearTotal in (m3)Total out (m3)Percent Discrepancy
2019105,402,119105,378,7350.02%
2020127,549,634127,525,7170.02%
Table 9. Estimested ponding water depth on the stream (alluvium).
Table 9. Estimested ponding water depth on the stream (alluvium).
MonthMonthly Rainfall (mm)Monthly Potential Recharge (mm)Rainfall Excess (mm)Excess Volume (m3)Ponding Water Depth (m)
Nov-187.856.870.98883,484.010.00675
Dec-180.050.050.0000
Jan-193.053.040.0112,318.800.00009
Feb-190.20.2000
Mar-19-----
Apr-19-----
May-1949.239.899.318,388,209.300.06413
Jun-190.350.35000
Jul-19-----
Aug-199.658.311.341,210,170.960.00925
Sep-191411.782.221,999,664.430.01529
Oct-1915.2512.782.472,226,530.370.01702
Nov-191.051.05000
Dec-19-----
Jan-209.458.151.301,173,872.410.00897
Feb-202.12.10000
Mar-20-----
Apr-20-----
May-20-----
Jun-201.51.5000
Jul-20-----
Aug-201.151.15000
Sep-2069.756.2613.4412,108,810.710.09257
Oct-2040.532.947.566,809,222.360.05206
Nov-200.650.650.000.000.00000
Dec-2031.3525.645.715,148,563.680.03936
Jan-210.40.400.000.000.00000
Feb-2116.8514.062.792,516,918.770.01924
Mar-21-----
Apr-215.254.790.46411,602.850.00315
Table 10. Yearly water balance for transient groundwater modeling Scenario 2.
Table 10. Yearly water balance for transient groundwater modeling Scenario 2.
YearTotal in (m3)Total out (m3)Percent Discrepancy
2019129,548,305.92129,517,995.900.02%
2020164,043,765.75164,008,080.530.02%
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Al-Amri, N.; Budiman, J.; Elfeki, A. Integrated Surface Water and Groundwater Modeling in Arid Environment, Al-Lusub Watershed, Saudi Arabia. Water 2022, 14, 3075. https://doi.org/10.3390/w14193075

AMA Style

Al-Amri N, Budiman J, Elfeki A. Integrated Surface Water and Groundwater Modeling in Arid Environment, Al-Lusub Watershed, Saudi Arabia. Water. 2022; 14(19):3075. https://doi.org/10.3390/w14193075

Chicago/Turabian Style

Al-Amri, Nassir, Jaka Budiman, and Amro Elfeki. 2022. "Integrated Surface Water and Groundwater Modeling in Arid Environment, Al-Lusub Watershed, Saudi Arabia" Water 14, no. 19: 3075. https://doi.org/10.3390/w14193075

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop