Study Area

The Mississippi Alluvial Plain region covers parts of seven States, extending as far north as Illinois and south to the Gulf of Mexico (fig. 1). The extent of the Mississippi Alluvial Plain covers about 44,200 square miles and can be divided into seven generalized subregions to allow for subregional comparisons of data (Ladd and Travers, 2019) (fig. 1A). The subregions are bounded by major streams or other groundwater divides, such as Crowleys Ridge. For this study, the southernmost subregion, the Deltaic and Chenier Plains, was expanded to include alluvial deposits along the Mississippi River between New Orleans and Lake Pontchartrain, La., and associated water-quality data (Killian and Knierim, 2023).

The alluvial aquifer is coextensive with the Mississippi Alluvial Plain region and overlies geologic units deposited by transgressive and regressive events during the Cretaceous and Tertiary Periods (Clark and Hart, 2009). The alluvial aquifer comprises Quaternary sands, silts, and clays deposited by the ancestral Mississippi and Ohio Rivers (Renken, 1998; Arthur, 2001). The thickness of the alluvial aquifer varies across the study area and is typically between 130 and 280 feet (ft) (Hart and others, 2008). In areas where the alluvial aquifer overlies subcropping aquifer units, the bottom of the alluvial aquifer is not well defined, and different sources provide different thickness estimates (Torak and Painter, 2019a; James and Minsley, 2021). A thin layer of shallow surficial confining material of silt, clay, and fine sand is present in the alluvial aquifer but may be absent in some areas (Renken, 1998). The surficial confining material is not uniform and is thought to impede recharge to the alluvial aquifer where present (Renken, 1998). The alluvial aquifer is above a large plunging syncline that dips toward the Gulf of Mexico. There are several subcropping aquifer hydrogeologic units that are Tertiary in age that have been mapped beneath the alluvial aquifer in the Mississippi embayment part of the study area (Hart and others, 2008).

Response Variable

Response variable data included spatially distributed measurements of SC collected from domestic wells and from public-supply, irrigation, and observation wells owned by several Federal, State, and local cooperators. SC water-quality data from 6,626 wells were obtained from the U.S. Geological Survey National Water Information System (U.S. Geological Survey, 2021). The Mississippi Department of Environmental Quality provided data from an additional 888 wells with SC water-quality measurements for a total of 7,514 wells (Killian and Knierim, 2023).

Wells that were less than 400 ft below the land surface were used in the response variable dataset, reducing the number of wells available for use in the models from 7,514 to 6,042 for the final dataset (Killian and Knierim, 2023) (fig. 3). Wells that were greater than 400 ft deep were not used because the depth of investigation of AEM survey data is generally 400 ft deep across the study area (James and Minsley, 2021). Wells completed in the alluvial aquifer and deeper Tertiary aquifers were included to allow the model freedom to predict without the constraint of aquifer designation because the base of the alluvial aquifer is not well defined in some areas. The median well depth was 114 ft below land surface (Killian and Knierim, 2023). A total of 987 wells did not have a defined bottom of well screen. Most of the wells were completed in the alluvial aquifer (number of samples [n] = 5,673), 98 wells had no aquifer designation, 49 wells were completed in terrace deposits, and the remaining wells were completed in deeper aquifers (n = 222).

Observed SC in the filtered dataset ranged from 4.7 to 23,400.0 microsiemens per centimeter at 25 degrees Celsius (µS/cm) (natural log 1.55 to 10.06) and was positively skewed (fig. 3) (Killian and Knierim, 2023). The median observation was 632 µS/cm (natural log 6.45). BRT models split the response variable data into trees that minimize variance, which will allow higher values to have a greater effect on variance. Because of how BRT models handle response variable data and because the data are positively skewed, the data were log-transformed using the natural log to reduce the relative importance of the higher values and increase the predictive accuracy for the bulk of the response variable data.

The caret package in R (Kuhn and others, 2019) was used to randomly split the response variable data into two groups: a training dataset and a holdout dataset. The training data consisted of 4,835 measurements (about 80 percent) and were used to tune the model; the remaining 1,207 measurements (about 20 percent) were held out to evaluate model performance.

Groundwater wells across the Mississippi Alluvial Plain region have varying numbers of samples ranging from one sample to continuous monitoring throughout a year or across several decades. If more than one observation existed for a well, the most recent measurement was used in model tuning and prediction. To increase the number of SC measurements for the model, samples from 1942 to 2020 were used (fig. 2). Censored values, or measurements where the value is only partially known (for example, below the detection limit of an instrument), were not included in the response variable dataset (n = 7). Some subregions of the Mississippi Alluvial Plain region—such as the Atchafalaya and the Deltaic and Chenier Plains (fig. 1A)—have not been sampled for several decades, whereas some subregions—such as the Cache and the Delta (fig. 1A)—are lacking historical measurements (fig. 2A). The number of SC measurements collected across the Mississippi Alluvial Plain region changed throughout time (fig. 2B). Because only the most recent SC measurement at each well was used in both models, information about changes in SC through time was effectively removed, and this modeling effort assumed no change in water quality over time at a well.

Explanatory Variables

Two BRT models were developed to understand the influence of geophysical explanatory variables on SC predictions: one model without geophysical variables (the “baseline” model) and one model with geophysical variables from an aquifer-wide AEM survey (the “geophysical” model). SC predictions were compared between the BRT models to quantify the change in model performance and the sensitivity of the models to geophysical variables (Knierim and others, 2020; Stackelberg and others, 2021).

Both models initially included 29 explanatory variables describing well construction information; surficial characteristics such as soils and geomorphology; positional information such as latitude, longitude, generalized subregions, and multi-order hydrologic position; depth to the water table; and year of SC measurement (table 1, at back of report). The geophysical model included the same 29 explanatory variables as the baseline model and an additional 19 geophysical variables (for a total of 48). Geophysical explanatory variables included bulk electrical resistivity (average electrical resistivity through the depth of investigation), which is sensitive to groundwater salinity and aquifer texture. Additional variables interpreted from the AEM survey included surficial confining unit top and bottom altitudes, the estimated base of the alluvial aquifer, and aquifer connectivity with streams (table 1, at back of report). More information about each explanatory variable, including name, description, dataset group, and source information, can be found in table 1 (at back of report) and Killian and Knierim (2023).

Explanatory variables were attributed to well points by using Python (Python Software Foundation, 2021) to execute a point extraction from the coincident raster cell, calculate averages for a buffer around the well, or extract a value from the depth coincident with the well screen. Point and buffer extractions were completed with the zonepy package (Clark and others, 2019). Point extraction was used to maintain high resolution of the source data, where applicable (table 1, at back of report), and for categorical explanatory variables. To account for spatial variation around a well, continuous explanatory variables were attributed using a 500-meter (m) buffer around the well. Well buffers also serve as surrogates for the contributing area (Johnson and Belitz, 2009). Two of the geophysical variables—aquifer resistivity and facies class—were attributed using the midpoint depth of the well screen and the corresponding coincident raster cell and depth slice from the AEM survey data. If well depth was unknown, resistivity and facies class were not assigned.

The year of SC measurement collection was included in the BRT models as an explanatory variable because the SC dataset spanned 78 years. Ultimately, time is inherently included in the model because of the decades over which SC measurements were collected. By explicitly including the year of measurement collection as an explanatory variable, the effect of time on SC predictions could be determined. Although the BRT models can make SC predictions for any specific year between 1942 and 2020, the models were not developed to predict SC through time. The explanatory variable dataset generally includes variables that capture conditions around 2018 (Killian and Knierim, 2023). To fully incorporate a time component in the BRT models and predict SC through time, explanatory variables that capture changes over time would also need to be included in model training.

Random noise was also included as an explanatory variable. The BRT models can determine the relative importance of each explanatory variable. Random noise was used as a benchmark to remove explanatory variables that were less influential than noise, which helped to improve model run time without impacting model performance. This method of including random noise for variable reduction is novel; however, similar approaches have been used (Ljunggren and Ishii, 2021). The random noise raster was generated with the ArcMap version 10.8 Create Random Raster tool available in the Spatial Analyst extension (Esri, 2020). The seed value was set to 2,640, and the Random Generator type was set to the default ACM599. Values in the random noise raster range from 0 to 1.

Geophysical Data Collection in the Mississippi Alluvial Plain Region and Incorporation in BRT Models

A unique aspect of the SC predictions from BRT models is incorporation of regional geophysical datasets, which provides information about the hydrogeology of the Mississippi Alluvial Plain region (Minsley and others, 2021). The geophysical data include an AEM survey, which measures resistivity signatures of water-bearing geologic materials (Smith and others, 2007). AEM survey data can help to map the thickness and extent of an aquifer. At the time of this study, the AEM survey data collected in the Mississippi Alluvial Plain region in 2019–21 covered 14,932 line-miles in the northern part of the study area and included shallow (less than about 330 ft) and deep (up to about 980 ft) estimates of inverted resistivity (mres) (James and Minsley, 2021; Minsley and others, 2021). Mres was further interpreted to estimate the bottom of the alluvial aquifer (mrvaBaseElv) and the presence or absence of surficial confining material (cnf_class). The AEM survey also included an CS-3 cesium vapor magnetometer (Scintrex Ltd., Concord, Ontario, Canada) to measure magnetic intensity (magint) and an RS-500 spectrometer radiochemical sensor (Radiation Solutions Inc., Mississauga, Ontario, Canada) to detect near-surface radioelements (predominantly potassium, uranium, and thorium) (Burton and others, 2021). These data provide three-dimensional information that has aided in the delineation of the base of the alluvial aquifer, locating where surficial confining material may be present or absent, and evaluating the degree of potential hydraulic connectivity between the alluvial aquifer and deeper subcropping aquifers (James and Minsley, 2021). At the time of this study, AEM coverage extended from the northern extent of the Mississippi Alluvial Plain region to approximately 31 degrees north latitude (figs. 3, 4), meaning that the AEM survey data did not cover the full Mississippi Alluvial Plain region. The authors sought to test how a BRT model would perform with an incomplete explanatory variable dataset.

Model Development

BRT models were developed in the R language (R Core Team, 2021) by leveraging the caret and gbm packages (Kuhn and others, 2019; Ridgeway, 2019). BRT is an ensemble-tree method that builds many trees with branches based on if/then statements that group data into similar homogenous groups with respect to the response variable (Kuhn and Johnson, 2013). The BRT method provides models that are more stable with better predictive performance than a single decision tree (Elith and others, 2008). The number and shape of trees in BRT models are controlled by five hyperparameters, which determine the model fit and complexity. Hyperparameters include the number of trees (nt), interaction depth (id; how deeply the trees split), shrinkage (sh; the proportion of each tree used as the model grows), minimum number of observations per node, and bagging fraction (bootstrap aggregation) (Elith and others, 2008; Kuhn and Johnson, 2013; Ridgeway, 2019).

Cross-validation (CV) tuning was used to tune the baseline and geophysical models and find the combination of hyperparameters with the best performance, or lowest root-mean-square error (RMSE) (Kuhn and others, 2019). Bagging fraction was held at the recommended value of 0.5 (Friedman, 2001) and minimum number of observations per node at 10. A tuning grid was constructed with the following ranges for the three other hyperparameters: nt (500–10,000 × 500), id (6–12 × 2), and sh (0.002–0.012 × 0.002). The tuning grid included 480 hyperparameter combinations (Killian and Knierim, 2023). Tenfold CV tuning was used, where 10 percent of the training data are used as a testing dataset to measure RMSE, and the process is repeated 10 times (Kuhn and Johnson, 2013). CV tuning was completed using the caret package (Kuhn and others, 2019) on a computer with 86 cores and was completed in approximately 20 minutes.

The BRT model with the lowest RMSE was considered the “best” model; however, the best model may be complex based on the hyperparameters, higher nt, id, and sh. Conversely, simpler models have lower nt, id, and sh. Models with high complexity may overfit to the training data and perform poorly when making new predictions. To find BRT models with similar predictive performance, models of lower complexity with RMSEs within one standard error (1SE) of the best model’s RMSE (Nolan and others, 2015) were inspected. Holdout data were used to evaluate 1SE model performance. Of the 480 hyperparameter combinations, there were 127 1SE baseline models and 172 1SE geophysical models. Because any 1SE model may provide reasonable predictions and could be selected as the “final” model, a “favorite” model was selected somewhat subjectively from the pool of 1SE models—one each for the baseline and geophysical models—favoring a simpler model (relatively low complexity controlled by the hyperparameters) and lower RMSE. Once the favorite model was selected, explanatory variables less important than the noise explanatory variable were removed on the basis of the assumption that random noise is insignificant to the model, and the model was retrained to obtain the “final” model. The final model was used to make SC predictions at wells and across the aquifer for both the baseline and geophysical models.

Ensemble-tree models, such as BRT, tend to underpredict high observations and overpredict low observations (Zhang and Lu, 2012; Belitz and Stackelberg, 2021). Empirical distribution matching (EDM) was used to correct for this bias and backtransform predictions from the natural log to original units (microsiemens per centimeter) (Belitz and Stackelberg, 2021). The EDM method calculates two empirical cumulative distribution functions for the training dataset: one for observed values and one for the ML predictions. Ordered pairs are matched between the two functions, which allow for linear interpolation between the points (Belitz and Stackelberg, 2021). As a result, if the training dataset includes the minimum and maximum observed values, then the range of EDM-corrected values will match the range of observations.

BRT models were inspected to better understand the importance and influence of explanatory variables for predicting SC. Variable importance was quantified using the gbm package in R (Ridgeway, 2019). Partial dependence plots (PDPs) quantify how a specific explanatory variable influences the response variable and were created using the pdp package in R (Greenwell, 2017). A PDP is generated from an average of individual condition expectation plots (Goldstein and others, 2015; Greenwell, 2017). These plots show the predicted value of the response variable (natural log of SC) on the y-axis if the explanatory variable of interest on the x-axis is changed over a range of observations and all other explanatory variables are kept at their observed or measured value (Elith and others, 2008; Greenwell, 2017; Knierim and others, 2020).

Prediction Mapping

SC predictions were made across the Mississippi Alluvial Plain region at multiple depths by using the final tuned BRT baseline and geophysical models and the “predict” function of the gbm R package (Ridgeway, 2019; R Core Team, 2021). Explanatory variable rasters were resampled to the 1-kilometer National Hydrogeologic Grid (Clark and others, 2018), and a flat file of the values (“raster stack”) was created. Predictions were made on the National Hydrogeologic Grid at 20-, 75-, 100-, 200-, and 300-ft depths below land surface. Prediction depths coincide with approximate depths to water at cones of depression and the depth of investigation for the geophysical data (Minsley and others, 2021; McGuire and others, 2023). Continuous prediction surfaces were generated for the entire Mississippi Alluvial Plain extent, and a water-level mask was applied to remove predictions that were made above the 2018 alluvial aquifer potentiometric surface (McGuire and others, 2020a). The year 2018 was chosen for the prediction year because it was coincident with the geophysical data collection (James and Minsley, 2021; Minsley and others, 2021) and the spring 2018 potentiometric surface (McGuire and others, 2020b) for the alluvial aquifer (spring is when water levels are assumed to be at maximum recovery, before the start of the main irrigation season). Predictions were bias corrected with the EDM method and backtransformed to units of microsiemens per centimeter.

Uncertainty

The 1SE models were used to quantify the range of SC predictions at the 200-ft prediction depth. Any of the 1SE models may provide a reasonable prediction; leveraging all of the 1SE model predictions provides a way to quantify the variability in SC predictions across the Mississippi Alluvial Plain region. This variability is used to express model uncertainty, though confidence intervals were not calculated. Prediction maps were generated for the 127 1SE baseline models and 172 1SE geophysical models. The rasters were then used to calculate the 5th percentile, 50th percentile (median), 95th percentile, and difference between the 95th and 5th percentiles of SC predictions for both models at each cell.

SC predictions at the 5th percentile, median, and 95th percentile were used to provide qualitative classifications of uncertainty relative to a threshold. TDS are used to define brackish and saline groundwaters, and SC and TDS show a highly correlated linear relation across the Mississippi Alluvial Plain region (Knierim and others, 2020). Therefore, TDS were calculated from SC predictions by using the correlation workflow in Knierim and others (2020). The TDS prediction rasters were then used to define where the models would very likely or likely predict TDS greater than (>) or less than (<) 1,000 mg/L, the brackish water threshold. The following logic was used:

If TDS-5th < 1,000 mg/L and TDS-median < 1,000 mg/L and TDS-95th < 1,000 mg/L, then TDS very likely < 1,000 mg/L.

If TDS-5th < 1,000 mg/L and TDS-median < 1,000 mg/L and TDS-95th > 1,000 mg/L, then TDS likely < 1,000 mg/L.

If TDS-5th < 1,000 mg/L and TDS-median > 1,000 mg/L and TDS-95th > 1,000 mg/L, then TDS likely > 1,000 mg/L.

If TDS-5th > 1,000 mg/L and TDS-median > 1,000 mg/L and TDS-95th > 1,000 mg/L, then TDS very likely > 1,000 mg/L.

BRT Model Results and Performance

The “final” BRT models underpredicted high observations and overpredicted low observations, which is typical of decision tree ML models (Zhang and Lu, 2012) (fig. 5A, B). When predictions were bias corrected with the EDM method, the predicted SC values matched the range of observed values because the training dataset included the minimum and maximum observations (fig. 5C, D). For the baseline model, EDM bias-corrected SC predicted at each well ranged from 4.7 to 23,400.0 µS/cm (natural log 1.55 to 10.06), with a median prediction of 634.0 µS/cm (natural log 6.45) (Killian and Knierim, 2023). For the geophysical model, the EDM bias-corrected SC predicted at each well ranged from 4.7 to 23,400.0 µS/cm (natural log 1.55 to 10.06), with a median prediction of 631.5 µS/cm (natural log 6.45) (Killian and Knierim, 2023). The EDM-bias correction method (Belitz and Stackelberg, 2021) provided better matches between the observed and predicted values for both models.

Overall BRT model performance, reported as coefficient of determination (R²) and RMSE, was similar for the baseline and geophysical models (fig. 5, table 2). The EDM bias-corrected baseline model training R² and holdout R² values were 0.896 and 0.617, and the EDM bias-corrected geophysical model training R² and holdout R² were 0.906 and 0.626 (fig. 5) (Killian and Knierim, 2023). There was a slight increase in model accuracy based on R² value of holdout data for the geophysical model (fig. 5). Training data predictions generally had higher accuracy (lower RMSE) than did holdout data predictions, which is typical for BRT models. Training data accuracy is typically higher than holdout data accuracy because training data were used to develop the model and holdout data were used to assess how well the model can predict using new data. Holdout RMSE values ranged from 0.42 to 0.45 (table 2) (Killian and Knierim, 2023).

Each Mississippi Alluvial Plain subregion had a different number of SC measurements, and the range and accuracy of SC predictions varied by subregion (fig. 6). Training R² values were similar among subregions, ranging from 0.79 to 0.92 (Killian and Knierim, 2023). Holdout R² values ranged from 0.26 to 0.88 (Killian and Knierim, 2023). The geophysical model had a higher R² value than the baseline model for the holdout data for all subregions except the Atchafalaya and Grand Prairie. Holdout R² values for both models were closest in the Delta subregion, where the models had the highest number of SC measurements. Holdout R² values were highest for both models in the Boeuf subregion. The Saint Francis subregion had the lowest holdout R² value for the baseline model, where SC is generally lower relative to the rest of the dataset. The Atchafalaya subregion had the lowest holdout R² value for the geophysical model. Holdout R² values differed the most in the Deltaic and Chenier Plains subregion, where SC is highest and the number of observations is lowest.

When the models are compared by depth slice, SC predictions tend to be higher and more accurate for the geophysical model compared to the baseline model. Holdout R² values generally increased with depth for both models (fig. 7). Both models also had holdout R² values less than 0.2 for wells in the 20-ft depth slice (0–20 ft below land surface). Wells less than 20 ft deep (that is, closest to the water table) were the least well defined (n = 78 wells). The geophysical holdout had a higher R² value compared to the baseline holdout for all depths except for wells in the 75-ft depth slice.

The geophysical model tended to predict higher maximum SC values at each depth slice compared to the baseline model, ranging from 30 to 200 µS/cm higher, although the median predictions were nearly identical (fig. 8A). Once the models were bias corrected using the EDM method, the difference in maximum predictions among depth slices between the baseline and geophysical models was less noticeable (fig. 8B). Therefore, the geophysical explanatory variables may provide additional information to the BRT model to predict high SC values prior to bias correction.

SC Raster Predictions and Uncertainty

Both the baseline and geophysical models predicted a general increase in SC from north to south, with the highest SC around the Gulf of Mexico (pl. 1). Additionally, both models predicted high SC in areas of known high salinity: around Iberville Parish, La. (Welch and Hanor, 2011); in southeastern Arkansas (Larsen and others, 2021); and around White County, Ark. (Bryant and others, 1985) (figs. 1, 3). The baseline model predicted higher SC near the Gulf of Mexico compared to the geophysical model; frequency of predictions increased between 5,000.1 and 28,000.0 µS/cm in the baseline model (pl. 1, A). The geophysical explanatory variables do not extend to the Gulf of Mexico (figs. 3, 4), so a true one-to-one evaluation of model predictions cannot be quantified. However, a one-to-one comparison of models with identical explanatory variable extents coincident with AEM survey data coverage was explored, and no significant impacts to model performance were discovered because of the nature of how BRT models function. The comparison between models with differing explanatory variable extents was sufficient to assess the impact of AEM survey data on the model, as well as to assess the ability of a BRT model to predict where explanatory variable data are lacking. The extent of AEM survey data coverage created a horizontal linear feature in the geophysical model for all depth predictions (pl. 1, B). The difference in predictions between the baseline and geophysical models increased at the 200- and 300-ft depth slices compared to shallower predictions (pl. 1, C).

The geophysical model predicted higher SC with increasing depth compared to the baseline model, especially above 31 degrees north latitude, where the geophysical datasets were collected (figs. 8, 9, pl. 1). For the 20-, 75-, and 100-ft depth slices, the baseline model predicted a larger interquartile range and higher median than did the geophysical model (fig. 9). The interquartile ranges for the 200- and 300-ft depth slices for both models were similar, but the geophysical model had higher medians and ranges of predictions (fig. 9). This difference between the baseline and geophysical model predictions is more noticeable in areas of known high SC where geophysical explanatory variables are complete (that is, no missing data), such as southeastern Arkansas and northeastern Louisiana (pl. 1).

Uncertainty in SC predictions was high where observed SC was high and where there was a greater range in SC predictions among all 1SE models in both the baseline and geophysical models (fig. 10A, B). High relative uncertainty where the response variable is high has been noted in other ML models (Knierim and others, 2020). The importance of quantifying uncertainty in model predictions depends on how the predictions will be used; for example, whether the modeling effort is focused on accuracy for low or high observed values. High groundwater salinity can limit drinking water and irrigation uses; therefore, qualitative uncertainty was mapped to show where the model is very likely or likely to predict greater or less than the brackish groundwater threshold of 1,000 ppm TDS (fig. 10C, D) (U.S. Geological Survey, 2013). Overall, there were few areas where either model does not predict fresh versus brackish groundwater with a high level of confidence. Confidence was determined by whether the 5th percentile (considered the lowest prediction), median, and 95th percentile (considered the highest prediction) all predicted SC greater or less than 1,000 ppm TDS.

Explanatory Variable Influence

Both the baseline and geophysical models had fewer explanatory variables in the final model compared to the potential explanatory variables used in initial model tuning (table 2). Twenty-two of the 29 explanatory variables were more important than noise and retained in the final baseline model (fig. 11A). The geophysical data added 19 potential explanatory variables. Thirty-three of the 48 explanatory variables, including 12 explanatory variables from the AEM survey, were more important than noise and retained in the final geophysical model (fig. 11B). The same six explanatory variables were found to be less influential than noise in both models, including variables that describe soil texture and the age of surficial sediments (fig. 11).

PDPs for continuous explanatory variables shared between both models show a similar pattern for predicting SC (fig. 12). For example, SC was predicted to be lowest at approximately −92 degrees longitude and highest at approximately −91 degrees longitude in both models (fig. 12E), despite longitude being the first and third influential for the baseline and geophysical models, respectively (fig. 11). The baseline model PDPs generally predicted a higher range of SC than did the geophysical model PDPs for the shared explanatory variables (fig. 12). The depth of the well screen (BOT) was one explanatory variable that showed a different relation between the two models; SC increased with depth in the baseline model and decreased with depth in the geophysical model (fig. 12B). Eleven continuous explanatory variables (fig. 12U–AE) were unique to the geophysical model, and the variable with the largest range of SC predictions was mres (fig. 12Z).

Missing Data Versus Structured Nulls in BRT Models

BRT models can make predictions even where explanatory variables are missing data because the decision trees are built with the missing data, or null values. Both the baseline and geophysical models were able to predict SC in the southern part of the Mississippi Alluvial Plain region, where there were fewer SC observations (fig. 3, pl. 1). In addition to the paucity of SC observations, geophysical data were not collected below approximately 31 degrees north latitude in the Mississippi Alluvial Plain region prior to or during model development (figs. 3, 4); therefore, there is a large area of nonrandomly distributed (or “structured”) nulls. Because overall model performance was not impacted by the presence of a structured null, as assessed by a one-to-one model comparison with identical extents not described herein, the ability of a BRT model to predict with a structured null as a result of differing explanatory variable extents was tested. For wells in this area, the geophysical model had fewer explanatory variables on which to train. The large structured null created a linear feature in SC predictions for the geophysical model at approximately 31 degrees north latitude (pl. 1, B). Although BRT models can make predictions where there are randomly distributed missing data in explanatory variables, the BRT models do not appear to compensate for structured nulls.

SC was predicted to be lower in the geophysical model relative to the baseline model in the southern extent of the study area (pl. 1), despite the general trend of increasing observed SC toward the Gulf of Mexico (fig. 3). SC predictions north of the linear horizontal feature created by the structured null for the geophysical model were similar between the two models up to approximately the 200-ft depth slice (pl. 1). It cannot be determined which model would predict higher SC in the southern part of the Mississippi Alluvial Plain region if geophysical data were available for the full extent. Although predictive performance is similar between the two models in the Deltaic and Chenier Plains subregion (fig. 6), SC predictions from the geophysical model are likely not reasonable because of the large structured null coincident with missing geophysical explanatory variables.

Position in the Groundwater System Influences SC

The most influential explanatory variables for both the baseline and geophysical models indicate that position in the groundwater system is an important predictor of SC across the Mississippi Alluvial Plain region (fig. 11). A similar result was observed in a previous study in the Mississippi Alluvial Plain region where surrogates for position along groundwater flow paths were important predictors of SC and chloride (Knierim and others, 2020).

The top three most influential explanatory variables for the baseline model (DEC_LONG_V, GW2018_ft, and region [subregion]) are indicators for position in the groundwater system (fig. 11A). SC was predicted to be highest at approximately −91 degrees longitude (fig. 13A), which is the approximate center of the Mississippi Alluvial Plain and where the alluvial aquifer tends to be thickest (Clark and others, 2011). The 2018 alluvial aquifer potentiometric surface (fig. 13C) also reflects position in the system because groundwater altitude is higher in the northern part of the Mississippi Alluvial Plain region, where land surface altitudes are higher compared to sea level. SC was predicted to be high where the potentiometric surface is low (fig. 13C), which occurs toward the Gulf of Mexico. A similar north-to-south trend of increasing SC toward the Gulf of Mexico is evident in SC observations and mapped predictions (fig. 3, pl. 1). SC was also influenced by the subregion where the SC measurement was collected (figs. 11A, 13). The southernmost subregion, the Deltaic and Chenier Plains, tends to have higher SC because of saltwater influence from the Gulf of Mexico (Welch and Hanor, 2011).

The geophysical model was also influenced by explanatory variables that describe position in the groundwater system. For the geophysical model, the most influential explanatory variable was the bottom altitude of the surficial confining material (cnf_bot) (fig. 11B). Like groundwater altitude, the altitude of the bottom of surficial confining material tends to mimic land surface altitude and decreases from north to south. SC was predicted to be higher where the cnf_bot altitude was low, between 0 and 20 m (fig. 13B), which occurs in the southern half of the Boeuf and Delta subregions. The second most influential explanatory variable for the geophysical model was the bottom of the well screen relative to land surface (BOT) (fig. 13D). Generally, SC would be expected to increase along groundwater flow paths because dissolved solids increase with greater rock-water interaction; depth can serve as a surrogate for flow path (Moore and others, 2019). However, the geophysical model predicted lower SC with increasing well bottom screen depth, which is opposite of the baseline model (fig. 12B) (see below for discussion).

Variation in Model Performance and SC Predictions With Depth

The baseline and geophysical BRT models had similar overall predictive performance despite being tuned with different explanatory variables (tables 1, 2, figs. 5, 7, 11). The predictive performance of both models was also similar to a previous study that predicted SC by using a BRT model in a smaller part of the study area (Knierim and others, 2020). The poor model performance at the 20-ft depth slice (n = 78) may be explained by the lack of SC observations near the water table and thus a lack of training data (fig. 7). Because of the poor model performance at this depth slice, the 20-ft depth prediction is likely inaccurate. More SC observations near the water table would be required to increase prediction accuracy at shallow depths; with the current training data, the BRT models should not be used to make predictions less than 75 ft deep.

Overall, the baseline model predicted similar SC at wells for all depth slices and across all raster depth slices, whereas the geophysical model predicted increasing median SC with depth at wells and raster depth slices (figs. 8, 9, pl. 1). Observed SC was not different among most depth slices, except for higher SC at the 300-ft depth slice (fig. 8), which likely includes wells below the alluvial aquifer. The baseline model included explanatory variables that predominantly describe surficial conditions across the Mississippi Alluvial Plain region (such as multi-order hydrologic position or soils) or two-dimensional position in the system (fig. 11A). BOT was the only three-dimensional explanatory variable in the baseline model and was an important predictor in both models (fig. 11). The geophysical model incorporated additional three-dimensional data including mres measurements that vary with depth within the system (figs. 4, 11B). With the addition of geophysical data, surficial variables decreased in importance in the geophysical model, but position in the groundwater system was still the predominant driving predictor of SC in the Mississippi Alluvial Plain region (fig. 11B). The addition of three-dimensional data in the geophysical model may explain why the median SC predictions at the 200- and 300-ft depth slices are higher than in the baseline model (fig. 9, pl. 1).

The PDPs for BOT show contrasting behavior between models in comparison to mapped predictions; SC was predicted to increase with well screen bottom depth in the baseline model, but it decreases in the geophysical model (fig. 12B). This contrasting relation may be related to the structured nulls in the geophysical model, which led to generally lower SC predictions in the geophysical model compared to the baseline model south of the geophysical data extent (pl. 1). Each PDP provides SC predictions if all other explanatory variables are held at their measured or observed values, except for BOT. The mapped predictions integrate spatial changes in all explanatory variables. Therefore, the BOT explanatory variable may modify the influence of other explanatory variables in predicting SC differently between the models. If AEM geophysical data were available for the entire Mississippi Alluvial Plain extent, the BOT PDP for the geophysical model may display a trend similar to that for the baseline model.

Geophysical Model Explanatory Variables and Sources of Salinity

Areas of elevated SC in the Mississippi Alluvial Plain region have been suggested to be controlled by several possible sources of salinity, including ocean water; salt domes; connection with deeper, more saline subcropping aquifers; or surficial concentration through evaporation in areas of low infiltration controlled by the presence or absence of surficial confining units (Pettijohn and others, 1988; Williamson and others, 1990; Kresse and Clark, 2008; Kingsbury and others, 2014). The north-to-south trend of increasing SC observations and predictions illustrates that ocean water influences groundwater SC in the southern part of the study area close to the Gulf of Mexico (pl. 1, callouts 3 and 6). Salt domes were not included as an explanatory variable and so were not tested as a source of salinity in the BRT models. The addition of geophysical explanatory variables added three-dimensional information about hydrogeology, including the location and thickness of confining units at the surface and with depth (James and Minsley, 2021; Minsley and others, 2021), that was unavailable to the baseline model. These additional variables allowed exploration of how surficial confinement or connection with deeper units may influence SC.

Relative Influences of Lithology and Salinity

Mres at the midpoint of the well screen was the seventh most influential explanatory variable in the geophysical model (fig. 11B). Mres captures both aquifer properties and porewater salinity, and previous analysis of the Mississippi Alluvial Plain region has shown that saline groundwater can influence resistivity measurements (Minsley and others, 2021). There is an inverse relation between SC and resistivity (fig. 14). Therefore, the geophysical model provides additional evidence that mres includes a resistivity signal driven, in part, by high SC (or saline) groundwater (fig. 14). Mres provides important three-dimensional information to the geophysical model and can be used to delineate contacts between the alluvial aquifer and deeper Tertiary confining units or aquifers (figs. 1, 4). However, because of the presence of saline groundwater in some areas of the alluvial aquifer (fig. 3), low mres (or high conductance) signals in the AEM survey data are not well constrained (Ward, 1990; Hilmi and Madun, 2022). Low resistivity can be the result of high salinity, the presence of fine-grained lithology, or a combination of the two (Zohdy and others, 1974). More work is needed to resolve the groundwater SC (salinity) signal within AEM-derived resistivity data across the Mississippi Alluvial Plain region. Removing the salinity signal could allow for better estimation of subsurface materials, thus allowing for more accurate estimates of available water resources.

Lack of Temporal SC Variations at Wells

The year an SC measurement was collected was included as an explanatory variable in both models and had a similar influence on SC predictions (fig. 11). The SC dataset includes sampling networks from six States and includes long-term State cooperative programs to monitor salinity in Louisiana and Arkansas (Killian and Knierim, 2023), but network sampling priorities have changed over time, and there are gaps in the spatial and temporal coverage of data. For example, the Deltaic and Chenier Plains subregion shows elevated SC (fig. 3) but lacks data compared to other subregions (fig. 2A). Most SC measurements in this subregion were made prior to the year 1980 (fig. 2A). In contrast, there was an increase in SC observations after 1980 in the Delta subregion (fig. 2A). There may also be a bias introduced from networks that have the goal of monitoring localized high-salinity areas. The importance of this explanatory variable is likely driven by differences in sampling networks, as opposed to the BRT models capturing systematic changes in SC conditions throughout the aquifer.

In addition to the variation in the sampling network over time, groundwater quality may change at a location as water is withdrawn for irrigation or the aquifer is recharged (Kingsbury and others, 2014; Jurgens and others, 2020). Because only the most recent sample was used if a well had more than one observation, the model was not given information about how SC changed over time at a well. The SC maps were predicted using the year 2018 because 2018 was coincident with data collection for several explanatory variables. Although SC prediction maps could be made for any specified year, it would not be appropriate to use the maps to evaluate SC change through time across the aquifer. To understand how SC changes through time and which explanatory variables drive the change, explanatory variables that capture these changes through time are needed in addition to time-series SC observations.

Uncertainty in High SC in the Mississippi Alluvial Plain Region

The baseline and geophysical models showed overall similar trends in predicting broad areas of low SC in the northern part of the Mississippi Alluvial Plain region and high SC from southeastern Arkansas to the Gulf of Mexico (pl. 1). There are also several areas of localized known high SC in the Mississippi Alluvial Plain region (fig. 3), and both models were able to predict elevated SC in these areas (pl. 1). However, the spatial extent and magnitude of SC predictions were different between the baseline and geophysical models (pl. 1). Depending on which model is used, the same location can be predicted with high confidence to have either fresh or brackish groundwater. For example, the Grand Prairie subregion (fig. 1) is predicted to have more brackish groundwater in the geophysical model than in the baseline model (fig. 10). Both models have similar overall accuracy (fig. 5) and similar accuracy among prediction depth ranges (fig. 7) and subregions (fig. 6). Therefore, both models have a similar magnitude of overpredicting or underpredicting SC across the Mississippi Alluvial Plain region (fig. 5).

The appropriate model and prediction maps to use depend on the goal and area for assessing SC conditions. Both models showed poor prediction accuracy for shallow SC conditions (fig. 7). SC predictions near the Gulf of Mexico from the geophysical model may not be realistic because of the structured nulls in training data (fig. 6, pl. 1). Additionally, the high SC predictions in northeastern Arkansas from the geophysical model may be spurious because of the lack of high SC observations in that area. In contrast, the baseline model was trained with no three-dimensional data other than well depth, so the geophysical model may be more appropriate to use when making deeper SC predictions (fig. 7). Both the baseline and geophysical models can make accurate SC predictions for the Mississippi Alluvial Plain region with the given input data, but model outputs should be selected and interpreted on the basis of the scientific goal.