Spatiotemporal Analysis of Available Freshwater Resources in Watersheds Across Northern New Jersey

1. Introduction

Water resource is critical to virtually every sector of human activities; from domestic water supply, agriculture, industrial use to other varied public water supply purposes [1,2,3,4]. Groundwater is a critical source of this water supply, as the United States utilized about 281 billion gallons of freshwater per day in 2015. About 30% of the total freshwater withdrawn was from groundwater sources. The same year, New Jersey withdrew about half a billion gallon of freshwater per day from groundwater supplies [1]. Since groundwater represents a vital component of the United States’ water supply, it is critical to ensure the long-term sustainability and quality of this valuable resource. Due to the huge demand, several studies across the United States are already observing the deteriorating quality of groundwater resources [4,5,6,7].

There are several controls driving the changes of groundwater quality. Lindsey et al. [8] studies of water chemical composition across the United States shows that increasing levels of Na and Cl ions serve as the main indicator for changes in groundwater chemistry. This was particularly observed in urbanized regions. In the United States' temperate regions, groundwater contamination poses a significant threat to water resources from varied sources. One major contamination source is urban runoff. During precipitation events, runoff carries pollution-ladened effluents into important water resources. Deicing salts applied on road surfaces during winter months infiltrate into groundwater, thereby compromising the water quality. Rural septic tank leaks are another source of contamination. Sewage waste from faulty septic systems contaminates groundwater, mostly impacting areas without centralized water and sewer infrastructure. The sewer system overflows carry every form of contaminant during heavy flooding events. Heavy storms and flooding cause combined sewer systems (such as found in New Jersey) to overflow, mixing with pollution-ladened runoff and seeping into groundwater. Aging sewer systems and inadequate infrastructure capacity often lead to overflows, leaks, and contamination. Generally, these combined sewer systems carry a mix of stormwater and sewage overflow during intense storms through manhole and line leaks, resulting in groundwater pollution [9].

Many studies in the United States have linked increased anthropogenic activities such as the use of deicing salt and rapid urbanization to continual deterioration of the groundwater quality [10,11,12]. Increasing groundwater contamination such as salinization has been observed to enhance the transport of contaminants from the land surfaces into the groundwater. Constituents like sodium ions from deicing salts application have been observed to release heavy metals from soil substrates into surface waters and groundwater [13,14]. The United States consumed about 58 million metric tons of salt in 2019, with about 43% of this salt (which are mainly in the form of sodium chloride) reported to be applied to Highway cleaning [11,15]. Roads in the urban regions of northern New Jersey consume about a third of a million metric tons of deicing salt annually. Reports from the New Jersey Department of Transportation (NJ DOT) reveal that New Jersey used about 10,000 kg of NaCl for every lane kilometer (0.621 mile) in 2017 [16]. The northern New Jersey's region covers about 18,000 miles of urban public roads [17]. This is the key reason for its much-needed application of deicing salt during wintry conditions.

Specific conductance serves as a valuable proxy for evaluating groundwater quality, especially in areas impacted by runoff from urbanization and agricultural activities, where increased levels of conductance indicate eutrophication and pollution [18]. Specific conductance indicates a broad insight into the levels of water quality, representing the various contaminant levels available in the water. Specific conductance values are largely used for extrapolating salinity levels in water [19,20,21]. Cartwright et al. [22] indicates that Specific conductance can be applied as a proxy technique for establishing salinity. The metric unit for specific conductance of water is in microsiemens per centimeter (μS/cm) and this provides a measure of salinity for various water sources; ocean water (55,000 μS/cm), distilled water (0.5 – 3.0 μS/cm), melted snow (2 – 24 μS/cm), portable water (30 – 1500 μS/cm), tap water (50 – 800 μS/cm) [23]. A level of 150 – 500 μS/cm is usually desired for general purposes [24,25]. United States Salinity Laboratory [26] classifies salinity with respect to specific conductance, as low salinity (100 – 250 μS/cm), medium salinity (250 – 750 μS/cm), high salinity (750 – 2250 μS/cm), and very high salinity (> 2250 μS/cm).

Increases in salinity have been observed in northern New Jersey, as several anthropogenic sources (such as the use of deicing salt) have been ascertained to influence the quality of water resources [12,27]. In recent times, rapid urbanization that occurred in New Jersey has seen an increase of about 23% during the short period of 1995 – 2000, compared to observation from 1986 to 1995. These land use modifications were observed all across the state of New Jersey, with Morris, Sussex, and Somerset counties suffering the most intense loss of forested areas in northern New Jersey [28].

Continuous assessment of groundwater resources is essential to provide decision-makers and stakeholders with critical information on freshwater availability, ensuring it meets the growing demand for supply. However, to ensure the public can easily understand the status of their groundwater quality, there is a need to present this information in a simplified and accessible manner. Such clear and concise data is vital during community engagement meetings, enabling informed discussions and the development of sustainable solutions to address water quality challenges. Visualizing these data would provide a better understanding of the current status and trend of the groundwater quality. It is necessary to know whether a particular region is heading in the right direction in terms of work that needs to be done to ensure continual availability of freshwater.

The power of visualization in GIS can be adapted to provide a characterization of groundwater quality across the study area. The use of GIS would provide spatial variability of the water quality in the region across varied land uses. This study focuses on the use of geospatial analysis to estimate the specific conductance as a measure of the groundwater quality in watersheds across northern New Jersey. The study uses GIS to provide spatial and temporal variability in specific conductance in the various WMAs, considering land use and Highways in the region. The study provides a mapped-out correlation between land use, Highways and spatial variability of specific conductance across different WMAs in northern New Jersey. The total extent of areas with available freshwater across northern New Jersey are outlined using freshwater metrics as indicated by United States Salinity Laboratory Staff [26]. The study methodology provides an estimate of salinity as an indication of freshwater in areas without measured specific conductance and indicates the various levels of water quality in the WMAs in the region. The aim is to develop a geostatistical model to predict salinity levels as a measure of groundwater quality, using specific conductance values already measured from specific locations to determine the areal coverage of freshwater. The analysis would reveal the salinity concentration distribution over time and indicate changes occurring over that period. A simple data visualization of the effect of deicing salts on groundwater across the region is provided. The kriging interpolation method in GIS is used to develop the geostatistical model which provides a representation of salinity in the area. This technique quantifies salinity in the entire region using data from measured locations to estimate data in unmeasured locations. Several studies have applied varied spatial interpolation techniques in assessing water quality using various parameters as indicators of suitability or state of quality [29,30,31,32,33]. Research indicates that the Kriging interpolation technique outperforms other interpolation methods [34,35,36]. This technique is known for its ability to generate accurate estimates at unsampled points using values from neighboring measured points [29,37], while minimizing error [38]. It is also known to provide a more accurate distribution of parameters across geographical areas [39]. Currently, there is no spatial analysis study that provides a spatial distribution of water quality across northern New Jersey, nor any research identifying areas in the region with available freshwater. Therefore, this study provides decision makers and concerned individuals with valuable insights on the current status and availability of freshwater, as well as potential changes in groundwater quality overtime in northern New Jersey.

1.1. Study Area

Geologically, New Jersey is partitioned into two parts by a fall line, the north and south parts. The non – coastal plains region is the north part of the fall line while the coastal plains region is south of the fall line. This study is focused on New Jersey non - coastal plains water region which is also referred to as the northern New Jersey water region. Three water regions are located in northern New Jersey, the Northwest, Northeast, and Raritan water regions. There are eleven WMAs in the three regions; Upper Delaware (WMA - 1), Wallkill (WMA - 2), Pompton, Wanaque, Ramapo (WMA - 3), Lower Passaic and Saddle (WMA - 4), Hackensack and Pascack (WMA - 5), Upper Passaic, Whippany, and Rockaway (WMA - 6), Arthur Kill (WMA - 7), North and South Branch Raritan (WMA - 8), Lower Raritan, South River and Lawrence (WMA - 9), Milestone (WMA - 10), and Central Delaware [40,41] (Figure 1a). The study area is characterized by a unique bedrock geology. The physiographic provinces include, “the Valley and Ridge, Highlands, and Piedmont provinces.” [41]. Detailed literature on the bedrock geology of northern New Jersey is found in studies by Dalton et al. [42], Drake et al. [43] and NJDEP [44]. Northern New Jersey comprises primarily of old consolidated geological units covering about 3,116 square miles in area and makes up two-fifths of New Jersey State [41,45]. The area covers 12 counties of New Jersey; Sussex, Passaic, Bergen, Hudson, Essex, Morris, Warren, Hunterdon, Somerset, Union, and parts of Middlesex and Mercer [45]. The fall line cuts through WMAs 9, 10 and 11 (Lower Raritan, Millstone and Central Delaware WMAs), thereby placing some of these WMAs’ groundwater wells in the south part of the fall. In this study, all the wells within WMAs 9, 10 and 11 in both the north and south parts, were analyzed as part of the northern New Jersey region. This was done to harness as much data as possible.

The northern New Jersey region comprises agricultural, undeveloped, and urban land use (Figure 1b). The Northwest water region has the largest agricultural area. The Northeast water region has no agricultural land use area, and it is the most urbanized water region. The Land use areal markings of the data used in this study were obtained from “1986 and 1995 land-use geospatial mapping” [46].

2. Methods

2.1. Data Collection

To spatially assess the groundwater quality of the watersheds of northern New Jersey, available data were analyzed and mapped out. The data was obtained from the achieves of the New Jersey Ambient Groundwater Quality Monitoring Network (AGWQMN) database, containing hydrogeological data from shallow wells across New Jersey. The data furnish information on “nonpoint-source contaminants” impacting the groundwater. The New Jersey AGWQMN data used in this study cover four sampling time cycles (cycle 1: 1999 – 2004, cycle 2: 2004 – 2008, cycle 3: 2009 – 2013 and cycle 4: 2014 – 2016). Two hundred and ninety-five (295) well observations from the water regions in northern New Jersey were obtained for this study. Groundwater wells used in this study were completed in unconfined aquifers to enhance insights on groundwater contamination from the near surface [41]. This study was evaluated based on four sampling cycles. Cycles 1, 2, and 3 included data from 76, 77, and 77 wells, respectively, covering all three land use types: agricultural, undeveloped, and urban. Cycle 4 was split into Cycle 4 and 'Adjusted Cycle 4,' with data from 65 and 77 wells, respectively. Unlike the earlier cycles, Cycle 4 excluded data from 12 wells of undeveloped land use area. The 'Adjusted Cycle 4' combines well data from Cycle 4 with a set of "well data from the undeveloped land use" taken from Cycle 3. This adjustment assumes a scenario where the well data from the undeveloped land use during the Cycle 3 sampling period remained unchanged and were observed again in Cycle 4. As a result, Cycle 4 does not include undeveloped land use well data, while 'Adjusted Cycle 4' incorporates well data from the undeveloped land use area.

The collected data were cleaned and prepared for analysis using R, an analytical software [47], along with QGIS [48] and ArcGIS [49] for geospatial analysis in this study. The choice of the geochemical parameter ‘specific conductance’ for the spatial modeling and estimation of the groundwater quality was based on well documented literature [18,22].

The analysis process is computationally intensive, following a workflow (Figure 2) to meet the goal of the study. The process can be divided into four major tasks: (1) data collection and processing – exploratory data analysis, (2) develop/fit semivariogram model, (3) Cross – validation and (4) Kriging interpolation. The kriged predictions were then mapped out and freshwater areas calculated.

2.2. Data Analysis

This study uses kriging interpolation to provide a clear understanding of the level of salinity and the overall quality of the groundwater in Northern New Jersey’s 11 WMAs, over a time period. Specific conductance values, which are a measure of salinity, were used to spatially evaluate the groundwater quality in the WMAs across the region, indicating varied water quality distribution over the several water sampling cycles.

2.2.1. Kriging Interpolation

Kriging is a geostatistical interpolation technique utilized for predicting values at unsampled locations based on measurements from known locations. Kriging interpolation is distinguished by its ability to assess estimation errors and to quantify uncertainty in the interpolated surfaces. Studies have indicated that Kriging provides more accurate predictions compared to other deterministic and stochastic methods, making it a preferred technique for spatial analysis [35]. Kriging interpolation is one of the most efficient and widely applied techniques in geospatial analysis of environmental variables [2,31,32,33,50]. The interpolation methods include universal kriging (UK), simple kriging (SK), ordinary kriging (OK), and cokriging (CK) [51,52,53]. Among these, ordinary kriging is widely recognized, and it is applied in this study due to its simplicity, minimal assumptions, and superior prediction accuracy relative to other kriging variants [32,54]. Ordinary kriging can be used for the interpolation of spatial distribution with trend [53]. This variant is particularly suitable when the primary objective is to produce a prediction map. In addition to the prediction map, kriging interpolation also produces uncertainty, which makes it a more reliable interpolation techniques for prediction of unknown points [55].

2.2.2. Exploratory Data Analysis

Exploratory Data Analysis is aimed at analyzing the statistical and spatial attributes of the data. The data is explored to indicate the distribution patterns, data uniformity and trends. A plot of histogram is used to determine the normality of the data distribution. In occasions when the data is not normally distributed as commonly found among environmental data, the data undergoes transformation (commonly applied is the log transformation) to make the data log-normal [32].

It is important to note that for a spatial distribution of data to undergo ordinary kriging interpolation, the data has to meet the assumption of normality and stationarity (constant variance), to provide accurate predictions. Data transformation is a quick fix to the issue of normality assumption, and also an applicable solution during instances when the data fails to meet the stationarity assumption. The Log transformation of skewed data sets during the exploratory data analysis draws the data closer to normality and satisfaction of the stationarity assumptions [54]. In order to minimize the deviation from the assumptions and/or avoid the violation of these assumptions, it is therefore important for skewed data to be transformed. The primary basis for analyzing the trend during geostatistical interpolation is to satisfy stationarity assumptions [53]. Since the log transformation of the data fixes the issue of normality assumption and satisfaction of stationarity assumption, it is inferred that there is no trend in the data distribution and the stationarity assumption is satisfied. Johnston et al. [53] and Chiu [56] indicate that the Voronoi map provides a deeper understanding of the distribution patterns of points and examines the trend and stationarity assumption.

2.2.3. Experimental Semivariogram

The Kriging method is based on the principle of spatial autocorrelation, which is central to geostatistical applications. To analyze the spatial distribution of data values and their parameters, semivariograms serve as a key descriptive tool. Spatial dependence between measured points is determined by calculating semivariance, which considers the distances between these points. The use of plotted semi-variogram indicates the relationship between the lag distance and the semi-variance, the variogram evaluates the sampling errors [35]. The lag distance refers to the separation between measurements of a given property. Typically, the semivariance values increase with decreasing lag distance, indicating stronger spatial autocorrelation [57]. Various models, such as Gaussian, Spherical, and Exponential, are then applied to fit the semivariogram and identify the most appropriate model for generating optimal interpolation weights [52]. Kriging is a highly versatile interpolation method that can produce both exact and smooth surfaces. It supports a range of output surfaces, including predictions, standard errors of prediction, and probabilities. The Kriging technique facilitates optimal, unbiased estimation of regionalized variables at unsampled locations by leveraging the properties of semivariograms and initial data points [35].

In geostatistics, the semivariance function is utilized to calculate and interpret spatial correlation patterns, particularly when the sampled data values follow a normal distribution, and the skewness falls within the range of -1 to +1 [35]. Semivariogram analysis is a key technique for visually indicating spatial correlations among neighboring data points. A semivariogram plot is generated by calculating semivariance values across various lags. The models—circular, spherical, exponential, and Gaussian—furnish insights into the spatial structure for Kriging interpolation. Key characteristics of the semivariogram (as shown in Figure 3) include range (a), sill (C0 + C), and nugget effect (C0) [32]. The nugget refers to the Y-intercept of the semivariogram model, while the sill is the semivariogram value at which the model stabilizes or levels off [58]. The selection of the appropriate lag size is calculated based on multiplying the number of lags times the lag size, which should be ‘half the farthest distance between any point’ [53] or about the ‘whole farthest distance between any point,’ as applied in this study. However, the range of the fitted semivariogram model is highly influenced by the lag size [53]. The equations below indicate how the semi-variance function, γ(h) and unmeasured values at a given distance, h are calculated [59,60]:

$$\mathsf{\gamma }\left(\mathrm{h}\right)={\displaystyle \frac{1}{2n}}{\sum }_{i=1}^n{(Z\left(x_i\right)-Z(x_i+h\left)\right)}^2$$

(1)

where: γ(h) = semi variance value at lag distance, h,

Z = regionalized variable,

Z ($x_i$) = measured value at point, $x_i$,

Z($x_i+h)$ = measured value at point, ($x_i+h)$,

n = number of pairs separated by distance, h.

$$Z\left(X_o\right)={\sum }_{i=1}^n{\lambda }_{i}Z\left(x_i\right)$$

(2)

where: Z (X_o) = estimate of unknown true value,

λ_i = weighted coefficient, n = number of neighboring observations.

2.2.4. Cross-Validation

Experimental semivariogram values were tested for different types of semivariogram models (Circular, Spherical, Exponential, and Gaussian) and were evaluated to identify the most suitable model for salinity estimation. The predictive performance of each fitted model is assessed using spatial cross-validation metrics. Various error metrics are calculated, these include “mean error (ME), mean square error (MSE), root mean square error (RMSE), average standard error (ASE), and root mean square standardized error (RMSSE).” Ideally, if the predictions are unbiased, the ME should be close to zero, though this statistical inference has limitations, as it is scale-dependent and does not fully capture inaccuracies in the variogram. Therefore, MSE, which builds on ME by squaring the errors, is often used, with an ideal MSE value approaching zero, indicating an accurate model [32].

As MSE should be close to zero to ensure prediction accuracy, RMSE also should be as minimal as possible [29,51,53]. If RMSE aligns with ASE and their values are close, it indicates that the prediction errors were correctly estimated. However, if RMSE value is smaller than ASE, the prediction is overestimated, while an RMSE greater than ASE suggests that the prediction is underestimated. Similarly, RMSSE values close to one indicate correct estimation of variability. Values above one suggests underestimation, while values below one imply overestimation of variability. After completing the cross-validation process, Kriging-generated maps were produced, providing a visual depiction of the spatial distribution of groundwater quality [31].

2.2.5. Kriging Map Generation

After cross-validation, the estimated kriging values of the groundwater for all sampling cycles are mapped out, and the maps are graded according to the level of salinity as indicated by USSL [26]. Available freshwater quality of the WMAs in northern New Jersey is evaluated, as the study uses 750 μS/cm value of specific conductance as the threshold for salinity, and any value higher than 750 μS/cm is not considered to be freshwater. Total area with available freshwater is assessed for the four sampling cycle periods – indicating the freshwater coverage and changes in freshwater per the study area coverage over time. Thereafter, layers of WMAs, major roads and highways in northern New Jersey are overlaid on the generated kriging maps and are assessed to observe their relationship and influences in the varied levels of groundwater quality across the region.

2.2.6. Limitation

It is important to recognize that actual spatial variations can significantly differ from the values predicted by spatial interpolation, which may present limitations in the Kriging method, particularly when data is sparse or unevenly distributed. Therefore, understanding both the number of data points and the geographical coverage of the region is crucial. The kriging interpolations are executed based on two primary assumptions: the stationarity of the target variable and its conformity to a normal distribution. Nevertheless, environmental data often deviate from these assumptions, exhibiting complex patterns and distributions. In response to this challenge, this study undertakes a data transformation to attempt to satisfy the assumptions [52]. One of the most challenging aspects in such cases is accurately estimating the variogram model, which becomes more difficult with fewer data points [31]. Additionally, ME statistical metric used for evaluating prediction performance has notable limitations, as it is scale-dependent and does not fully account for variogram inaccuracies. To address this, MSE which builds on ME by squaring the errors, is often used, with an ideal MSE value close to zero indicating an accurate model [32].

3. Results and Discussion

Rapid urbanization and industrialization across the world are degrading freshwater resources, as anthropogenic activities continue to alter land use and introduce contaminants into the environment. Kriging interpolation provides visual insight into the state of groundwater quality, clearly indicating the freshwater coverage in any region. Groundwater specific conductance value of 750 μS/cm or lesser is used as a marker for freshwater.

3.1. Exploratory Data Analysis

First an exploratory analysis process, which is crucial for spatial relationship of the data and appropriate selection of model is performed for all sampling cycles. The histogram in Figure 4a reveals that the groundwater data is highly skewed, does not have a bell-shape, and therefore, unsuitable for kriging interpolation. Figure 5a shows the use of normal Q-Q plots to provide an analysis of univariate normality. Data for all the sampling cycles indicated to be asymmetric, and most of the points are drifting from the line of best fit. The Normal Q-Q plots and histograms indicate that the data is not normally distributed and needs to be transformed before further analysis. After the log transformation of the data, the histograms and normal Q-Q plots show that the distributions are close to normal (Figure 4b and 5b). The log normal statistics (Table 1) reveal close values of the means and medians, the skewness is near zero (falls within the range of -1 to +1), and the kurtosis values are close to 3. These statistics of the log transformed data show that the data are normalized and meet at least one of the assumptions for kriging interpolation – the normality assumption.

This log transformation ensures equal variance across the data, thereby possibly satisfying the other assumption (stationarity). The Voronoi maps (Figure 6) of the transformed data provide a visual and qualitative assessment of trends, as well as an analysis of the spatial variability of the data [56,61]. The Voronoi maps indicate that there is no clear existence of linear global trend in the data and provide satisfaction towards the stationarity assumption. Figure 6 has an Illustration of trend occurrence in spatial data, as the Voronoi map shows the occurrence of a visually noticeable linear trend.

The removal of trend from data during kriging interpolation should be based on a proper justification [53]. When spatial coordinates indicate an observed or suggested trend in the data [49], as seen in this study where land use significantly influences water quality, it would make sense to retain the trend.

Furthermore, Jarmołowski et al. [62], in their study on "drawbacks of trend removal during kriging interpolation," provide rationales for retaining trends in heterogeneous spatial data, such as data significantly influenced by land use type. Their rationales include (1) maintaining the trend of a spatial data can help capture the underlying spatial structure more effectively, (2) unbiasedness of ordinary kriging must be ensured, and retaining the global trend can help achieve this by providing a more stable reference for predictions, (3) In scenarios where data is sparse, keeping the trend can provide a more reliable framework for interpolation. Following the logarithmic transformation of the data in this study, a deliberate decision was made to refrain from detrending the data. This approach was adopted to minimize complexity and ensure reliability of the results. Considering that groundwater quality in this study is influenced by varied land use patterns and non-point pollution sources, which are spatially heterogeneous and linked to specific locations, detrending may mask these spatial relationships. Therefore, not removing the trend (if any) in the data is justified. This study relies on the predictive performance metrics of the kriging interpolation to measure the performance of the kriging interpolation.

3.2. Experimental Semi-Variogram

The plotted semi-variograms (Figure 7) show the relationship between the lag distance and the semi-variance. The variogram evaluates the sampling errors across the spatial data for all the sampling cycles. The selection of the appropriate lag size for this study is evaluated based on multiplying the number of lags by the lag size, which is about the ‘whole farthest distance between any point.’ This is in contrast with the approach suggested by Johnston et al. [53], that the ‘golden rule’ for the selection of lag size is based on half of the “farthest distance between any points.” Applying ‘about the whole distance’ instead of ‘half the distance’ in this study (410,000 – 450,000 meters) showed to support more prediction accuracy as reflected in the prediction performance metrics (MSE, RSME, etc.). Using the established lag size and lag number, the experimental semivariograms are fitted for various types of semivariogram models (Circular, Spherical, Exponential, and Gaussian). Nugget and sill parameters were calculated using the Geostatistical Analyst tool for each model, and no adjustments were deemed necessary. The nugget represents the Y-intercept of the semivariogram model, while the sill refers to the semivariogram value at which the model stabilizes or levels off [58]. The gaussian model showed to be the best fit for all the sampling cycles, and the statistics metrics (Table 2) were used to evaluate the performance of each model.

3.3. Cross-Validation

Cross-Validation facilitates model evaluation by estimating prediction accuracy at unmeasured locations, with ideal model performance characterized by near zero MSE, minimal RMSE, standardized RMSE near 1 and equivalent values of RMSE and ASE. Table 2 shows the model performance metrics for all the sampling cycles, indicating that the selected models are of optimal performance level. The metrics show that RMSE and ASE have equivalent values, both close to zero, while MSE and RMSSE are near 1, suggesting that the model was neither overestimated nor underestimated.

3.4. Kriging Map Generation

After evaluating the model performance metrics, the estimated kriging values of specific conductance of the groundwater for all sampling cycles are mapped out, and the maps are graded according to the level of salinity as indicated by USSL [26]. The kriging interpolation produced a raster map of salinity indicating various levels of groundwater quality. The kriging interpolation provided a salinity distribution representing the groundwater quality of northern New Jersey WMAs, from sampling cycling 1 to 4, plus sampling cycle 4 – adjusted. This study uses 750 μS/cm value of specific conductance as the threshold for salinity, and any water with salinity higher than 750 μS/cm is not considered to be freshwater. The kriged maps (Figure 9) provide an insightful visualization of the groundwater quality various across the land use types in the study area. The legend of the kriged map shows color patterns from blue to red, with blue indicating the lowest levels (best water quality), and red the highest levels of salinity. Several layers were overlaid on the kriged surface to provide better understanding with varied influences from anthropogenic activities. These layers include highways, all major roads, sampling points, and 11 WMAs of northern New Jersey.

The major highway and WMA layers are overlaid onto the land use layer map (Figure 1b) to visually illustrate the region’s urbanized area and their relationship with varied levels of groundwater quality. The highway overlay suggests potential contributions to pollution, from major highways, likely from plying vehicle and the application of deicing salts during icy and snowy conditions. The land use map revealed that the Northeast water region is the most urbanized, followed by Raritan water region. The Northwest region has the most coverage of undeveloped and agricultural land use areas. The urbanized Northeast region is characterized by a dense network of roads, including highways and major roads. When these roads were overlaid on the kriged map, they nearly covered the entire Northeast region, making it difficult to clearly recognize the patterns (colors) on the kriged map beneath (as shown in Figure 8). This high road density highlights the prevalence of impermeable surfaces, which serve as areas where contaminants accumulate. During precipitation events, runoff carries these contaminants to areas not originally associated with such pollution, exacerbating water quality issues linked to non-point pollution sources. Also, the larger the impervious surfaces (such as roads), the larger the volume of deicing salts applied during snowy and icy conditions, indicating that areas with heavily dense roads will have their water resources more impacted by the applied deicing salts. The overlaid sampling points reveal that the sampling points are sparse and evenly distributed, covering an area of about 3,400 square miles. The sampling points are considered sparse because certain WMAs have very few points such as WMA 03 – Pompton WMA in the Northeast region. This WMA has only two wells (in an area of about 240 square miles) for monitoring the ambient water quality. The Pompton WMA is mostly undeveloped forested area (as indicated in the land use map), that is observed to have better water quality relative to other land use areas. That the location of WMA is mostly forested area, may be one reason the AGWQMN monitoring network has only two wells. However, with the rapid urbanization occurring as indicated by Lathrop [28], more wells may need to be developed in such regions so as to capture any possible water quality deterioration that could result from anthropogenic activities. The rapid change in land use towards urbanization could cause the groundwater to deteriorate in quality.

Figure 9 provides a kriging interpolation prediction map of groundwater quality in northern New Jersey across four ambient groundwater sampling cycles (1999 – 2004, 2004 – 2008, 2009 – 2013 and 2014 – 2016). The map provides a prediction of specific conductance values, which serves as an indicator of salinity levels, and an overall metric for the quality of the groundwater in the region. The map aims to quantify the coverage of available freshwater in northern New Jersey, while also highlighting the WMAs that are exhibiting increasing salinity and/or characterized as ‘non-freshwater’. On the kriged map, colors ranging from blue to yellow indicate freshwater. The yellow color indicates a specific conductance value of 750 μS/cm, which serves as the threshold for freshwater in this study. In contrast, brownish or reddish areas on the kriged map indicate high salinity levels, reflecting poor groundwater quality. As shown in Figure 9a, the spatial distribution of salinity reveals a distinct pattern. The Northwest water region, predominantly blue, exhibits low salinity, while the Raritan region shows low to moderate salinity (ranging from green to yellow). In contrast, the Northeast region displays the highest salinity, predominantly red, during the first sampling cycle. The distinctive levels of kriged salinity in the study area is highly associated with the varied land use of the region. The most urbanized area exhibits the worst water quality, while regions with less urbanization and more undeveloped land use showed better water quality with respect to salinity. This expression of salinity is also associated with the density of road network across the study area. The Northeast water region, which has the most densely concentrated road network, exhibited the highest salinity. This is likely attributed to the widespread use of deicing salts, facilitated by the region's large areas of impermeable surfaces [5,63] Aside from the influences of deicing salts, various studies [5,63,64] have revealed that wastewater contributes to increasing salt ion levels in water resources. This heavily urbanized Northeast water region is serviced by old sewer collection systems (combined sewer system) that combine sanitary sewer waste and stormwater in one line. During heavy storm events, overflow through the manholes of these combine sewer systems could bring up heavy salt-laden wastewater to the subsurface, thereby infiltrating into the groundwater system and contaminating water resources.

Figure 9b illustrates the kriging prediction for the second sampling cycle. The spatial distribution of salinity across the study area closely resembles that of the first sampling cycle, with the northeast water region again showing the highest salinity. This is consistent with the region's heavy urbanization, as previously discussed. However, the kriged map reveals an overall increase in salinity across all water regions. Notably, areas that were previously blue, green, or even brown color in the first sampling cycle now exhibit more yellow coloration, indicating fluctuating salinity levels across the region. The Northwest region shows more of the yellow color in this sampling cycle. The increasing salinity in this region may be attributed to use of agricultural fertilizer, an activity that is prominent in the Northwest region’s land use. Kaushal et al. [5] and Marghade et al. [64] indicate that agrochemicals contribute to freshwater salinization, however, the Northwest region still has a good coverage of freshwater compared to other water regions.

During the third sampling cycle, the kriged map (Figure 9c) reveals that the spatial distribution of the salinity across the study area remains similar to those of the two previous sampling cycles. The Northeast water region has a larger coverage of poor groundwater quality compared to other water regions in northern New Jersey. The water quality in the water regions were chiefly influenced by their varied land use activities. There is one visually observable difference in sampling cycle 3 within the Northwest region. In this region, there appears to be a reduction in salinity, as indicated by areas that were yellow during cycle 2, now turning green in cycle 3, suggesting a decrease in salinity. For instance, an area with a specific conductance value of 750 µS/cm (yellow) during cycle 2 may now measure 700 µS/cm (green) in cycle 3. While this change demonstrates a decline in salinity, it does not necessarily imply an increase in freshwater coverage. This is because, from the cycle 2 to cycle 3, any area with specific conductance value of 750 µS/cm or less is considered to be freshwater. Thus, the observed reduction in salinity is not accompanied by clear evidence of expanded freshwater availability in the region. The reduced salinity in the Northwest region could be attributed to adequate freshwater flushing in the groundwater system, as reported by Rossi et al. [63] and Sun et al. [65]. Their studies stated that adequate freshwater flushing of the groundwater system from rainfall can influence groundwater chemistry by diluting the contaminated effluent entering the groundwater system.

Figure 9. Northern New Jersey Salinity for all Sampling Cycles.

Figure 9. Northern New Jersey Salinity for all Sampling Cycles.

The kriged prediction in sampling cycle 4 (Figure 9d) is unique, showing a concentration pattern that is different from the earlier three sampling cycles. The only noticeable similarity with the previous sampling cycles is that anthropogenic activities have a major influence on the groundwater quality. The kriged map shows that the Northeast is still the region with the most coverage of poor water quality. The uniqueness of sampling cycle 4 is associated with the number of observation wells monitored during this sampling period. This sampling cycle did not include 12 observation wells (about 15 % of the total observations) that were used in the previous three sampling cycles. Sampling cycle 4 had only 65 observation wells as compared to 77 wells in previous sampling cycles. These 12 wells are located within the undeveloped land use (forested) areas that have been shown to have better water quality compared to the other land uses. Nine of these wells are located in the upper part of the study area and two are in the lower part. The omitted wells from the undeveloped land use area in the region usually have low specific conductance values and some of them may be regarded as outliers. As indicated by Johnston et al. [53], such outliers can highly influence kriging estimation values. The kriged prediction map reveals that the Northwest region which had the best groundwater quality coverage in the previous sampling cycles, now shows increased salinity, although most of the specific conductance values are not beyond the 750 μS/cm threshold for freshwater. The omitted observation wells highly influenced the overall spatial autocorrelation of the study area, and this has a major influence on the kriging prediction. The kriging interpolation is not based entirely on the distance between the observed and the estimated points like in other spatial interpolations, but it is based mostly on the weighted average of all the observed points [53]. The kriged map looks different because of the dissimilarity between the measured points and the overall spatial arrangement in the previous sampling cycles. It is known that additional data enhances the performance of kriging predictions. The kriged maps of sampling cycles 1, 2, and 3, showed minimal changes of salinity in the ‘undeveloped (forested) land use’ areas over time. Assuming that this minimal change in salinity will continue into the 4th cycle, it was decided in this study to incorporate the ‘undeveloped land use’ well data from sampling cycle 3 into the sampling cycle 4. This modified sampling cycle is referred to as the ‘Adjusted sampling cycle 4’.

The kriged map from the ‘Adjusted sampling cycle 4’(Figure 9e) shows that spatial distribution of groundwater quality across northern New Jersey is highly related to land use and rapid urbanization. The analysis indicates a general rise in groundwater salinity in the Northeast and Raritan water regions, while the northwest region experienced a decline in salinity. Notably, the Raritan water region displayed a notable increase in salinity for the first time in the adjusted cycle 4. Similarly, the Northeast water region saw some increase in groundwater salinity, making it the water region with the poorest water quality during all the sampling periods.

This study quantifies the coverage of freshwater across the WMAs in northern New Jersey. All the WMAs in northern New Jersey have about 3,400 square miles coverage representing the total coverage of groundwater in the study area. Quantifying the extent of the freshwater coverage in the region enhances the understanding of available freshwater resources, offering insights that go beyond the visual interpretation of the kriged maps. To determine the coverage of freshwater from the kriged map, the raster surfaces were converted to polygons using ArcGIS geostatistical analytical tool and the area of polygons with 750 μS/cm and lesser were calculated to determine the area of freshwater in the study area. During sampling cycle 1 (1999 – 2004), about 80 % of the total coverage of the study area was freshwater. Non-freshwater coverage was observed in sections of WMA – 03 (Pompton), WMA – 04 (lower Passaic), WMA – 05 (Hackensack), WMA – 06 (upper Passaic), and WMA – 07 (Authur Kills). Among these, WMA-04 (Lower Passaic) and WMA-05 (Hackensack) exhibited the highest levels of salinity impact, marking them as the most affected WMAs. Sampling cycle 2 (2004 – 2008) had more coverage of freshwater, about 82.5 % of the total coverage of the study area was freshwater. This is an increase of about 2.5% coverage of freshwater in the study area. Similar to sampling cycle 1, in cycle 2, non-freshwater coverage was also observed in parts of WMA – 03 (Pompton), WMA – 04 (lower Passaic), WMA – 05 (Hackensack), WMA – 06 (upper Passaic), and WMA – 07 (Arthur Kills), with poorer groundwater quality observed in WMA – 04 (lower Passaic) and WMA – 05 (Hackensack). Sampling cycle 3 (2009 – 2013) has similar freshwater coverage as cycle 2, with about 82.4 % of the total coverage of the study area indicating freshwater. Non-freshwater coverage was observed in WMA – 03 (Pompton), WMA – 04 (lower Passaic), WMA – 05 (Hackensack), WMA – 06 (upper Passaic), and WMA – 07 (Arthur Kills), with the poorest groundwater quality observed in WMA – 04 (lower Passaic) and WMA – 05 (Hackensack). WMA – 03 (Pompton) had more freshwater coverage during sampling cycle 3, while there was a decrease in freshwater coverage in Arthur Kills WMA in the same sampling cycle. Sampling cycle 4 (2014 – 2016), which has sparse distribution and influenced by the data’s spatial autocorrelation indicated that about 51% of the total coverage of the study area had freshwater. This is about 30% less freshwater coverage compared to sampling cycle 3. Though statistics of the sampling cycle 4 indicated the highest ‘minimal, maximum and average” value of specific conductance, the spatial prediction of the data was heavily influenced by unavailability of important undeveloped land use data. This concern brought about the ‘adjusted sampling cycle 4’ which includes undeveloped land use data from cycle 3. For the ‘adjusted sampling cycle 4’, the freshwater coverage was about 74%, resulting in a decrease of about 8% in the groundwater quality when compared to sampling cycle 3. Non-freshwater coverages were observed in parts of WMA – 03 (Pompton), WMA – 04 (lower Passaic), WMA – 05 (Hackensack), WMA – 06 (upper Passaic), WMA – 07 (Arthur Kills), WMA – 08 (north Branch Raritan), and WMA – 10 (Millstone), with poorer groundwater quality observed in WMA – 04 (lower Passaic) and WMA – 05 (Hackensack) while WMA – 08, and WMA – 10, indicated a generalized increase in salinity. The kriging interpolation reveals a disturbing trend after four sampling cycles: a steady decline in freshwater coverage. The Northeast water region particularly WMA – 04 (lower Passaic) and WMA – 05 (Hackensack) experienced the most significant impact of freshwater depletion.

4. Conclusions

Kriging interpolation technique proves effective and invaluable for groundwater resources management. This study employed kriging interpolation to spatially analyze groundwater in WMAs across northern New Jersey. Kriging effectively modeled salinity distribution by developing models for the prediction of specific conductance values in locations lacking direct observation. Kriging also facilitated the identification of fresh groundwater coverage in the study area. The spatial and temporal distribution of groundwater quality was assessed.

The findings provide critical insights into freshwater availability and groundwater quality trends over four sampling cycles (1999–2016). Using kriging and GIS tools, the study quantified freshwater distribution, revealing that the Northeast water region had significantly less freshwater coverage compared to the Northwest and Raritan regions. Specific problem areas, such as WMA-04 (Lower Passaic) and WMA-05 (Hackensack) in the Northeast region, consistently exhibited poor groundwater quality due to high salinity, highlighting the need for good water management practices in the area.

The study reveals that anthropogenic activities have had a significant impact on groundwater quality, particularly in the Northeast water region. The spatial analysis identified strong influences of road network density and land use activities on the groundwater quality. The rapid change in land use towards urbanization and the application of deicing salts on the large surface area of the road network appears to be a major control of the deteriorating groundwater quality.

This spatial investigation of groundwater quality in northern New Jersey provides a great insight into groundwater quality in the region. It serves as a beneficial guide to decision-makers and stakeholders, highlighting the need for effective water management strategies in the region. This study is a prelude to more research in the study area. A regionalized and localized study is needed to provide a better understanding of the poor groundwater quality revealed in parts of the northeast water region. It is pertinent to know whether the groundwater contamination is entirely localized or regional. A comprehensive groundwater transport study could indicate if salinity contributions in the Northeast water region originate from activities from surrounding watersheds. It will also be important to know whether the rivers (surface-water bodies) around the watershed are saline, and whether they are losing water into and contaminating the surrounding groundwater.