The GRACE satellites have been monitoring the Earth’s gravitational pull since their 2002 launch to detect variations in mass. GRACE is comprised of a pair of satellites that orbit the Earth at 450–500 km apart and share an identical orbit. When Earth’s gravity field shifts, the satellites’ velocities adjust accordingly. As the pair of satellites reaches a positive mass anomaly, such a mountainous region, the gravitational pull increases, causing the leading satellite to speed up and increase the gap between them (Houborg et al., 2012). It has been found that variations in TWS are the primary cause of mass anomalies in the absence of tectonic activity (Swenson et al., 2008).

This research made use of three different time-variable gravity mascon solutions from the GRACE processing center. Mascon products include datasets from the Jet Propulsion Laboratory (JPL), the Center for Space Research (CSR) at University of Texas, and the Goddard Space Flight Center (GSFC) at NASA all comprise gravity solutions. When compared to spherical harmonic solutions, the spatial resolution and accuracy of these mascon solutions are superior; they are able to detect all signals below the GRACE noise threshold. Not even a little de-striping or smoothing is necessary. These products may not even need the scaling factor (Watkins et al., 2015; Luthcke et al., 2013; Save et al., 2016; Wiese et al., 2016).

In the JPL-RL06M v02 solution, the TWS anomalies were determined for each of the three spherical cap mascon blocks with the same area (Watkins et al., 2015; Wiese et al., 2016). The TWS fluctuations are sampled in the final data set to 0.5° × 0.5° grids. By employing the scaling factors provided by JPL, leakage signals caused by 3° mascon blocks can be recovered, presuming the JPL-RL06M v02 has a real 3° resolution. The new JPL-RL06M v02 data improves upon the original by allowing for greater differentiation between land and sea mascons (JPL-RL05M). For each mascon that crosses a coastline, the Coastal Resolution Improvement filter was applied as a last step in the process (Watkins et al., 2015).

The CSR RL06M v02 solution that was computed from an equal region at a resolution of 1, displays TWS variations in 0.25° × 0.25° grids (Save et al., 2016; Sultan and Save, 2020). This new grid separates the shoreline’s hexagonal tiles into ocean and land sections to reduce interference between the two. The scaling factor is unnecessary for such tiny mascons. We also utilized the mascon solution produced by the GSFC-RL06v1.0 for each 0.5° equal-area square mascon.

The missing monthly data was interpolated using the cubic-spine approach. In this work, the mascon solutions were used significantly. We determined the slope values of the TWS trends. In the next step, we estimated the errors in the trend values.

The choice of using GRACE mascon data in this research is well-justified based on its unique capabilities and advantages. GRACE, consisting of two satellites orbiting the Earth, is designed to monitor variations in Earth’s gravitational field caused by changes in mass. These variations can result from factors like total water storage (TWS) changes, making GRACE a valuable tool for assessing groundwater fluctuations.

Since data from gauge stations in the research area were unavailable or insufficient during the time period under study, we rely on estimates of soil moisture, extracted from the Global Land Data Assimilation System (GLDAS) model (Rodell et al., 2004). Those results were put to use in order to break down the TWS into its constituent parts. This is because GRACE does not understand the various TWS parts. In this research, we compared two different NSAS GLDAS versions (VIC and CLM). The same time period as GRACE data was covered by the monthly data, which has a spatial resolution of one degree. The use of GLDAS data can be justified due to the absence of gauge station data, the need to dissect TWS into its components, the reliability of GLDAS, the appropriate spatial resolution, and the temporal consistency with GRACE data. This choice enables a more comprehensive and accurate analysis of groundwater dynamics and their contributing factors in the study area.

Due to a lack of available rainfall data from the ground, the TRMM (Kummerow et al., 1998) data was used instead. The time span covered by TRMM data, which has a spatial resolution of 0.25° × 0.25°, overlaps with that of GRACE and GLDAS. Through yearly averaging of all monthly rainfall, TRMMdata was utilized to create a map of the average annual rainfall (AAR) (Figure 3).

Based on the weather patterns and the beginning of the 2007 drought, two separate time periods have been determined (Trigo et al., 2010). The first period, from 04/2002 to 12/2006, experienced a higher than average rate of precipitation, whereas the second period, from January 2007 to July 2017, experienced a lower than average rate of precipitation. Each period’s AAP rate was determined by averaging that period’s monthly rainfall data over a full year. For this study, we used the AAP to look into how precipitation affected GWS in various time periods. The use of TRMM data not only fills a critical data gap resulting from the lack of ground-based rainfall measurements but also provides valuable information on long-term precipitation trends. This data source enabled to explore the connections between precipitation and groundwater storage variations over specific time periods, which is essential for the study’s objectives and analyses.

The period of study in this investigation spans from the launch of the GRACE satellites in 2002 to the conclusion of data analysis in 2022. The integration of GRACE, GLDAS, and TRMM data sources, along with various processing steps, allowed for a comprehensive examination of groundwater dynamics and their relationship to precipitation patterns over these years. ArcGIS 10.6.1 was used to prepare all the maps in this study.

Figure 2 presents a graph displaying the monthly rainfall rates, revealing a pattern of higher values from December to May and lower values from June to November. This climatic data facilitates the identification of three distinct time periods for the entire Makkah region, attributed to the substantial area with varying precipitation levels. Precipitation diminishes as we move from the coastal zone towards the central part of Saudi Arabia, as depicted in Figure 3. These three time intervals were determined based on the typical annual rainfall pattern, as shown in Figure 4. The Average Annual Precipitation (AAP) ranged from a minimum of 59.7 mm in 2008 to a maximum of 157.8 mm in 2018. Period I (2002–2011) had the lowest AAP at 92.8 mm, while period III (2017–2019) recorded the highest AAP at 101.8 mm. Period II (2012–2016), falling between the other two, had an AAP of 48.7 mm (refer to Table 1). Focusing on Wadi Fatimah and its nearby coastal area, the AAP exhibited slightly higher values, ranging from 130.7 to 131.1 mm throughout the study period, compared to the lower AAP of 99.3 mm for the entire Makkah region. Given that Wadi Fatma and its surrounding area represent the lowlands near the coast within the study area, the AAP for Wadi Fatma displayed a slightly higher positive trend, calculated at 3.20 mm. This is in contrast to the lower positive trend of 1.07 mm observed over the entire region during the entire study period. These precipitation rates continue to underscore the arid nature of the study area, with most rainfall events and storms occurring between the months of November to May.

Time-series plots depicting the variations in ΔTWS across the study area are displayed in Figure 5, combining data from two different solutions and their average. Notably, a strong correlation of 0.82 was observed between the two mascon datasets, encompassing the entire region. Analyzing the average trends in TWS from CSR and GSFC, positive values are evident in the southwestern portions of the research area, while negative trends prevail in the northeastern regions (Figure 6). Conducting a linear regression analysis on the averaged ΔTWS solutions illustrated in Figure 5 reveals three distinct trends across the entirety of the Makkah region throughout the study period (Figure 7; Table 1). The first trend (period I) exhibits a slightly negative slope, estimated at −0.49 ± 0.051 cm/yr (Table 1). Period II displays a very slight positive trend, estimated at −0.012 ± 0.115 cm/yr. The third trend (period III) takes a notably steep descent, with a slope of −0.64 ± 0.15 cm/yr. Interestingly, these trend variations do not align with the Average Annual Precipitation (AAP) time series, which displays a positive increase from 92.8 mm in Period I to 116.4 mm in Period III (Table 1; Figure 4).

Shifting the focus to Wadi Fatma, the temporal patterns in GRACE-derived TWS reveal slightly negative trends, estimated at −0.18 ± 0.039, −0.19 ± 0.019, and −0.18 ± 0.023 cm/yr, as determined by the mascon solutions from CSR, GSFC, and their averages (Table 2; Figure 7). The smaller area surrounding Wadi Fatma exhibits TWS patterns that align somewhat with those of Wadi Fatma (Figure 5). The TWS extracted through GRACE illustrates trend values estimated at −0.19 ± 0.036, −0.19 ± 0.017, and −0.19 ± 0.022 cm/yr when considering the GAREC mascon CSR, GSFC, and their averaging solutions (Table 2; Figure 7). As depicted in Figure 5, the seasonal fluctuations of TWS for Wadi Fatma exhibit relatively higher values from October to April and lower values from May to September. Importantly, these TWS peaks are noticeably delayed by 1–3 months compared to the peaks in rainfall. This temporal discrepancy may be attributed to the gradual response of TWS to changes in precipitation, potentially influenced by surface runoff within the study area. Additionally, the presence of surface basement rocks might hinder the infiltration of water, partly due to the moderately steep slope of the ground surface. Analyzing the long-term trends in TWS variations in Wadi Fatma reveals negative TWS values at the upstream part of the Wadi, gradually decreasing to slightly negative values as we move downstream (Figure 6).

Because GRACE lacks the vertical resolution to distinguish the various components of total water storage (TWS), such as surface water storage (SWS), groundwater storage (GWS), and soil moisture storage (SMS), it becomes challenging to attribute mass changes in a given geographical area to specific storage types. This limitation arises from the fact that the gravitational signal detected by GRACE results from the combined impact of mass changes in the atmosphere, oceans, and terrestrial water storage, as noted by Wahr et al. (2006). To address this challenge, we employed a range of methodologies, including the utilization of climate models, to identify and differentiate the contributions of each storage component. While it is worth noting that traditional soil moisture measurements have often been characterized by a lack of reliability and consistency, recent advances in both satellite and ground-based observational data products have significantly improved the accuracy of soil moisture (SM) data.

One valuable system in this regard is the Global Land Data Assimilation System (GLDAS), as developed by Rodell et al., in 2004 and 2009, as well as Moiwo et al. (2011). GLDAS integrates various sources of data to provide comprehensive information about land surface conditions, including soil moisture. These datasets play a crucial role in diverse applications, including climate prediction, hydrological studies, and water resource management.

To assess changes in SMS throughout the study period, we leveraged three GLDAS datasets, namely, VIC, CLM, and Noah. Despite differences in SMS trends, the variations in SMS ranged from +0.008 ± 0.002 mm/yr for the entire Makkah region to 0.046 ± 0.026 mm/yr for the smaller area surrounding the Wadi. This variation in SMS trends is influenced by factors such as local rainfall patterns. In the case of Wadi Fatma, characterized by slightly higher rainfall data and its position downstream from surrounding regions, the SMS exhibited a relatively higher trend, calculated at +0.41 ± 0.1 mm/yr over the entire study period.

The capabilities of GRACE to distinguish between different components of total water storage (TWS), including soil moisture, surface water, and groundwater, have limitations. To address this challenge, we employed GLDAS model results to estimate the non-groundwater components, as illustrated in Figure 8. Changes in surface water storage (SWS) were found to be insignificant due to the lower rates of surface runoff.

We utilized Eq. 1 to compute groundwater storage (GWS) by subtracting the GLDAS-derived soil moisture storage (SMS) from the GRACE-derived total water storage (TWS). Figure 9 illustrates the variation in GWS, which, much like TWS, exhibits higher values between October and April and lower values from May to September. Δ TWS = Δ GWS + Δ SMS (1)

The time series of ΔGWS, derived using Eq 1, is presented in Figure 10. The analysis of this figure reveals that the average GWS across the entire Makkah region is declining at a rate of −0.38 ± 0.017 cm/year throughout the study period. Period I displays a negative GWS trend, estimated at −0.49 ± 0.049 cm/yr. In Period II, the GWS trend remains relatively stable, computed at −0.0172 ± 0.11 cm/yr, while Period III indicates a significant decrease in GWS, with a trend of −0.607 ± 0.15 cm/yr (Table 1; Figure 10).

For Wadi Fatma, employing the adjacent mascon solutions from GARCE CSR and GSFC, we determined a slightly negative trend in GRACE-derived GWS at −0.18 ± 0.022 cm/yr. The GWS patterns in the surrounding area closely resemble those of Wadi Fatma. Using this averaged mascon solution yielded a trend similar to that of Wadi Fatma, calculated at −0.185 ± 0.022 cm/yr (Table 2; Figure 10). Notably, Wadi Fatma experiences a slower rate of groundwater depletion compared to the entire Makkah region, attributed to its proximity to the coast and slightly higher precipitation rates. Furthermore, the structural influence of the southwest-trending fault zone plays a role in shaping the behavior of groundwater storage in Wadi Fatma, distinguishing it from the broader Makkah region.

We obtained data on sediment thickness from the website of the National Oceanic and Atmospheric Administration’s National Geophysical Data Center (Divins, 2003). The sedimentary layers exhibit variations in thickness within the Wadi area. The eastern part of the Wadi has sediments, approximately 100 m thick, whereas the central part has sediment thickness exceeding 700 m. As we move towards the western part of the Wadi, the sediment thickness reduces to around 600 m. In a small section in the southwestern part, notably, there are elevated values of sediment thickness, reaching up to 1800 m (as shown in Figure 11).

The significant sediment thickness observed in the central and southwestern portions of the Wadi, along with the slightly negative total water storage (TWS) signals in this area, suggests the presence of a substantial groundwater reservoir. The eastern part of the Wadi receives substantial precipitation, which is likely to drain into the central and downstream regions, contributing to the recharge of groundwater aquifers. However, it is worth noting that groundwater levels are declining at a rate of −0.18 ± 0.022 cm per year over the entire study period, primarily due to the heavy exploitation of groundwater resulting from the growing population and increased agricultural activities.

The topographical information provided by Figure 12, which depicts the elevation variations in the Wadi Fatimah basin, can be crucial in understanding and integrating water storage studies using GRACE data.The elevation differences observed in the upstream section, ranging from 1000 to 2000 m above sea level, are indicative of the landscape’s topographic complexity. In areas with significant elevation changes, water storage dynamics can exhibit unique patterns and responses to various hydrological factors.

Furthermore, the presence of rugged igneous and metamorphic rock formations, as well as the network of narrow faulted wadis or channels, can influence groundwater flow and storage in the region. These geological features may affect the permeability of subsurface layers, impacting the movement and retention of water. The small wetland areas, responsible for channeling surface water towards the central and southwestern parts of the region, are essential components of the local hydrological system. These areas can act as natural conveyors of water, potentially contributing to groundwater recharge in the adjacent aquifers.

The periodic runoff events and the saturation of alluvial deposits above the weathered basement following rainfall occurrences are critical aspects of the water cycle in the region. These processes can significantly influence water storage variations, particularly in the middle and downstream segments of the Wadi.