Irrigation-Intensive Groundwater Modeling of Complex Aquifer Systems Through Integration of Big Geological Data

Vahdat-Aboueshagh, Hamid; Tsai, Frank T.-C.; Bhatta, Dependra; Paudel, Krishna P.

doi:10.3389/frwa.2021.623476

ORIGINAL RESEARCH article

Front. Water, 06 April 2021

Sec. Water and Hydrocomplexity

Volume 3 - 2021 | https://doi.org/10.3389/frwa.2021.623476

This article is part of the Research Topic Modeling-based Approaches for Water Resources Problems View all 7 articles

Irrigation-Intensive Groundwater Modeling of Complex Aquifer Systems Through Integration of Big Geological Data

$\nHamid Vahdat-Aboueshagh$ Hamid Vahdat-Aboueshagh¹

Frank T.-C. Tsai²^*

Dependra Bhatta³

Krishna P. Paudel³

¹Department of Civil and Environmental Engineering, Louisiana State University, Baton Rouge, LA, United States
²Department of Civil and Environmental Engineering, Patrick F. Taylor Hall, Louisiana State University, Baton Rouge, LA, United States
³Department of Agricultural Economics and Agribusiness, Louisiana State University and LSU Agricultural Center, Baton Rouge, LA, United States

This study identifies hydrogeologic characteristics of complex aquifers based on constructing stratigraphic structure with large, non-uniform well log data. The approach was validated through a modeling study of the irrigation-intensive Chicot aquifer system, which is an important Pleistocene-Holocene aquifer of the Coastal Lowlands aquifer system in the southwestern Louisiana. Various well log types were unified into the same data structure, prioritized based on data sources, and interpolated to generate a detailed stratigraphic structure. More than 29,000 well logs were integrated to construct a stratigraphy model of 56 model layers for the Chicot aquifer system. The stratigraphy model revealed interconnections of various sands in the system, where 90% of the model domain is covered by fine-grained sediments. Although the groundwater model estimated a slight groundwater storage gain during 2005–2014 for the entire region, groundwater storage in the agricultural area was depleted. Nevertheless, the quick groundwater storage recovery during the non-irrigation seasons suggests that the Chicot aquifer system is a prolific aquifer system. The groundwater modeling result shows that the gulfward groundwater flow direction prior to pumping has been reversed toward inland pumping areas. The large upward vertical flow from the deeper sands indicates potential saltwater migration from the base of the Chicot aquifer system.

Introduction

Groundwater resource management of coastal aquifers has become a focal center of interest for many studies (Bocanegra et al., 2010; Custodio, 2010; Christelis and Mantoglou, 2016; Quintana et al., 2018). Overdrafting coastal aquifers has reportedly caused significant change in pre-development flow patterns leading to saltwater intrusion (Barlow and Reichard, 2010; Azizi et al., 2019; Parizi et al., 2019) and land subsidence (Eggleston and Pope, 2013; Guo et al., 2015). Hence, sustainable management of the aquifers plays an important role in counteracting the negative impacts (Abd-Elhamid and Javadi, 2011; Sreekanth and Datta, 2015). A key to understand the sustainability of aquifer systems is the analysis of influxes and outfluxes, which is elaborated in the concept of sustainable yield (Bredehoeft, 2002; Alley and Leake, 2004). As previous researches have emphasized, measuring the contributions of different inputs and outputs to a groundwater system is not a straightforward task for complex aquifer systems and flow modeling could help to overcome this issue (Loucks, 2000; Cao et al., 2013). The Chicot aquifer system in the southwestern Louisiana, as a very important groundwater resource in the Coastal Lowlands aquifer system (Weiss, 1990), has been exploited mostly for irrigation during the past decades (Sargent, 2011). Understanding groundwater dynamics in the Chicot aquifer system is the primary step for its future sustainable development.

Understanding the groundwater dynamics in complex aquifers needs the assistance of stratigraphy models as well as groundwater models with different hydrogeologic properties (Zhang et al., 2006; Lee et al., 2007; Berg and Illman, 2015). Selecting an appropriate stratigraphic modeling method mainly depends on data availability, computational resources, and adequacy of results reflecting the reality. Geostatistical methods have been extensively applied to developing stratigraphy models (Johnson and Dreiss, 1989; Soares, 1990; Marinoni, 2003; Kostic et al., 2005). However, computational expenses with a very large dataset are among the limitations of using kriging methods (Pardo-Iguzquiza and Dowd, 1998; Tartakovsky et al., 2007; Cheng, 2013). Non-geostatistical interpolation methods such as the natural neighbor interpolation method (Watson, 1999) and the triangulated irregular networks (Clevis et al., 2006) avoid covariance matrix calculations and are computationally efficient.

One of the challenges associated with constructing stratigraphic structure for large-scale models is processing a huge volume of geological data. Several attempts have been made in organizing well log records and geological information by using Geographic Information System (GIS) tools (Chang and Park, 2004; Ross et al., 2005; Chesnaux et al., 2011). However, how to process raw well log data has not received enough attention. Prioritizing data from different well log sources based on their reliability (Ross et al., 2005) and selecting the most reliable piece of information are solutions to reduce errors in stratigraphy models. Another way to increase the model reliability is to consider an effective depth, to which a given well log includes proper stratigraphic information (Danielsen et al., 2007).

Model calibration is very challenging for an expensive groundwater model. Due to high dimensionality of parameters in a large-scale, complicated aquifer system, it is impossible to tune parameters without the aid of automatic calibration through optimization. Global search optimization methods such as the genetic algorithm (Kumar et al., 2010), the simulated annealing algorithm (Aarts and van Laarhoven, 1987), the ant colony optimization (Dorigo and Di Caro, 1999), and the particle swarm optimization (Poli et al., 2007) have the advantage of being derivative-free (Rios and Sahinidis, 2013) over derivative-based local search optimization methods, but require a large number of model runs. The Covariance Matrix Adaptation Evolutionary Strategy (CMA-ES) (Hansen et al., 2003) has proven to be an efficient method to calibrate expensive groundwater models with parallel computing (Elshall et al., 2015). This study adopted the CMA-ES for model calibration.

The goal of this study is to address the two aforementioned modeling challenges: (1) building a complex stratigraphy structure with large, non-uniform well log data, and (2) calibration of a high dimensional groundwater model, while investigating the groundwater dynamics of the Chicot aquifer system in Southwest Louisiana. The Chicot aquifer system has been known for its intensive use of groundwater for irrigation (Konikow, 2013) and its complex aquifer structure is largely ignored. This study adopted a lithofacies modeling method and a parallel optimization method for model calibration to improve groundwater modeling. Moreover, to the authors' best knowledge, there is no such a detailed large-scale groundwater model been built to advance hydrogeological knowledge in the Gulf of Mexico region and other coastal areas. The study presents an important worked example to bridging the gap between research and professional practice. In section Chicot Aquifer System, the Chicot aquifer system is introduced. Section Data and Methods discusses data and methods to develop a Chicot groundwater model. Section Results and Discussions shows results and discussion. Section Conclusions draws the conclusions.

Chicot Aquifer System

The Chicot aquifer system is part of the Coastal Lowlands aquifer system (Weiss, 1990). It covers Southeast Texas and Southwest Louisiana (Weiss, 1990). Nominally, the Chicot aquifer system only refers to its extent in Southwest Louisiana (Jones et al., 1956; Nyman et al., 1990; Sargent, 2004). The annual precipitation is about 1,473 mm (58 inches). The Chicot aquifer system is the most heavily pumped aquifer in Louisiana, accounting for 41% of the state's groundwater use (Sargent, 2011). Seventy-two percentage of groundwater withdrawals are for the agriculture sector (including aquaculture). The Chicot aquifer system is also a sole source aquifer (USEPA, 2019), providing the major water supply to the public in the region. Figure 1 shows the study domain, major rivers, lakes, the Chenier Plain, and 2,709 pumping wells. The Chicot aquifer boundaries are determined by the USGS (Smoot, 1986). About 87% of total pumping wells have depth <120 m, which are generally irrigation wells. Deep wells are usually public supply wells and industrial wells. More than 2,400 irrigation wells are in the east-central area. Monthly pumping rates for irrigation were estimated based on cultivation area, crop type, growing seasons of crops, and pumping capacity of wells. The pumping wells in the west-central area are mostly industrial wells. The major public supply wells are located in or around municipalities. The pumping rates for the industrial use, public supply, and power generation were obtained from the U.S. Geological Survey (USGS).

FIGURE 1

Figure 1. Chicot aquifer domain (red line). Black dots are pumping wells. Blue lines show the major rivers in the Chicot aquifer domain. The pie chart shows the fractions for depth ranges of pumping wells. The coordinates are in NAD83 / UTM zone 15N.

The Chicot aquifer system is the uppermost hydrogeological unit in southwest Louisiana and encompasses seven aquifers (Nyman et al., 1990). The aquifer system gently dips southward at an angle of 0.06° (Jones et al., 1956). Figure 2 displays the location of well logs and geological fault traces. The fault traces are part of the Tepetate fault system (McCulloh and Heinrich, 2013). Well log data include 29,107 drillers' logs and 811 electrical logs that were collected from the Louisiana Department of Natural Resources (LDNR) and the USGS. Well logs are sparse in the Chenier Plain at the south and in the Atchafalaya River Basin at the east. The drillers' logs are mostly shallow (<120 m), accounting for 93% of total well logs. The electrical logs are primarily deep. Seventy percentage of electrical logs have depth deeper than 120 m. The well logs show about 40% of lithological information to be sand facies, indicating a sizable amount of permeable formations for storing groundwater.

FIGURE 2

Figure 2. Location of well logs and fault traces in the study area. The pie chart shows the fractions for depth ranges of well logs.

Data and Methods

Well Logs for Stratigraphy Model Development

This study employed two types of well logs to develop a stratigraphy model: electrical logs and drillers' logs. Electrical logs show the electrical properties of strata along a borehole, where in general high resistivity and low resistivity indicate sand facies (high-permeability facies) and clay facies (low-permeability facies), respectively (Archie, 1942; Galloway, 1977). The gamma ray and spontaneous potential curves when available were also used to assist sand picks in the electrical logs. The 811 electrical logs were analyzed. Each electrical log was interpreted as a sequence of sand facies and clay facies. Since electrical logs usually start at a deeper depth below land surface and contain stratigraphic information for deeper formations, they were employed to build the deeper portion of the stratigraphy model.

The drillers' logs contain useful lithologic descriptions along boreholes. Since drillers' logs normally begin at land surface, they provide important lithologic information for the top portion of the stratigraphy model. Lithologic descriptions in all 29,107 drillers' logs were interpreted as a sequence of sand facies and clay facies. Logs with ambiguous descriptions, depths and locations were disregarded. This study used lithologic descriptions of drillers' logs down to 152 m (500 ft). Since some of electrical logs and drillers' logs may be available at the same location, we prioritized geological information of electrical logs over drillers' logs. Moreover, screen intervals of pumping wells were assumed at permeable formations and provided sand facies interval information. Figure 3 shows a composite well log data out of an electrical log, a drillers' log, and two well screens at the same location.

FIGURE 3

Figure 3. Integration of different well log data.

Lithofacies Modeling

This study adopted the lithofacies modeling technique in Vahdat-Aboueshagh and Tsai (2021) to construct a complex stratigraphy model shown in Figure 4. Lithofacies modeling involves domain translation and facies interpolation. In general, strata in Louisiana have changing dips (Vahdat-Aboueshagh and Tsai, 2021). However, the Chicot aquifer system appears to be homoclinal and dip south. Therefore, the lithofacies model was constructed with the same dip direction and dip angle.

FIGURE 4

Figure 4. A schematic of lithofacies modeling with domain translation and facies interpolation: (A) well logs in a dipping domain, (B) well logs in a non-dipping domain, (C) lithofacies modeling in the non-dipping domain, and (D) the stratigraphy model in the dipping domain.

Consider a cross section in the coordinate system (x,y,z) with a changing dip angle and a number of well logs that are interpreted and labeled as a sequence of sand facies and clay facies (Figure 4A). The well logs in the dipping cross section is translated to a non-dipping domain (Figure 4B) with respect to the pivot plane in a new coordinate system (x,y,w):

\begin{array}{l} f_{(x, y, z)} \to g_{(x, y, w)} & (1) \end{array}

where f_{(x, y, z)} is the facies label data from well logs in the dipping domain, and g_{(x, y, w)} is the facies label data of well logs in the non-dipping domain. The updip area with respect to the pivot plane is translated down and stretched from the projection top to the projection base. The downdip area with respect to the pivot plane is translated up and squeezed between the projection top and the projection base. An indicator interpolation method is performed to generate a non-dipping facies model (Figure 4C):

\begin{array}{l} ĝ_{(x, y, w)} = \sum_{i} λ_{i} g_{i (x, y, w)} & (2) \end{array}

where ĝ_{(x, y, w)} is the estimated facies label in the non-dipping domain, g_{i (x, y, w)} is the facies label data from well log i in the non-dipping domain, and λ_i in the weighting coefficient for an interpolation method. Finally, the non-dipping facies model is translated back to the dipping domain (Figure 4D).

\begin{array}{l} ĝ_{(x, y, w)} \to {\hat{f}}_{(x, y, z)} & (3) \end{array}

where ${\hat{f}}_{(x, y, z)}$ is the estimated facies label in the dipping domain.

An indicator natural neighbor interpolation method (Tsai and Li, 2009) was utilized to construct a stratigraphy model using the 29,918 well logs, the dip direction of 270° (southward), and the structure dip angle of 0.06° (Jones et al., 1956). The indictor value for sand facies is 1 and for clay facies is 0 (Pham and Tsai, 2017). A cutoff value was suggested based on the proportion of sand facies in well logs (Elshall et al., 2013). In this study, the cutoff value was 0.4 for the Chicot aquifer system. The stratigraphy model started at the land surface and extends to a depth of about 305 m (1,000 ft). The study area was discretized into 56 layers. Each layer had 40,755 cells with a resolution of 1 km with varied layer thickness. A MODFLOW-2005 grid (Harbaugh, 2005) was developed based on the stratigraphy model. Readers are referred to Pham and Tsai (2017) for the detailed MODFLOW-2005 grid generation techniques.

Chicot Groundwater Model Development

Different MODFLOW-2005 packages were set up in order to develop a MODFLOW model with monthly stress periods. A River package was used to include river reaches and temporal river stages. The location of USGS stream gages is shown in Figure 5. A Recharge package was created based on a fraction of monthly surficial recharge rates provided by the USGS from January 2004 to December 2014 (Reitz and Sanford, 2019a). The USGS recharge rates were estimated as the residuals between water supply and the rest of water budget components (Reitz and Sanford, 2019b). A Well package was developed to include agriculture, industry, and public supply wells. A Horizontal Flow Barrier package was developed to simulate the Tepetate fault system. The lakes were not included in the groundwater model because of very limit hydraulic connection to the Chicot aquifer system (Vahdat-Aboueshagh and Tsai, 2021).

FIGURE 5

Figure 5. Locations of USGS stream gages and groundwater observation wells. Measured and simulated groundwater levels in the wells with well name are shown in Supplementary Figure 1.

Historical groundwater data were prepared for the simulation period from January 1, 2004 to December 31, 2014. There were 211 USGS wells with 2,161 groundwater level data. The location of USGS observation wells is shown in Figure 5. The groundwater data of January 2004 were employed to determine the initial condition for the groundwater model by spatial interpolation. Groundwater observation wells close to the model boundaries were used for the time-variant specified-head boundaries. The values of groundwater level from observation were interpolated spatially and temporally to obtain groundwater levels at different layers of the model. The temporal river stage data of Sabine River and Atchafalaya River were also used to estimate the boundary values at the stream portions. A total of 132 monthly stress periods were used in the MODFLOW model.

USGS water wells and Borrok and Broussard (2016) indicate that salt water is limited to the southern part of the Chicot aquifer system. The majority of the study region is not impacted by denser groundwater flows. Therefore, the effect of density-dependent flows was assumed insignificant. A future study may include the density effect in the vicinity of the Gulf of Mexico.

Model Calibration With Parallel Computing

Model parameters to be estimated were hydraulic conductivity (K) of sand, specific storage (S_s) of sand, hydraulic conductivity of clay, specific storage of clay, river conductance for 4 rivers, a fraction of USGS surficial recharge rates, and fault hydraulic characteristic (fault permeability per unit width of a fault). This study considered heterogeneous K and S_s for sand. Other parameters were considered to be homogeneous for the purpose of simplicity. The study assumed that the USGS recharge data is only applicable to surficial sand facies that directly link to deep aquifers. Recharge to surficial clay is negligible due to very low K. We applied the pilot point approach (Doherty et al., 2010) to estimate K and S_s of sand at 57 locations of the USGS groundwater observation wells. The natural neighbor interpolation method (Watson, 1999) was used to derive spatial distribution of K and S_s of sand. In total, there were 124 model parameters to be estimated through optimization. The ranges of model parameter values in optimization are shown in Table 1.

TABLE 1

Table 1. Ranges and estimated parameter values for the Chicot groundwater model.

The Chicot groundwater model was calibrated using the parallel CMA-ES (Elshall et al., 2015) approach to minimize the root mean square error (RMSE) (Chai and Draxler, 2014) between the 2,161 observed and simulated groundwater levels:

\begin{array}{l} min_{p} RMSE = \sqrt{\frac{1}{L} \sum_{i = 1}^{L} {[h_{i} (p) - h_{i}^{o b s}]}^{2}} & (4) \end{array}

where h_i and $h_{i}^{o b s}$ represent the i^th simulated groundwater level and the i^th observed data, respectively, L is the number of observed data, and p is a set of model parameters.

The Nash–Sutcliffe model efficiency coefficient (NSE) was used to evaluate the model calibration result.

\begin{array}{l} NSE = 1 - \frac{\sum_{i = 1}^{L} {(h_{i} - h_{i}^{o b s})}^{2}}{\sum_{i = 1}^{L} {(h_{i}^{o b s} - {\bar{h}}^{o b s})}^{2}} & (5) \end{array}

where ${\bar{h}}^{o b s}$ is the mean of the observed groundwater data.

Parallel computation was performed on a supercomputer that has 380 nodes. Each node has 16 processors. A population size of 1,240, a ten-fold size of unknown parameters recommended by Pham and Tsai (2017), was used. The computation time for each model run was about 50 min on a 2.6 GHz Xeon 64-bit Processor. To maximize the computational efficiency, this study ran the parallel CMA-ES code on 1,240 processers. The data-model processes are presented in the flowchart, Figure 6.

FIGURE 6

Figure 6. A flowchart of groundwater model development with stratigraphy data, hydrogeologic data and hydrologic data.

Results and Discussions

Stratigraphy Construction

A stratigraphy model of 56 non-uniform model layers was constructed for the Chicot aquifer system. The layer thickness varied from 2.2 to 10.4 m. The distributions of surficial sand facies and surficial clay facies at land surface are shown in Figure 7A. Around 7% of the model domain is covered by the surficial sand facies. The surficial sand facies appear to the northern area, which is the essential recharge zone for the Chicot aquifer system (Jones et al., 1956). The surficial sand facies are also seen in the vicinity of the Sabine River in the northwest.

FIGURE 7

Figure 7. The Chicot stratigraphy model: (A) distributions of surficial sands and surficial clays, (B–D) north-south cross sections, and (E–G) east-west cross sections. The cross sections are magnified 50 times in vertical direction. The vertical datum is NAVD 88.

The north-south cross sections AA′, BB′, and CC′ (Figures 7B–D) show the surficial clay thickening toward the Gulf of Mexico. Most of the surficial sand facies in the north connect to the Chicot aquifer system, where precipitation and streams recharge the aquifer system. The Chicot aquifer system encompasses seven aquifers (Nyman et al., 1990): the Upper Chicot and the Lower Chicot aquifers in the east, the Undifferentiated sand in the central area, and the Shallow sand, the “200-foot” sand, the “500-foot” sand, and the “700-foot” sand in the west. The late Pleistocene-Holocene Shallow sand in the west (AA′ and EE′ cross sections) is underlain by the Pleistocene “200-foot” sand, the “500-foot” sand, and the “700-foot” sand (Nyman et al., 1990). The Shallow sand is sparsely interbedded in the top confining clay layer and has limited spatial extent, implying that the Shallow sand is not as productive as the other sand aquifers. The three sands grade from a fine-medium sand to coarse gravel (Nyman et al., 1990). They are separated by interbedded clays in a few areas, but are connected in the most areas.

The “200-foot” sand, the “500-foot” sand, and the “700-foot” sand (EE′ and FF′ cross sections) merge into the late Pleistocene Undifferentiated sand in the center of the aquifer system (Nyman et al., 1990). The Undifferentiated sand links to the Upper Chicot aquifer and Lower Chicot aquifer in the central and south of the aquifer domain. The Undifferentiated sand is much thicker than other sands in the aquifer system. The Pleistocene Upper Chicot and Lower Chicot aquifers are separated by clay beds as thick as 30 m (100 ft) in a few areas. The clay beds are not extensive (Williams and Duex, 1995). In general, sands in the north are relatively shallow and often interbedded with clays. Sands thicken and clays pinch out from the north to the south.

Model Calibration Result

The model calibration was completed in 300 h with the parallel CMEA-ES code and 1,240 processors. The estimated model parameter values are shown in Table 1. In general, the simulated groundwater levels show good agreement in terms of value, trend, or pattern with respect to the observed groundwater levels (Figure 8). The root mean square error is 2.14 m. The Nash–Sutcliffe model efficiency coefficient (NSE) of 0.98 indicates a good calibration result. However, it is understood that data uncertainty and model error do exist. The groundwater model could not perform satisfactorily at few observation wells. Comparisons of observed and simulated groundwater levels at selected wells (Figure 5) with long records are shown in Supplementary Figure 1. The model calibration result shows that 30% of the USGS surficial recharge amounts (Reitz and Sanford, 2019b) became groundwater recharge to the Chicot aquifer system. The mean monthly surficial recharge rates to the outcrop zone for each quarter in a year are shown in Supplementary Figure 2.

FIGURE 8

Figure 8. Scatter plot of observed groundwater levels vs. calculated groundwater levels in 2004-2014. The groundwater level datum is NAVD 88.

Temporal Water Budget Analysis

The water budget of the Chicot aquifer system was analyzed by investigating the storage change from January 2005 to December 2014 with respect to the initial storage in January 2005. The model data prior to 2005 was not analyzed in order to skip the influence from the initial condition. Figure 9A shows the cumulative storage change, the monthly storage change, and the monthly pumpage. The average monthly net gain in aquifer storage was estimated about three million m³. By the end of 2014, groundwater storage gain was estimated 350 million m³ with respect to the initial storage in January 2005. However, this storage gain was not significant, on average about 1-mm groundwater level increase across the entire domain per year. The storage increase is likely attributed to groundwater recharge increase and pumping decrease. Groundwater recharge to the Chicot aquifer in 2004–2015 (Reitz and Sanford, 2019b) was estimated about 2.4 times greater than that in 1951–1980 (Wolock, 2003). In addition, groundwater pumping decreased from year 2000 to year 2010 by 18% according to Sargent (2011). The model shows slight decrease in storage during the March-July period due to pumping activities in the rice growing season. The aquifer system is quickly replenished during the non-irrigation season owing to the high transmissivity in the pumping areas (Hartono, 2005). The analysis of annual net groundwater inflows and outflows (Figure 9B) shows that the effect of pumping is virtually counteracted by inflows from boundaries, surficial recharge and rivers.

FIGURE 9

Figure 9. Estimated groundwater budget for the Chicot aquifer system: (A) monthly pumpage (million m³), monthly storage increase/decrease (million m³), and cumulative groundwater storage change (million m³) with respect to January 2005, and (B) annual net groundwater inflow/outflow (million m³) from the model boundaries, the pumping wells, the rivers, and the surficial recharge along with annual groundwater storage change (million m³).

The concept of safe yield or sustainable yield, where groundwater withdrawals do not exceed surficial recharge rate, has been criticized (Alley et al., 2002) and may not be applicable to the Chicot aquifer system. Instead, the concept of “capture” may be proper to explain the sustainability of the Chicot aquifer, where all the contributing factors to the storage change should be taken into account (Bredehoeft, 2002; Zhou, 2009).

Groundwater Level Analysis

The contour map of groundwater levels in 4-m intervals at the end of December 2014 is shown in Figure 10A. The values were calculated by averaging groundwater levels vertically along different model layers. The groundwater level varies from 40 m above the vertical datum NAVD 88 in the north to 20 m below the NAVD 88 in the industrial area. There are two cones of depression. One is in the agricultural area in the east and the other is in the industrial area in the west. The groundwater level distribution pattern in Figure 10A aligns well with the USGS reports (Fendick and Nyman, 1987; Lovelace et al., 2004).

FIGURE 10

Figure 10. Groundwater level map of December 31, 2014 for the Chicot aquifer domain: (A) with pumping, and (B) without pumping. The groundwater level datum is NAVD 88.

In order to better understand how the pumping activity has altered the groundwater flow in the predevelopment condition, the groundwater model was run without pumping. Figure 10B shows the contour map of groundwater levels at the end of December 2014. The groundwater level pattern resembles the predevelopment groundwater levels according to Martin and Whiteman (1989). The differences between Figures 10A and 10B imply that the pumping activity has reversed the groundwater flow from a gulfward direction to a landward direction. As Jasechko et al. (2020) has shown, the saltwater intrusion in Southwest Louisiana can be connected to groundwater level decline in inland areas. However, the origin of the salt water is not the Gulf of Mexico as they argued. Indeed, since the sand beds of the Chicot aquifer system dip at a greater angle than the Continental Shelf slope, a direct connection between the aquifer system and Gulf of Mexico is not possible (Martin and Whiteman, 1989). The salt water in the Chicot aquifer system is of ancient origin when the sea level rose and resulted in landward deposition of marine environments (Weiss, 1990).

To further analyze the groundwater flow pattern in the Chicot aquifer system, the average groundwater levels were obtained in three zones shown in Figure 11. The zones were defined based on the stratigraphy model and groundwater level variations in different model layers. Zone 1 includes the “200-foot” sand, the “500-foot” sand, and the “700-foot” sand in the west. Zone 2 is for the Undifferentiated sand in the central area. Zone 3 includes the Upper Chicot aquifer and Lower Chicot aquifer in the east.

FIGURE 11

Figure 11. Delineation of three zones, where Zone 1 includes the “200-foot” sand, the “500-foot” sand, and the “700-foot” sand, Zone 2 includes the Undifferentiated sand, and Zone 3 includes the Upper Chicot aquifer and Lower Chicot aquifer. Blue dots are pumping wells.

Monthly average groundwater levels in 2005-2014 in different sand systems are shown in Figure 12A. The average groundwater level in the Undifferentiated sand is always less than that in the Upper Chicot and Lower Chicot aquifers, indicating that groundwater generally flows from the Upper Chicot and Lower Chicot aquifers into the Undifferentiated sand. The heavier pumping in the Upper Chicot sand made the average groundwater level in the Upper Chicot aquifer lower than that in the Lower Chicot aquifer. The seasonal groundwater level decline and recovery in the Upper Chicot aquifer, the Lower Chicot aquifer, and the Undifferentiated sand correspond to the seasonal agriculture pumping. Figure 12A shows fast groundwater level recovery during the non-irrigation season. The low groundwater level variation in the “200-foot” sand, the “500-foot” sand, aind the “700-foot” sand shows that groundwater levels in the west are mildly affected by the agricultural pumping. Figure 12A shows downward gradients from the “200-foot” sand to the 500-foot sand and slight upward gradients from the “700-foot” sand to the “500-foot” sand. This verifies the fact that the “500-foot” sand is the most pumped sand in the west.

FIGURE 12

Figure 12. (A) Mean monthly groundwater level (m) in different sands of the Chicot aquifer system from January 2005 to December 2014. The groundwater level datum is NAVD 88, and (B) estimated cumulative groundwater storage change (million m³) with respect to January 2005 for Zone 1, 2, and 3, and monthly pumpage (million m³).

Groundwater depletion in the Chicot aquifer system (Konikow, 2013) mainly occurs in the agricultural zones as shown in Figure 12B, which reflects the groundwater level decline in Figure 12A. By the end of 2014, groundwater storage depletion was estimated about 45 million m³ in the Undifferentiated sand and about 35 million m³ in the Upper Chicot aquifer and the Lower Chicot aquifer with respect to the beginning of 2005. Groundwater storage depletion by the industrial pumping is insignificant during 2005-2014.

A schematic view of flow directions and average flowrates between different sands in 2005 to 2014 is shown in Figure 13. The highest vertical groundwater flows occur between the Upper Chicot aquifer and the Lower Chicot aquifer, and between the “500-foot” sand and the “700-foot” sand. The dominant vertical flow directions indicate that potential salt water may vertically migrate from deep brackish formations to the Upper Chicot in the east and to the “500-foot” sand in the west. The highest horizontal groundwater flows are westward from Zone 3 to Zone 2. The overall horizontal flow is westward.

FIGURE 13

Figure 13. Flow directions and mean flow rates (m³/day) in 2005–2014 between different sands in the Chicot aquifer system.

The intra-annual monthly average groundwater levels in different sands are shown in Figure 14. The average groundwater levels drop in the Undifferentiated sand, Upper Chicot aquifer, and Lower Chicot aquifer during the rice irrigation season, which accords with the fact that most of the wells in these sands are for rice irrigation. However, the relatively stable average monthly groundwater levels in the “200-foot” sand, the “500-foot” sand, and the “700-foot” sand were due to the year-round continuous groundwater pumping by the industries.

FIGURE 14

Figure 14. Intra-annual monthly groundwater level variation for different sands in the Chicot aquifer system from 2005 to 2014. The groundwater level datum is NAVD 88.

Conclusions

The methodology in this study is able to take into account structural dip information, build realistic groundwater models using big well log datasets and reveal aquifer complexity. The methodology was successfully applied to discovering Chicot aquifer complexity and developing a detailed Chicot groundwater model for the agriculture-intensive southwestern Louisiana, where the Chicot aquifer system is under high stress by irrigation pumping. More than 29,000 well logs in Southwest Louisiana were interpreted and prioritized. Prioritizing electrical logs, well screens, and drillers' logs in this order helps reduce errors in the lithofacies modeling as electrical logs have better quality than drillers' logs. The developed Chicot model is one of few models that reveals hydrogeological complexity in the Coastal Lowlands aquifer system.

The Chicot stratigraphy model reveals that the aquifer system consists of highly interconnected sand units and outcrops in the north. Central and southern areas are covered by clays thickening toward the Gulf of Mexico. Sand units are relatively thin in the north and are frequently interbedded by clays. Sand units are thick in the central and southern areas. The Upper Chicot aquifer and the Lower Chicot aquifer are interconnected in the east. The “200-foot” sand, the “500-foot” sand, and “700-foot” sand are interconnected in the west. The Undifferentiated sand is thick in the central area. These aquifers are laterally connected.

The groundwater modeling result shows that the storage loss due to groundwater pumping in the Chicot aquifer system is offset by inflows from surficial recharge, rivers and boundaries. The two large cones of depression created by the agricultural pumping in the east and by the industrial pumping in the west represent the key feature in the Chicot aquifer system. The storage loss occurs during the rice irrigation season. The aquifer system can be quickly replenished by the inflows during the non-irrigation season. The storage of the Chicot aquifer system was estimated to slightly increase during the period 2005–2014.

The groundwater model shows that groundwater levels in the Lower Chicot aquifer, the Upper Chicot aquifer, the Undifferentiated sand are heavily impacted by the seasonal irrigation activities and result in groundwater storge depletion. On the contrary, groundwater levels in the “200-foot” sand, the “500-foot” sand, and “700-foot” sand have less variability. In other words, groundwater levels in the west are not heavily impacted by the irrigation activities.

The groundwater model indicates large amounts of upward vertical flows from the Lower Chicot aquifer to the Upper Chicot aquifer and from the “700-foot” sand to the “500-foot” sand. As salt water was reported at the base of the Chicot aquifer system, potential saltwater intrusion due to vertical saltwater migration from the deeper sands is likely to occur in the future.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

Author Contributions

HV-A and FT contributed to conception and design of the study. DB and KP contributed to the pumping database. FT organized the database. HV-A performed the model development. HV-A and FT performed analysis. All authors contributed to manuscript preparation, read, and approved the submitted version.

Funding

This study was funded in part by the Louisiana Sea Grant College Program as part of Program Project ID R/CRM – 02 under NOAA Award NA18OAR4170098, the U.S. Geological Survey under Grant/Cooperative Agreement No. G16AP00056, and the U.S. National Science Foundation under Award No. 2019561.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The Louisiana State University High Performance Computing (HPC) is acknowledged for providing a supercomputer facility for this study. USGS is acknowledged for providing industrial and public supply pumping data.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/frwa.2021.623476/full#supplementary-material

References

Aarts, E. H., and van Laarhoven, P. J. (1987). “Simulated annealing: a pedestrian review of the theory and some applications,” in Pattern Recognition Theory and Applications. NATO ASI Series (Series F: Computer and Systems Sciences), Vol. 30, eds Devijver P.A., Kittler J, (Berlin: Springer). doi: 10.1007/978-3-642-83069-3_15