Influence of Land Use and Land Cover Changes and Precipitation Patterns on Groundwater Storage in the Mississippi River Watershed: Insights from GRACE Satellite Data

Abstract

Growing human demands are placing significant pressure on groundwater resources, causing declines in many regions. Identifying areas where groundwater levels are declining due to human activities is essential for effective resource management. This study investigates the influence of land use and land cover, crop types, and precipitation patterns on groundwater level trends across the Mississippi River Watershed (MRW), USA. Groundwater storage changes from 2003 to 2015 were estimated using data from the Gravity Recovery and Climate Experiment (GRACE) satellite mission. A spatiotemporal analysis was conducted at four scales: the entire MRW, groundwater regimes based on groundwater level change rates, 31 states within the MRW, and six USGS hydrologic unit code (HUC)-2 watersheds. The results indicate that the Lower Mississippi region experienced the fastest groundwater decline, with a Sen’s slope of −0.07 cm/year for the mean equivalent water thickness, which was attributed to intensive groundwater-based soybean farming. By comparing groundwater levels with changes in land use, crop types, and precipitation, trends driven by human activities were identified. This work underscores the ongoing relevance of GRACE data and the GRACE Follow-On mission, launched in 2018, which continues to provide vital data for monitoring groundwater storage. These insights are critical for managing groundwater resources and mitigating human impacts on the environment.

1. Introduction

Groundwater accounts for up to 30% of the fresh water available on earth for human use. Over the past few decades, dependence on groundwater sources for daily activities has increased significantly. Consequently, the estimation and monitoring of groundwater resources have become critical due to the growing pressure to meet water demands. While developed nations have monitoring well networks in place, determining groundwater storage on a regional scale using a limited number of wells fails to provide comprehensive data as groundwater occurrence is heterogeneous both spatially and temporally [1,2]. This challenge is further complicated by data processing inconsistencies, missing spatio-temporal data, human method-produced and systematic errors, and a lack of metadata for converting groundwater head measurements to groundwater storage [3]. To bridge this data gap, researchers have begun developing and testing satellite-based remote sensing systems to capture the spatio-temporal variations in groundwater storage [4]. One such satellite-based system is the Gravity Recovery and Climate Experiment (GRACE; [5]).

GRACE is a twin satellite-based remote sensing system that detects gravitational anomalies near the Earth’s surface while orbiting in a near-polar trajectory at an altitude of approximately 485 km. Jointly managed by NASA and the German Aerospace Centre, the GRACE mission was launched in 2002 and has provided valuable measurements of terrestrial water storage anomalies worldwide. After the original mission concluded in 2017, the GRACE Follow-On (GRACE-FO) mission was launched in May 2018 to continue this critical work. GRACE-FO builds on the legacy of its predecessor, continuing to monitor Earth’s gravitational field anomalies and offering ongoing insights into changes in groundwater storage, ice sheets, and other key environmental indicators. Numerous studies have demonstrated that both GRACE and GRACE-FO are valuable for understanding hydrology and the global water cycle [4,6,7,8,9,10]. Monthly terrestrial water storage anomalies derived from GRACE can be decomposed into various components. By isolating the contributions of individual components to temporal mass variability through auxiliary observations and numerical models, it is possible to estimate changes in groundwater storage over large regions [11].

GRACE missions have demonstrated great promise for estimating groundwater anomalies in various regions [3,12,13,14,15,16,17,18,19,20]. For instance, Rodell et al. [3] used GRACE data from April 2002 to July 2004 to estimate groundwater storage changes in the Mississippi River Watershed, demonstrating the efficacy of GRACE in capturing groundwater storage variability. Further investigations into the spatial variation in groundwater storage, such as those by Tiwari et al. [21] and Rodell et al. [12], revealed large-scale mass loss in North India due to excessive groundwater extraction. Additionally, stress indicator maps of groundwater levels in the Upper Ganges were generated, showing a decline in groundwater levels observed through a combination of numerical hydrological models and nationwide water use statistics [22].

These studies effectively depicted the net rates of groundwater decline but often did not fully explore the underlying reasons for fluctuations. Moreover, most previous research focused on single spatial-scale analyses, which may not provide a comprehensive understanding of temporal or spatial variability. Analyzing groundwater trends at multiple spatial scales is essential to capture the spatial and temporal heterogeneity of groundwater resources, revealing localized patterns that may be overlooked in broader analyses. A multi-scale approach offers more detailed insights into the influence of local factors, making the findings more actionable for both regional and localized groundwater management.

This study utilized GRACE data to analyze groundwater trends within the Mississippi River Watershed (MRW) through a comprehensive “four-scale spatiotemporal analysis”. GRACE data from 2003 to 2015 were examined to identify trends in groundwater levels alongside trends in land use and land cover changes, precipitation patterns, and cropping systems. This analysis was conducted at four distinct scales: the macro scale (entire MRW), the meso scale (specific groundwater regimes based on the rate of groundwater level change), the regional scale (United States Geological Survey (USGS) hydrologic unit code (HUC)-2 watersheds), and the local scale (state level). This study evaluated how land use and land cover changes, precipitation, and cropping patterns influence groundwater levels, revealing the factors driving these trends across the four scales.

2. Study Area

The study area for this research is the Mississippi River Watershed (MRW) (Figure 1). The Mississippi River, the largest river in the United States, serves as a critical water source for irrigation, urban use, the transportation of goods, and other human activities in the USA. The MRW is the fourth largest watershed globally, covering 40% of the contiguous United States and supporting some of the nation’s most important agricultural regions. The MRW spans 31 states, including Alabama, Arkansas, Colorado, Georgia, Illinois, Indiana, Iowa, Kansas, Kentucky, Louisiana, Maryland, Michigan, Minnesota, Mississippi, Missouri, Montana, Nebraska, New Mexico, New York, North Carolina, North Dakota, Ohio, Oklahoma, Pennsylvania, South Dakota, Texas, Tennessee, Virginia, West Virginia, Wyoming, and Wisconsin. It also extends into two Canadian provinces, Alberta and Saskatchewan.

The United States Geological Survey (USGS) divides the United States into numerically coded hydrologic unit code (HUC) watersheds, which are part of a four-level hierarchical hydrologic unit system. This system classifies geographic areas based on the drainage basins of large rivers or aggregated drainage areas. HUC-2 watersheds, represented by two-digit codes, are the largest hydrologic regions. The MRW encompasses six USGS HUC-2 watersheds: Ohio (05), Tennessee (06), Upper Mississippi (07), Lower Mississippi (08), Missouri (10), and Arkansas-White-Red (11).

3. Materials and Methods

GRACE-based estimates of terrestrial water storage (ΔTWS) were obtained for the period of 2003–2015 from GRACE-Tellus (Land release RL05 version), the data archive of the GRACE mission managed by NASA. The spatial resolution of the obtained data was 1° × 1°, and the temporal resolution was one month. The GRACE satellite system measures the gravitational anomalies and records these observations as mass distributions near the Earth’s surface. The ΔTWS consists of contributions from soil moisture (ΔSM), snow water equivalent (ΔSWE), groundwater storage equivalent (ΔGWE), canopy water equivalent (ΔCWE), and other sources. However, the contribution from other sources is negligible compared to the ΔSM, ΔSWE, ΔCWE, and ΔGWE. Therefore, groundwater storage can be numerically modeled using the following equation:

ΔGWE = ΔTWS − ΔSM − ΔSWE − ΔCWE

(1)

Terrestrial water storage was obtained from GRACE data, and associated scaling factors were applied. The components subtracted from the ΔTWS were derived from Global Land Data Assimilation Systems (GLDAS), which provides land surface parameters simulated by land surface models. The NOAH land surface model, with a spatial resolution of 1° × 1° and a temporal resolution of one month, was used to obtain these components. Additionally, monthly well log data from 500 wells across the entire MRW for the period of 2003–2015 were obtained from the USGS and used to validate the groundwater trends identified from the GRACE data.

The methodology adopted in this study is illustrated in the flowchart shown in Figure 2. The mean monthly groundwater anomaly image, calculated for the entire MRW from 155 images obtained between 2003 and 2015 using Equation (1), was classified using the Natural Jenks algorithm to determine groundwater regimes. The Natural Jenks algorithm was chosen because it identifies natural groupings inherent within the dataset and optimally divides the data by minimizing variance within classes while maximizing the deviation between class means. Through this method, the MRW was segmented into three distinct groundwater regimes based on annual groundwater level changes (Figure 3): (i) a regime with an increasing trend in groundwater levels, (ii) a regime with a stable groundwater level, and (iii) a regime with a declining trend in groundwater levels.

Trends for the states, the HUC-2 watersheds, and the regimes were derived using the Mann–Kendall trend test, and the direction of trends was determined through Sen’s slope estimates. The Mann–Kendall test is a non-parametric method used to identify trends in a time series. If the p-value from the test was less than the significant value of 0.05 at 95% confidence interval, the null hypothesis was rejected, indicating that a significant monotonic trend existed. The tau value from the Mann–Kendall test indicated whether the trend was increasing or decreasing; a positive tau value signified an increasing trend, while a negative value indicated a decreasing trend. Sen’s estimator of slope, a non-parametric robust method, was used to calculate the slope of the trend by determining the median of all lines passing through data points in a series. A positive slope indicated an increasing trend, and a negative slope indicated a decreasing trend.

Annual land use and land cover (LULC) data with a spatial resolution of 30 m were obtained from the United States Department of Agriculture (USDA) NRCS Geospatial Data Gateway, specifically the USDA Cropland Data Layer, to assess the influence of LULC trends on groundwater levels in the watersheds. Different classification schemes were used for the LULC data from 2008 to 2010 and from 2011 to 2015, requiring a combination of some classes in the data from the 2008–2010 period to align with the classification scheme used for the data from the 2011–2015 period. The USDA Cropland Data Layer LULC datasets included 255 distinct classes. For each year, these classes were consolidated into major LULC categories following the classification scheme proposed by Anderson et al. [23]. Subsequently, the proportions of LULC types—agriculture, forests, rangelands, and urban areas—were calculated as fractions of the total area covered by each LULC type.

To calculate the factors affecting groundwater trends across four scales—(1) the entire MRW watershed, (2) the three groundwater regimes, (3) the 31 states, and (4) the six HUC-2 watersheds—an integrated set of tasks was followed. LULC trends were derived by calculating how the fractional land coverage of specific LULC types increased or decreased over time. The fraction of land covered by a particular LULC type was determined by calculating the proportion of pixels occupied by that LULC in the raster data for each year, allowing for the creation of a time series over the specified period.

In addition to LULC, the groundwater levels in a region are also influenced by local climatic conditions. To evaluate the impact of precipitation on groundwater levels, trends were analyzed across the four scales alongside the spatiotemporal distribution of LULC and groundwater. Monthly precipitation data for this study was derived from the NOAH land surface model in GLDAS with a spatial resolution of 0.125° × 0.125°. Furthermore, this study examined the influence of cropping patterns on groundwater levels to determine whether specific crop types affected the rate of groundwater abstraction.

4. Results

The results and insights derived from the trend analyses of groundwater levels, LULC, precipitation, and cropping patterns across the four aforementioned scales are presented below.

4.1. Groundwater Trends

The groundwater anomalies for each month were used to plot the time series of the mean equivalent water thickness together with the LULC classes of agriculture, forests, rangeland, and urban areas for the first scale, i.e., the entire MRW (Figure 4). The mean equivalent water thickness in each month indicates changes in groundwater levels with respect to the previous month. Positive values represent an increase in groundwater levels, while negative values indicate a decrease. The error bars in the plot represent one standard deviation from the data over the period of 2003–2015. Simultaneously, the LULC trends for agriculture, forests, rangeland, and urban areas were plotted with the groundwater trends (Figure 4). Across the MRW, forests and rangelands showed a declining trend, while agriculture and urban areas exhibited an increasing trend, with the exception of anomalous behavior in urban areas during the 2008–2010 period. As stated earlier, this exception is attributed to differences in classification schemes used in the data. These trends suggest that since 2008, there has been a marked increase in urbanization and the conversion of forests and rangelands into agricultural and urban areas across the MRW.

The loss of forests can significantly impact groundwater levels, as demonstrated in other studies using GRACE data [24]. This creates a form of a positive feedback loop: the reduction in forest cover, combined with increased groundwater extraction for water-intensive crops like soybeans, can severely deplete groundwater storage over time. Additionally, a sharp decline in the mean equivalent water thickness was observed across the MRW in the 2012–2013 period, aligning with the fact that 2012 was a year of prolonged drought in the United States [25]

For the second scale, the states within the MRW with the highest proportion of land covered by (i) agriculture, (ii) forests, (iii) rangeland, and (iv) urban areas were identified using the 2015 LULC raster data. Iowa had the largest fraction of its land dedicated to agriculture, West Virginia to forests, Wyoming to rangelands, and Ohio to urban areas. Time series plots of the mean equivalent water thickness were generated for four states within the MRW, each chosen for having the largest area covered by a specific LULC type among the 31 states (Figure 5).

A time series is said to be weakly stationary if it meets the following conditions: (i) the mean remains almost constant over time, indicating no trend in the series; (ii) the variance remains stable, indicating no seasonality; and (iii) the covariance and correlation between the time series and its inherent lags remain almost constant. The time series for states with maximum forest cover (West Virginia), rangeland cover (Wyoming), and urban area coverage (Ohio) were weakly stationary. However, Iowa, with its extensive agricultural land, exhibited an increasing trend with a positive tau value of 0.3391 and a Sen’s slope estimate of 0.0692 (

Table 1. Mann–Kendall test results and Sen’s slope estimates in cm/year for mean equivalent water thickness for states in Mississippi River Watershed (no shading: decreasing trend; light grey shade: trend is not significant; dark grey shade: increasing trend). States are sorted as per Sen’s slope estimates.

State	Mann–Kendall Test p-Value	Mann–Kendall Test Tau Value	Sen’s Slope Estimate
Mississippi	0.0000664	−0.21609	−0.1053
Louisiana	0.000000936	−0.26569	−0.07216
Arkansas	0.0000376	−0.22329	−0.05038
Oklahoma	0.000000671	−0.26921	−0.04676
Colorado	0.002327	−0.16498	−0.01302
Tennessee	0.453942	−0.04064	−0.00782
Kentucky	0.970379	−0.00209	−0.00028
North Dakota	0.5967	0.028739	0.002648
Kansas	0.356451	0.050021	0.004662
West Virginia	0.381176	0.047507	0.005164
Indiana	0.174315	0.073649	0.009286
Missouri	0.254803	0.061751	0.009759
Ohio	0.023127	0.123083	0.016796
Wyoming	0.004164	0.155258	0.017525
Montana	0.000277	0.196984	0.021292
Pennsylvania	0.001361	0.173523	0.02424
Illinois	0.000769	0.182237	0.028055
South Dakota	0.0000000356	0.298534	0.042364
Wisconsin	0.0000254	0.228152	0.043546
Minnesota	0.00024	0.198995	0.045663
Nebraska	0.00000011	0.287641	0.050749
Iowa	0.000000000385	0.339087	0.069212

). In 2012, the groundwater levels in Iowa declined by 15 cm, likely due to excessive groundwater withdrawal for agriculture. Iowa, being heavily agriculture-dependent, was among the hardest-hit states during the 2012 drought [25].

For Tennessee, Kentucky, North Dakota, Kansas, West Virginia, Indiana, and Missouri, the p-values from the Mann–Kendall trend test were greater than 0.05, indicating no significant monotonic trend in the time series (Table 1). However, the p-values and tau values for Mississippi, Louisiana, Arkansas, Oklahoma, and Colorado suggest a declining trend in their respective time series. The rate of decline was the most pronounced for Mississippi, followed by Louisiana and Arkansas with slopes of −0.1053, −0.07216, and −0.05038, respectively (Table 1).

Similarly, the HUC-2 watersheds (third scale) within the MRW with the highest proportion of land covered by (i) agriculture, (ii) forests, (iii) rangeland, and (iv) urban areas were identified. The generated time series plots exhibited similar characteristics of weak stationarity as observed previously (Figure 6). According to the results from the Mann–Kendall trend tests, no significant trend was detected in the time series for the Tennessee and Ohio HUC-2 watersheds (

Table 2. Mann–Kendall test results and Sen’s slope estimate cm/year for mean equivalent water thickness for HUC-2 watersheds in Mississippi River Watershed (no shading: decreasing trend; light grey shade: trend is not significant; dark grey shade: increasing trend). HUC-2 Watersheds are sorted as per Sen’s slope estimates.

HUC-2 Watersheds	Hydrologic Region	Mann–Kendall Test p-Value	Mann–Kendall Test Tau Value	Sen’s Slope Estimate
08	Lower Mississippi	0.000118356	−0.208546292	−0.07067408
06	Tennessee	0.430061102	−0.042815249	−0.00734257
05	Ohio	0.474720228	0.038793465	0.006503669
10	Missouri	0.0000807	0.213573523	0.024593365
07	Upper Mississippi	0.00000034	0.276246334	0.046603873

). In contrast, an increasing trend was observed for the Upper Mississippi and Missouri watersheds (Table 2). For the Tennessee HUC-2 watershed, a decreasing trend was identified with a Sen’s slope of −0.007 cm/year for the mean equivalent water thickness, although the p-value was greater than 0.05, indicating that the trend was not statistically significant. For the Lower Mississippi HUC-2 watershed, a significant declining trend was observed, with the highest rate of groundwater depletion, reflected by a Sen’s slope of −0.07 cm/year for the mean equivalent water thickness.

Among the groundwater regimes, Region 1 exhibited a slightly positive rate for groundwater level change, Region 2 showed an almost constant rate, and Region 3 experienced a significantly negative rate of change between 2003 and 2015 (Figure 7). Interestingly, both Regions 1 and 3 with a slightly increasing trend and a pronounced declining trend, respectively, were marked by a rise in agricultural land use. This contrasting behavior illustrates how an increase in agricultural areas can lead to differing groundwater trends across regions. In some areas, agriculture promotes groundwater recharge through deeper infiltration, while in others, intensive irrigation and water-intensive crops can lead to the over-extraction and depletion of groundwater resources. Therefore, the effects of agricultural expansion on groundwater trends are highly dependent on local environmental conditions and water management practices.

The states and HUC-2 watersheds with the fastest and slowest rates of groundwater change were also identified (Figure 8). Among the 31 states in the MRW, the state of Mississippi exhibited the highest rate of groundwater change (Figure 8a), while Kentucky had the lowest (Figure 8b). At the HUC-2 watershed level, the Lower Mississippi HUC-2 watershed experienced the fastest rate of groundwater change (Figure 8c), whereas the Ohio River HUC-2 watershed showed the slowest (Figure 8d). Changes in agricultural practices in the Lower Mississippi HUC-2 watershed have contributed to a significant decline in groundwater levels, which was observed by interpreting the LULC trends.

4.2. Land Use and Land Cover Trends

Land use and land cover (LULC) trends for the groundwater regimes are shown in Figure 9. After ignoring the initial inconsistencies due to different classification schemes, it was observed that Region 3, which exhibited a high rate of groundwater level change, showed increasing trends in both urban areas and agricultural land. This suggests that the expansion in farming and urbanization, coupled with the decline in rangeland coverage in this region, contributed to the decrease in groundwater levels.

LULC trends were also analyzed at the state and HUC-2 watershed levels for areas with the fastest and slowest rates of groundwater level change (Figure 10). Additionally, LULC trends were derived for the states and HUC-2 watersheds with the highest and lowest groundwater level changes per unit area. The groundwater level change per unit area was calculated by normalizing the number of pixels occupied by each LULC class with the total number of pixels. The names of the identified states and HUC-2 watersheds are given in

Table 3. The states and HUC-2 watersheds with the fastest and the slowest groundwater change.

	Fastest Rate of Groundwater Change	Slowest Rate of Groundwater Change	Highest Groundwater Per Pixel	Lowest Groundwater Per Pixel
States	Mississippi	Kentucky	South Dakota	Mississippi
HUC-2 watersheds	Lower Mississippi	Ohio	Missouri	Lower Mississippi

. In the Mississippi region (the state of Mississippi and the lower Mississippi HUC-2 watershed), groundwater levels are among the lowest and are declining at a much faster rate compared to other areas. The fraction of rangeland, which naturally promotes rainwater infiltration, is also decreasing in this region. In contrast, the Missouri watershed has the highest fraction of rangeland and, consequently, the highest groundwater per pixel value. Furthermore, the expansion of urban and agricultural areas is placing additional pressure on local groundwater resources, as indicated by the negative correlation between the groundwater anomalies and both urban and agricultural area cover.

4.3. Precipitation Trends

The precipitation trends in the groundwater regimes were investigated, with Figure 11 depicting trends for the regions of our interest: Regions 1 and 3. Overall, no significant trends were observed in precipitation for these; however, seasonality is present in the precipitation time series. Precipitation trends were also analyzed at the state and HUC-2 watershed levels (Figure 12). States and HUC-2 watersheds with the largest areal coverage of agriculture, forests, rangelands, and urban areas were identified. The time series of precipitation rates for these areas suggest that areas with maximum forest and urban coverage exhibit a declining precipitation trend (Figure 12). Recent studies indicate that vegetation, particularly tree cover, has a greater influence on rainfall patterns than commonly assumed [26]. Over the past decade, forest cover has declined across the MRW, which has, in turn, impacted the region’s climate and rainfall patterns. Since rainfall is somewhat dependent on forest cover, areas with extensive forests are experiencing this dual effect of declining tree cover and reduced precipitation. In urban areas, studies show that anthropogenic activities can lead to regional precipitation declines [27]. These human activities alter greenhouse gas and ozone levels, changing the overall atmospheric composition, which indirectly affects atmospheric water vapor levels and, consequently, rainfall.

4.4. Cropping Pattern Trends

The distribution of the top ten major crops in 2015 within the MRW was determined based on their areal coverage, i.e., the number of pixels representing each crop type (Figure 13). Time series analyses were conducted for the top five crops at the entire MRW level (Figure 14) as well as at the state and HUC-2 watershed levels (Figure 15 and Figure 16). Since these time series were derived from LULC data, the analysis period was selected as 2010–2015. At the watershed level, the five major crops analyzed were corn, soybeans, winter wheat, alfalfa, and hay/non-alfalfa. The variation in the fraction of the area covered by these crops relative to the total area is shown in Figure 14. A significant drop in crop production observed between 2012 and 2015 was largely attributed to the 2012 drought [25]. Some crops, such as corn, experienced a continuous decline in production, while others saw a sharp decrease followed by recovery to normal levels (Figure 14).

Agriculture was identified as the dominant LULC class in states such as Arkansas, Indiana, Iowa, Louisiana, Mississippi, and Minnesota. Similarly, in the HUC-2 watershed level, the Upper Mississippi, Lower Mississippi, Missouri River, and Arkansas-White-Red River watersheds were found to have agriculture as their major LULC class. These states and HUC-2 watersheds, identified from the 2015 LULC raster data, had more than 50% of their land devoted to agriculture. Time series plots for the top three crop types (winter wheat, soybeans, and corn) were generated for these states and HUC-2 watersheds and are shown in Figure 13 and Figure 14, respectively. Crop production across the entire MRW was significantly impacted by the 2012 drought. However, soybean production in the Lower Mississippi region continued to increase despite the drought. Similar crop trends were observed in the HUC-2 watersheds. An increasing trend for soybeans is clear and eminently visible in the crop patterns of the Lower Mississippi watershed and exhibited a clear negative correlation with groundwater levels during the study period (Figure 8c). This supports the earlier conclusion that groundwater levels are influenced not only by LULC and precipitation but also by the types of crops grown, particularly when agriculture is the dominant LULC class.

4.5. USGS Well Log Data

Time series plots were generated for the mean depth to the water table for each well. An inverse relationship was observed between changes in groundwater levels (derived from GRACE data) and the mean depth to water below the surface (obtained from well log data). As groundwater levels decline, the depth of the water table increases. This relationship was evident in approximately 70% of the wells across the region, providing validation for the GRACE data against observed water levels in the wells. Validation plots at the state level are shown in Figure 17 and at the HUC-2 watershed level in Figure 18. Although the validation plots demonstrate a reasonable level of accuracy for GRACE data (70%), some errors and uncertainties persist, largely due to spatiotemporal under-sampling of data points. While unbiased errors tend to average out across the entire region, spatial under-sampling remains a concern, especially given the significant heterogeneity of aquifers, even at smaller scales [3].

5. Discussion and Conclusions

This study provides a comprehensive, watershed-wide analysis of groundwater changes in a large watershed and investigates the various factors contributing to these changes. Miller et al. [28] examined the impact of urbanization on storm runoff in a peri-urban catchment and found a 400% increase in flow due to urbanization. Similarly, Schilling et al. [29] quantified the effects of land use and land and cover (LULC) changes on discharge in the Upper Mississippi River and discovered that increased soybean acreage led to higher discharge. Our study adds a unique perspective by specifically evaluating the influence of urbanization and changes in cropping patterns on groundwater levels in the MRW. It fills a critical knowledge gap by examining how LULC affects groundwater levels at multiple spatial scales. In MRW, groundwater levels are declining due to anthropogenic activities. Urbanization has significantly expanded in the Lower Mississippi watershed over the past decade (Figure 10e), and much of the region’s agricultural production, predominantly soybeans, relies heavily on groundwater. This rising demand has negatively impacted groundwater levels. Our findings validate this, as soybean production trends show a strong negative correlation with groundwater anomalies in the region (Figure 8c and Figure 16c). The USDA National Agricultural Statistics Service also reported a steady increase in soybean production from 2008 to 2017, both in Mississippi and across the nation [30]. However, groundwater levels in the upper parts of the MRW were not observed to be declining. States like Indiana and Iowa, where agriculture is the dominant LULC type (Figure 5a), were somewhat affected, but overall, groundwater trends in these regions remain slightly positive. During the drought year of 2012, many anthropogenic activities increased reliance on groundwater sources, which explains the significant decline in groundwater levels across the entire watershed from 2012 to 2014.

There have been numerous investigations in the field of hydrogeology and groundwater studies, but most focus primarily on the first-order effects of groundwater depletion in a region ([31,32,33]). This study introduces a four-tier approach for analyzing groundwater trends using terrestrial water storage data derived from GRACE, along with associated changes in land use and land cover (LULC), precipitation, and cropping patterns, to explain the groundwater trends observed at four spatial scales. Although the boundaries of watersheds and subsurface aquifer systems do not always align, given the low spatial resolution of GRACE data and the large scale of the Mississippi River Watershed (MRW), the use of surface-water watersheds in this study is appropriate. Furthermore, leakage error is a type of signal distortion that can occur in GRACE data when filtering time-variable gravity field models. Applying appropriate scaling factors helps minimize these errors. However, some leakage may remain uncorrected by these scaling factors; this residual leakage can still impact GRACE-based estimates. We plan to quantify these uncertainties in our future work. Nevertheless, the close correlation between groundwater estimates from GRACE and observed well levels further validates the robustness of our approach. This multi-scale study draws meaningful correlations between each scale, and the regional analysis serves to validate the broader findings. For future studies, other climatic factors such as air temperature should also be considered, as previous research has shown that temperature can significantly influence groundwater levels [7,34,35].

Although many developed nations have long monitored groundwater levels through well networks, GRACE measurements are invaluable for monitoring and estimating groundwater levels at larger spatial scales. This method’s full potential can be realized particularly in data-deficient regions of the world or areas where access to data is limited. The methodology in this study incorporates data only up to 2015 due to the degradation in data quality following the failure of one of GRACE’s twin satellites’ accelerometer. Despite this setback, GRACE continued operating until the end of 2017. The GRACE Follow-On (GRACE-FO) mission, a joint effort by NASA and the Helmholtz Centre Potsdam German Research Centre for Geosciences (GFZ), now continues GRACE’s legacy [6]. The approaches used in this study can be easily adapted for GRACE-FO data, which would greatly benefit data-limited regions by enabling the analysis of groundwater trends. These findings can then inform models designed to manage land use in watersheds and help mitigate groundwater depletion.