Exploiting Chloride Conservative Tendencies as Contaminant Surrogates in Groundwater Transport Modeling in a Typical Hydrogeological Environment of Northern New Jersey

1. Introduction

Urbanization and industrialization in New Jersey have come at a cost: the steady decline of water quality across the state. As land use changes and human activities intensify, groundwater resources face growing contamination from both natural processes and essential societal practices [1]. This trend mirrors a national challenge, with studies documenting widespread groundwater degradation [2,3,4,5].

Two major contributors stand out. Agriculture, while critical for food production, introduces fertilizers that often seep into aquifers as harmful contaminants. Equally significant—and less visible—is the use of chloride-based deicing salts during winter months. These salts, though vital for road safety, infiltrate groundwater systems, driving salinization and threatening freshwater supplies [2,6]. These activities—though essential for everyday life—are not without consequences. Natural events, such as heavy precipitation and flooding, compound the problem by dispersing contaminants across vast areas. In urbanized regions like New Jersey, the risks are magnified. Aging combined sewer systems, which merge stormwater and municipal sewage, often overflow during extreme flooding events, releasing untreated effluents through manholes into critical water resources. These discharges inevitably seep into groundwater, exacerbating contamination [2].

Together, these scenarios highlight the dual role of point and non-point pollution sources. Among them, deicing salt applications emerge as a dominant driver of groundwater salinization [2,6,7]. Supporting this, Oyen and Ophori [8] found that New Jersey’s Northeast water region—particularly areas with dense road networks—experienced disproportionately high freshwater depletion and salinity spikes, directly linking winter salt use to groundwater degradation.

Despite its classification as a secondary contaminant, chloride is the primary driver behind groundwater salinity [6,9,10,11]. Notably, chloride concentrations often surpass the EPA-recommended 250 mg/L threshold for drinking water [9,12], with several studies documenting widespread exceedances [2,7,13]. Beyond salinity itself, elevated chloride levels can mobilize heavy metals [14,15] and hinder remediation of organic compounds such as Trichloroethylene (TCE) through salting-out effects [16,17] — a process where high ionic strength reduces contaminant solubility and prolongs their persistence in aquifers.

Elevated groundwater salinity does more than degrade water quality—it actively reshapes contaminant dynamics. Deicing salts introduce sodium ions (Na⁺) that trigger cation exchange in soils, displacing heavy metals like lead and zinc. These metals, once bound to soil particles, become mobile and migrate into both surface waters and groundwater [14,15]. The impacts extend beyond metals. McNaboe et al. [18] documented how deicing salts in Connecticut drive Na⁺ infiltration deep into soils, mobilizing radium (Ra) and radon (Rn) upward. Their studies reveal that highly saline groundwater systems are most likely to have Ra concentrations four times the recommended standard concentration required in groundwater—a clear health risk for groundwater-dependent communities.

This phenomenon is not isolated. Across the United States, deicing salts and urban sprawl have become combined drivers for groundwater salinization [10,19,20]. The ecological toll is equally pronounced: rising salinity favors invasive species over native ecosystems, disrupting habitats and biodiversity [2,19,21].

Compounding these risks is the “salting-out” effect—a process where elevated ionic strength reduces the solubility of dissolved constituents in groundwater. In principle, salting-out processes may effectively retain per-and polyfluoroalkyl substances (PFAs) fluids in deep waste disposal wells, where saline conditions and depth may isolate them from drinking water sources [22]. However, Newell et al. [22] warn that salting-out also enhances PFAs retention in saline groundwater, complicating remediation—especially in fractured systems like the Brunswick aquifer, where contaminants can persist for extended periods. In fresh groundwater with high salinity, PFAs can become persistent in the environment. Reduced solubility and heightened sorption to aquifer materials allow these contaminants to linger. Even after remediation, PFAS trapped in fracture pores can diffuse back into groundwater—a cyclic recontamination scenario that undermines cleanup efforts [16,17].

Deicing salts application represent a paradox: a measure for public safety that doubles as a silent environmental threat. Their pervasive use demands urgent scrutiny, particularly given their role in groundwater salinization and co-contaminant mobilization. Groundwater modeling offers a powerful solution to track the salt constituent (chloride) transport, with three key advantages: (1) provide a deeper understanding of advective chloride transport, (2) spatial quantification of dispersed chloride concentrations, and (3) enable chloride to serve as a surrogate for contaminants. By treating chloride as a conservative tracer, models can predict the fate of broader contaminant plumes. Groundwater modeling tools such as MODFLOW and MT3DMS are widely used for simulating groundwater systems, capable of reconstructing flow paths and contaminant dispersion. Their outputs inform critical decisions—from resource allocation to remediation strategies—across all phases of groundwater management.

These models translate subsurface processes into actionable insights, enabling data-driven resource management [13,23,24]. The model validation relies on the calibration of hydraulic parameters (such as heads, recharge) to match observed conditions—a process demonstrated as robust across diverse studies [23,24]. Statistical metrics like root mean square error (RMSE) and mean error (ME) further validate model performance, where lower values generally indicate higher accuracy. However, even models with optimal RMSE/ME may not be appropriate if calibrated parameters (such as hydraulic conductivity, recharge, etc.) deviate from actual field observations. Thus, researcher judgment remains essential to ensure realistic simulations [23].

Groundwater modeling has become an essential tool for characterizing aquifers across diverse geological settings, with particularly valuable applications across New Jersey [13,25,26,27]. These studies have significantly advanced our understanding of groundwater systems, especially in northern New Jersey, where complex hydrogeological features present unique management challenges compared to the south. Barry et al. [25] used the MODFLOW code to simulate groundwater flow in the Central Passaic River Basin (CPRB) of New Jersey, generating detailed recharge maps that revealed critical insights about transient capture zones in historically exploited aquifers. Building on this work, Ophori et al. [13] developed numerical models to examine groundwater–surface water interactions during historical pumping regimes in the CPRB. Meanwhile, Serfes et al. [27] employed chemical equilibrium modeling to investigate Arsenic(III) mobilization from pyrite oxidation in the Lockatong Formation, confirming earlier propositions about this contamination source. Lewis-Brown et al. [26] specifically targeted fractured-rock systems, using advanced modeling techniques to precisely quantify contaminant velocities, travel times, and flow paths in the Lockatong Formation. Collectively, these studies demonstrate how groundwater modeling can be adapted to address New Jersey’s diverse hydrogeological challenges, from resource management to contaminant transport prediction.

While prior studies [13,25] have successfully modeled groundwater flow in northern New Jersey’s aquifers, they predominantly treated fracture networks as homogeneous conduits, overlooking the dual-domain dynamics that govern contaminant retention and release. Our work bridges this gap by integrating fracture-matrix mass transfer—a process critical for predicting long-term salinity and co-migrating contaminant persistence—demonstrating how contrasting transport mechanisms (fracture networks vs. matrix diffusion) govern contaminant (chloride) migration disparities in urbanized fractured aquifers. This dual behavior, previously unquantified in the Brunswick aquifer, redefines remediation strategies for urban fractured systems.

This study develops a dual-porosity model to simulate chloride transport in the fractured Brunswick aquifer, a system characteristic of northern New Jersey's hydrogeological complexity. Using MODFLOW and MT3DMS codes, we examine (1) the dominance of fracture networks in short-term chloride migration to local streams, and (2) the role of matrix diffusion in sustaining long-term regional contamination. By resolving these processes, our work provides a framework for assessing salinization risks in similar hydrogeologic settings worldwide.

1.1. Study Area

The 105 km² study area encompasses part of the Brunswick aquifer in New Jersey's Newark Basin, located within the Lower Passaic watershed. The area is delineated by topographic divides and surface water boundaries, with groundwater flow converging toward local flow systems (e.g., Peckman River) and a regional discharge area (Passaic River). The area is bounded by the Peckman River (west), Passaic River (north and east), and the upper Montclair community (south). As shown in Figure 1, the predominantly urbanized study area includes Montclair Township, home to a commuter university, and extends into parts of Clifton, Paterson, Cedar Grove, Little Falls, Glen Ridge, Bloomfield, Belleville, and Nutley. The area covers parts of Essex and Passaic counties. Montclair State University (MSU), a commuter institution, serves over 20,000 students. The region is also supplied with water by three municipal wells; the Glenfield, Lorraine, and Rand wells located in Montclair Township, which provide water to Montclair and surrounding communities. Two of the three wells, the Lorraine well and Rand well, are situated within the study area. The study area represents a characteristic example of northern New Jersey's typical hydrogeological environment, having the region's key hydrogeological features. In this study, the MSU area, which lies on parts of the First Watchung Mountain, is identified as an area of critical importance. This designation results from its topographically high elevation and the daily gathering and commuting of large groups of individuals. Such factors are considered to intensify environmentally detrimental impacts.

While large gatherings and frequent commuting are often linked to increased environmental degradation, substantial efforts and actions are undertaken to ensure safe conditions within the university community. However, some of these actions, such as the application of deicing salts during wintry conditions to ensure safety and normalcy, may inadvertently harm the environment. These deicing salts, primarily composed of chloride compounds, infiltrate into groundwater, leading to increased groundwater salinity. These unique characteristics of the university community, along with other factors, influenced the selection of the study area. Another critical factor is the presence of three public supply wells around the study area, one of which has been identified as contaminated. A 2022 water report rated all three wells as highly susceptible to volatile organic compounds (VOCs), radionuclides, inorganic contaminants, and other pollutants [28]. Currently, the Glenfield and Lorraine wells operate with carbon absorbers, while the Rand well is not in operation. According to Staff [29], Federal funding was secured by the Montclair Township to rehabilitate the Rand municipal supply well, aiming to address concerns linked with elevated levels of perfluorooctane sulfonate (PFOS) and perchlorate contaminants. This planned rehabilitation may explain why the Montclair Water Bureau [28] reports that the Rand well is currently inactive.

1.2. Geologic Setting

The area is underlain by the Brunswick aquifer (Early Mesozoic), a hydrogeologic unit within New Jersey's Newark Basin. The aquifer features a complex geological history and consists primarily of consolidated sedimentary and igneous formations dating to the Late Triassic and Early Jurassic periods, including the Passaic Formation, Orange Mountain Basalt, and Feltville Formation. This diverse stratigraphy makes the Brunswick aquifer a true microcosm of northern New Jersey’s geology. The Passaic Formation is predominantly composed of red mudstone and siltstone, with smaller amounts of gray and black shale. It is distinguished by its fine-grained facies, which transition laterally into sandstone and conglomerate along the basin's margins. This formation can attain a maximum thickness of approximately 11,800 feet, making it one of the most substantial units in the region. The Orange Mountain Basalt consists of dark greenish-gray to greenish-black basalt and is part of a series of interlayered basalt and clastic sedimentary rocks from the early Jurassic period, with layers reaching up to 597 feet in thickness. It includes multiple basalt flows fed by diabase sheets and dikes, which intruded and thermally transformed older sedimentary deposits. The Feltville Formation, part of a sequence of sedimentary rocks deposited in environments such as fluvial-deltaic settings, is interlayered with basalt flows, reflecting a dynamic geological history. This formation consists of interbedded brownish-red to light-grayish-red, fine-to-coarse-grained sandstone, as well as gray and black coarse siltstone and silty mudstone, with a thickness of about 510 feet [30,31,32].

1.3. Hydrogeological Setting

The Brunswick aquifer, which forms part of the fractured rock aquifer system in the Newark Basin, consists of sandstone, conglomerate, siltstone, and shale. Fractures within the aquifer serve as both storage spaces and pathways for groundwater, which is typically fresh, non-corrosive, and dominated by calcium-bicarbonate. The aquifer has a median yield ranging from 100 to 250 gallons per minute (gpm) [33]. It exhibits heterogeneous properties, with significant spatial variations in permeability and transmissivity across its layers. Depths below 320 feet generally show low permeability (transmissivity less than 1 gpd/ft), while shallower formations present high transmissivity (up to 30,000 gpd/ft). Groundwater flow within the aquifer is primarily driven by local topographic gradients and the stratigraphic configuration of its units [31].

Within the Brunswick aquifer, groundwater flows through fractures across various strata [27]. The aquifer's hydrogeological characteristics are fundamentally shaped by its complex fracture network, which comprises several types of water-bearing features (WBFs). These WBFs include bedding planes, fracture planes, and linear intersections between them.

The bedding planes formed along the sedimentary rock layers are common throughout the Brunswick aquifer, where they serve as important conduits for groundwater movement. Fracture planes, typically oriented vertically or near-vertically, intersect with these bedding planes and enhance overall aquifer permeability by facilitating increased groundwater flow [31]. This fracture-dominated flow system exemplifies the hydrogeologic complexity of northern New Jersey's Brunswick aquifer, where heterogeneous permeability governs contaminant transport pathways. Urban stressors like road salts amplify contamination risks through dual-domain dynamics - where rapid groundwater transport through fractures coexists with matrix diffusion into low-permeability rock. These processes are particularly pronounced in recharge zones like the First Watchung Mountain area, where anthropogenic inputs enter the groundwater system.

Linear intersections occur where bedding planes and fractured planes cross, creating additional pathways that further promote groundwater movement. These intersections significantly enhance hydraulic connectivity within the aquifer system [31]. The aquifer's groundwater storage occurs predominantly within its extensive fracture networks, while long-distance transport is made possible through bedding planes and highly permeable zones located between confining layers. This dual mechanism of storage and transport, mediated by the diverse water-bearing features, defines the Brunswick aquifer's unique hydrogeological characteristics.

To effectively model aquifers like the Brunswick aquifer, Anderson et al. [23] recommend employing a dual porosity model approach for solute transport in fractured rock aquifers, as this method improves model predictive capability. The solute transport code MT3DMS [34] incorporates this dual porosity option. Unlike traditional single-porosity models, our approach incorporates mobile (fracture) and immobile (matrix) domains with distinct transport mechanisms (advection vs. diffusion), calibrated using field-measurements. This enables realistic simulation of ‘legacy’ contamination from matrix back-diffusion—a process omitted in prior northern New Jersey studies [13,25]. This modeling approach enables effective simulation of solute transport processes by specifically accounting for solute movement between fractures and the surrounding rock matrix.

1.4. The Typical Hydrogeological Environment of Northern New Jersey

The hydrogeological setting of northern New Jersey is distinguished by its unique interplay of topography, climate, vegetation, geology, and land use. The region’s terrain varies from the gently rolling hills of the Piedmont to the elevated Highlands and the Ridge and Valley provinces [35]. The landscape features a series of peaks and depressions that create a roughly sinusoidal surface, which influences groundwater flow by shaping the water table. Regional topographic slopes, measured from 1:250,000-scale maps, average approximately 0.02.

The area experiences a humid continental climate, marked by cold, snowy winters and warm, humid summers. Average annual precipitation reaches 40 inches, while the frost-free period ranges from 163 to 217 days [36]. The area is predominantly covered by hardwood forests, with species such as Sugar Maple, Red Maple, Beech, Yellow Birch, Sweet Birch, and Red Oak playing a key role in slowing surface runoff and enhancing groundwater recharge [37].

Geologically, the region is largely underlain by the Newark Basin, which includes the Brunswick aquifer. This basin consists of conglomerates, siltstones, and shales, with extensive fractures that store and transmit groundwater. The Brunswick aquifer, formed during the early Mesozoic era [31], is a highly fractured system composed of late Triassic to early Jurassic sedimentary and igneous rocks. It is a productive aquifer, yielding between 100 and 250 gpm [33]. Urban development, particularly in the northeastern part of northern New Jersey, has led to environmental concerns, including contamination from deicing salts applied to roads in winter [20].

Due to all of these characteristics, the Montclair State University community and adjoining areas are identified as a “typical hydrogeological environment” of northern New Jersey. The persistent influx of large crowds (students, faculty, and visitors), combined with heavy daily commuting patterns and related activities, imposes significant environmental stress on the local hydrogeological system.

2. Conceptual Model Design and Calibration

This study utilizes groundwater flow modeling techniques to analyze the flow conditions within the study area. By integrating field observations and well data, a three-dimensional numerical groundwater flow model was developed to simulate the groundwater flow dynamics in the study area. The simulation of steady-state groundwater flow for the study area was conducted using MODFLOW [38], a finite difference code embedded within the Groundwater Modeling System (GMS). For further details on the application and theoretical foundations of MODFLOW, refer to Harbaugh et al. [38]. The model serves as a conceptual framework for understanding the groundwater flow system in the region. MODFLOW requires the integration of surface topography and hydrogeological features, like head boundaries, sinks, and rivers to simulate the groundwater flow. The integration of recharge values, and observed heads enhances model accuracy, yielding an output map that provides valuable insight into the flow directions, velocities, and flux [39].

The geological features of the study area were incorporated into the GMS software and served as the base map for creating a polygon layer to define the model domain. The model domain was constructed around this polygon layer, with rivers bordering the study area assigned as general head boundaries. A topographic map of the region (Figure 2a) was incorporated to provide ground surface elevation data at locations corresponding to the wells used in the study.

Within the model coverage, the study area was divided into zones based on hydraulic conductivity and recharge patterns, which were determined by factors such as land use, geology, and topography. Hydraulic conductivity values were assigned using observed data where available; in cases where such data were unavailable, standard literature values relevant to the area’s geology were used. Recharge values were estimated as fractions of precipitation in the study area and adjusted according to variations in land use and land cover. The model layers were structured from top to bottom, beginning with an unconsolidated water-bearing formation at the surface and ending with a no-flow boundary at the bedrock layer below. Given that rock permeability typically decreases with depth, the model assumes a corresponding decrease in hydraulic conductivity across the layers, with values reducing by an order of magnitude as depth increases.

The model’s dual-porosity framework accounts for fracture-matrix mass transfer, where contaminants diffuse into low-permeability rock matrices—a process critical to predicting long-term contaminant persistence. This approach captures the heterogeneous permeability observed in the Brunswick aquifer, where fractures act as preferential flow paths and matrices as diffusive sinks. The western portion of the study area is partially defined by fracture networks filled with hydraulically conductive materials.

Water level data were incorporated into coverage layers to serve as observational tools. The river network within the study area, including the Third River, was defined in the river coverage, with river stages and riverbed conductance specified. Hydraulic heads in the study area were calculated using observed water levels, elevation data, recharge values, and hydraulic conductivity values. The top elevation values were assigned as the initial heads for the top layer. The hydrogeologic characteristics of the Brunswick aquifer have been well-established through multiple studies [5,30,31,32]. Building upon these existing studies, the groundwater flow model was conceptualized with ten layers to adequately represent the aquifer system. The observation wells incorporated in the model are completed across multiple geologic units, including claystone, shale, siltstone, and sandstone.

Once the groundwater system was sufficiently represented within the modeling coverages, a grid framework was established to facilitate the groundwater model. The MODFLOW model was configured with a grid dimension of 12,500 meters by 17,500 meters and set up as a steady-state simulation. The study area was divided into zones for hydraulic conductivity and recharge, guided by land use, geology, and topography. Groundwater flow across these boundaries was assumed to be negligible. Stream networks were incorporated in the river coverage, with stages assigned according to local hydrologic patterns.

The calibrated groundwater flow simulation model reveals that groundwater heads in the study area vary from approximately 140 meters in elevated regions, such as the First Watchung Mountain, to below 5 meters near the Passaic River to the east (Figure 2b). Calibration of the model involved adjusting hydraulic parameters to closely align simulated hydraulic head values with observed data. Initial recharge values were refined using local annual precipitation data [36] to better capture spatial variability. Recharge values, typically ranging from 5% to 25% of annual precipitation, depend on soil type and land cover [40]. In some zones, initial recharge values were reduced by an order of magnitude to achieve a more spatially consistent distribution across the study area. Initial hydraulic conductivity values were based on hydrogeological properties reported in literature [23,25] and were manually adjusted during calibration. The final calibrated hydraulic conductivity values align closely with field-measured values for the Brunswick aquifer [41] and fall within the expected range for similar rock types as documented in literature [23,25], ranging from 5.59x10^-6 - 1.96x10^-3 meters per second.

Observational wells [42] within the study area were used for model calibration. Parameters such as hydraulic conductivity and recharge were manually adjusted to minimize discrepancies between simulated and observed head values. Figure 2 b shows the output of the groundwater flow simulation, demonstrating that all monitoring wells were adequately simulated to match the field observed values. The calibrated groundwater flow model shows that a significant portion of the groundwater flow in the study area originates from or passes through the First Watchung Mountain area. This area serves as the primary recharge area for much of the study area, contributing groundwater to the two municipal wells located within the study area.

The model indicates that the highest hydraulic heads are located in the southwest part of the study area, while the eastern region exhibits a lower hydraulic head pattern. The analysis reveals a distinct flow regime, consisting of both local and regional groundwater systems. In the local flow system, groundwater diverges westward to the Peckman River and eastward to the Third River from the recharge area near the First Watchung Mountain. The regional flow system, originating from the same recharge area, extends further to the distant Passaic River, located at the far eastern end of the study area. Details on the regional groundwater flow and advective transport on a typical northern New Jersey aquifer are fully documented by Oyen and Ophori [43].

3. Chloride Contaminant Transport Modeling

Following the establishment of a MODFLOW solution, contaminant transport codes are employed for this study, MT3DMS is used with the established MODFLOW solution. This study uses MT3DMS to investigate chloride transport and predict solute behavior during varied travel times. The MT3DMS code is a finite difference upstream solution that effectively captures processes like degradation, dispersion, and diffusion. MT3DMS [34] is a three-dimensional multispecies transport model, simulating advection, dispersion, and chemical reactions of contaminants in groundwater systems [44]. For groundwater modeling that involves significant dispersion and/or chemical reactions, the basic advection – dispersion transport model may be inadequate for establishing the transport processes in highly heterogeneous or fractured groundwater systems. To improve the model accuracy, a dual porosity approach is employed, incorporating a term that describes solute exchange between fractures and the surrounding matrix [23]. Zheng and Wang [34] highlight the advantages of adopting a dual porosity approach when modeling solute transport in complex media. By distinguishing between primary and secondary porosity, this approach accounts for a distinct transport mechanism, where one domain characterized by mobile water and high hydraulic conductivity facilitates advection-dominated transport, while the other domain (secondary porous media), marked by immobile water and low hydraulic conductivity, relies on diffusion-driven transport. Heterogeneity expressed in groundwater systems significantly influences the localized accumulation of constituents within specific geological horizons [31]. The presence of contaminants within fractures generates a pronounced concentration gradient between the fracture fluid and the adjacent intact rock matrix. This gradient drives the diffusive transport of contaminants from the fracture into the surrounding porous matrix, effectively sequestering a portion of the contaminant mass from the fracture's advective flow field. As the contamination plume migrates deeper into the matrix over time, its extent and intensity will be influenced by the processes of diffusion and advection. Ultimately, the removal/cessation of the contamination source will trigger a gradual back-diffusion of the stored contaminants into the fracture openings, the contaminants slowly move back into the fracture, facilitated by fresh water flushing through the fracture network [45]. This phenomenon observed in fractured rocks can cause a recurring challenge during remediation of organic compounds.

According to several studies, chlorides occur naturally, but increased anthropogenic activities have increased their availability in the environment, contributing to the salinization of freshwater resources [2,6]. Dissolved chlorides exhibit conservative properties that actively deter it from full term degradation and sorption processes [46], making it a suitable tracer in hydrologic investigations. The conservative attributes which make chloride to be very soluble in solutions may hold back chlorides in the immobile solution in pores of aquifer materials, while the solute transport continues in the mobile media of the fracture. Therefore, variability in aquifer porosity in a fractured rock aquifer system indicates varied impacts on solute path, necessitating the need for varied input of porosity values during the groundwater modeling process. Chloride concentrations exceeding background levels originate from diverse sources; deicing salt application and degradation of organic compounds (such as TCE) are potential sources of elevated chloride concentrations. Chlorides from anthropogenic sources can serve as tracer and/ or surrogate for organic compounds and contaminants in general [47]. While the US EPA [9] and New Jersey Groundwater Quality Standards [12] establish 250 mg/L as the maximum recommended chloride concentration for drinking water, multiple studies [2,7,20] have documented widespread exceedances of this threshold in groundwater systems. For this study, a chloride ion concentration of 250mg/l and lower is indicated as a threshold for available freshwater resources. Therefore, Chloride ion concentrations above 250mg/l are deemed to be below the groundwater quality standard.

The transport and fate of chloride during varied travel times were studied to understand the solute transport dynamics in the Brunswick aquifer system. The study created a scenario which involves the introduction of chloride salt mass of 1000 kg into the MSU university community, to observe the varied output concentrations and the plume behaviors over varied travel times. For this study, only hypothetical chloride concentrations introduced at MSU community are considered, chloride concentrations from rainfall or other processes such as runoffs are not within the scope of this study. The concentration levels and travel times of chloride were observed at important landmark resources, including, the nearby water reservoirs (Cedar Grove and Great Notch reservoirs), surface-water bodies (Third River, Peckman River, and Passaic River) and the municipal supply wells for various travel times. The chloride concentration level of 250 mg/L is indicated as a reference point during the travel times through the various landmarks. To determine the chloride concentration introduced/infiltrated into the groundwater, chloride salt mass (1000 kg) is converted to mass per volume (per Mihelcic & Zimmerman [48]), the volume is determined from an area coverage of MSU with a supposed snowfall of 2 inches. The MSU area covers about 900,000m², and the study estimates that roughly one-tenth of the chloride concentration infiltrates into the groundwater. For a salt mass of 1,000 kg applied over this area with a snowfall of 2 inches (0.0508 meters), the chloride concentration is approximately 21,870 mg/L. The transport analysis assumes that one-tenth of the concentration (that is 2,187 mg/L) from the applied salt infiltrates the groundwater.

Following the MODFLOW transport analysis, few parameters are required to build on the MODFLOW model to run the MT3DMS transport model. With the aid and setup of the advection, dispersion, and chemical reaction packages, the MT3DMS model simulation was run to provide various concentration gradients and varied travel times. The dispersivity parameter, horizontal transverse dispersivity was assumed to be orders of magnitude lower than the longitudinal dispersivity. The value of horizontal transverse dispersivities was assumed to be higher than that of vertical transverse dispersivities by a magnitude of one order [40,49]. The longitudinal, horizontal and vertical transverse dispersivities of 20, 0.1, and 0.01 m were used, respectively. Dual domain mass transfer (without sorption) was used for the transport model, with mobile and immobile porosity of 0.2 and 0.4 respectively (porosity of aquifer materials with similar properties as the Brunswick aquifer as indicated by Anderson et al. [23] guided the selection of porosity values) and mass transfer coefficient of 0.0005 per second. Chloride concentrations (as stresses) at MSU were set at the start of the simulation and remain constant throughout the simulation. The stress period length (simulation period) for this study is 450 years; in order to capture the dynamic fluctuations of chloride concentrations, the stress period is partitioned into a series of time steps.

For this study, contaminant transport modeling is done for salt constituents (chlorides) acting as surrogates for contaminants found in groundwater. The established model is applied for travel distances and chloride concentrations at important landmark resources during varied travel times. Figure 3a,b show the contaminant transport from the time of introduction into MSU community to 450 years after, a complete simulation period at which the contaminant traveled through the study area.

Figure 4 shows the arrival of chloride at the Third River. The model indicated that after about 5 years of travel, chloride contaminants in the local discharge area of the Third River increase from 0mg/L before the introduction of the contaminants at study site to more than the recommended level of 250mg/L. At this stage, the chloride contaminant has traveled a distance of 450m northeast along local flow paths to reach the Third River and accrued that much concentration during the 5 years period. The Third River is a headwater stream that has possible interactions with groundwater (including groundwater recharge) along the path of the river and the arrival of chloride in the river may possibly increase the contaminant’s concentration levels in certain areas via processes beyond the scope of this study.

Figure 5 illustrates the arrival of chloride at the Cedar Grove reservoir. The model indicates that after approximately 15 years of travel, chloride concentrations in the area around the reservoir increased from 0 mg/L (prior to contamination) to levels exceeding the recommended threshold of 250 mg/L. During this period, the chloride contaminants traveled a distance of 600 meters west within local flow systems to reach the Cedar Grove reservoir, accumulating significant concentrations over the 15-year timeframe.

Figure 6 depicts the arrival of chloride at the Great Notch reservoir. According to the model, after approximately 32 years of travel, chloride concentrations in the reservoir increased from an initial level of 0 mg/L to levels surpassing the recommended limit of 250 mg/L. During this period, the chloride contaminants traveled 900 meters north along local flow paths to reach the Great Notch reservoir, accumulating significant concentrations over the more than three-decade span. By this stage, the contaminant plume had traveled approximately 1,700 meters in a northwest direction from MSU. The Great Notch reservoir, owned and managed by the Passaic Valley Water Commission, plays a critical role in meeting the water needs of Passaic County communities. As an open reservoir, its water quality can be influenced by various sources and processes beyond those discussed in this study.

Figure 7 depicts the arrival of chloride at the Peckman River. According to the model, after about 50 years, chloride levels in the river rose from 0 mg/L (pre-contamination levels) to concentrations surpassing the recommended limit of 250 mg/L. During this time, the chloride contaminants traveled 2,200 meters northwest within a local flow system from the regional recharge area (MSU) to reach local discharge points (Peckman River), building up substantial concentrations over the five-decade period. The Peckman River flows northward, eventually joining the Passaic River to the north of the study area.

Figure 8 indicates the arrival of chloride on the Passaic River from the MSU community. The model shows that after about 79 years of travel time, chloride contaminants in the Passaic River north of MSU would have concentration levels beyond the set recommended level of 250mg/l, from an initial chloride concentration value of 0mg/l. At this stage, the chlorides traveled about 2700m northward within a regional flow system to arrival at regional discharge points (Passaic River) and accumulated chloride concentrations surpassing recommended levels during the 79-year period. During this time, the contaminant plume had traveled about 3000m NNW and 3500m NE direction away from MSU along local and regional flow paths. The Passaic River is a regional discharge area for most of the region and possibly accumulates contaminants of varied sources via processes such as groundwater transport, runoffs, among others.

After the 450-year simulation period, landmarks like the Loraine well, and Rand well, did not show chloride concentrations beyond the recommended limits, as these landmarks were not in the groundwater transport pathway of the chloride contaminants introduced at the MSU community (Figure 9). Therefore, there are limited chances that possible contamination of the Loraine well, and Rand well are from groundwater processes originating from MSU area of the First Watchung Mountain. At this stage, the chlorides traveled approximately 7,000 meters within a regional flow system, moving northeast from the First Watchung Mountain area (MSU) toward the Passaic River and accumulated chloride concentrations of about 150 mg/L, which remained below the recommended threshold.

3.1. Implications of Chloride Transport for Groundwater Management

This hypothetical contaminant transport analysis, using a dual-porosity approach to represent the fractured Brunswick aquifer, reveals significant implications for groundwater quality and management, particularly concerning the transport of deicing salts (chlorides) and their associated impacts. The findings demonstrate that chloride contaminants that may be introduced at the MSU community from deicing salts can travel through the groundwater system and reach critical water resources, such as the Third River, Cedar Grove reservoir, Great Notch reservoir, Peckman River, and Passaic River, within decades. These contaminants (e.g. chlorides) would often exceed EPA-recommended levels (250 mg/L), posing risks to water quality, ecosystems, and public health. The potentially rapid accumulation of chlorides in surface water bodies, such as the Third River within 5 years and the Cedar Grove reservoir within 15 years, highlights the urgency of addressing anthropogenic sources of salinity, including deicing salts, fertilizers, and sewage overflows. The local groundwater flow systems appear to be more susceptible to contamination, and the travel times indicated in the model may be shorter, as water in discharge areas (such as the Third River) within the study region could originate from entirely different sources. As a result, significant variations in chemical quality are likely to occur [50].

Increasing salinity in groundwater due to anthropogenic activities can lead to several challenges. One significant issue is the "salting-out effect," where elevated salinity enhances the retention of per- and polyfluoroalkyl substances (PFAS) and other contaminants in groundwater [22]. This phenomenon complicates remediation efforts, as contaminants trapped in the pores of aquifer materials can slowly diffuse back into the water, re-contaminating remediated wells [16,17], especially in the Brunswick kind of fractured aquifer. Additionally, salinity is known to mobilize certain contaminants [14,15], further exacerbating groundwater quality issues. These processes demonstrate the need for comprehensive strategies to mitigate rising salinity and its cascading impacts on groundwater remediation. These impacts are quantified through the transport modeling results presented in this study.

The transport modeling results demonstrate that chloride contaminants migrate through the fractured aquifer via two distinct flow systems: (1) rapid transport through local flow paths (reaching discharge points like the Third River within 5 years) and (2) slower migration through regional flow systems (reaching the Passaic River after 79 years). This disparity reflects fundamental hydrogeologic controls—fracture networks enable fast contaminant transport in local systems, while longer flow paths and matrix diffusion processes dominate regional transport.

Chloride concentrations exceeded the 250 mg/L threshold 15 times faster in local streams compared to the regional base-level river, highlighting the aquifer’s dual-domain transport behavior. These findings align with observations of 'salting-out' effects, where matrix diffusion in regional systems prolongs contaminant retention, complicating remediation efforts in fractured aquifers [16,22].

The findings also emphasize the critical need to protect regional recharge areas (such as those around the First Watchung Mountain, where contaminants migrate most efficiently) and vital water resources like the Great Notch reservoir, which serves as a crucial water supply for Passaic County communities. The model shows that chloride concentrations in the reservoir exceeded recommended levels within 32 years, highlighting the vulnerability of such resources to contamination. This stresses the need for proactive measures in urbanized fractured aquifers, where reducing salt inputs and managing stormwater at recharge points can mitigate long-term salinity buildup. Key strategies include targeted deicing salt reduction, improved stormwater management, and land-use regulations to minimize contaminant infiltration into groundwater.

Furthermore, the study reveals that while some wells, like the Loraine and Rand wells, were not directly affected by chloride transport from the MSU community, the broader regional flow system still poses risks. The transport of chlorides over long distances (e.g., 7,000 meters toward the Passaic River) demonstrates the potential for widespread contamination, necessitating regional-scale monitoring and management efforts. While our model focuses on northern New Jersey, the framework applies generally to fractured aquifers under anthropogenic salinization stress. The rapid (<5-year) contamination of local streams highlights the urgency of regulating winter salt use in recharge zones, whereas delayed regional transport warns of long-term ‘legacy’ salinity. These insights are critical for cities relying on fractured bedrock aquifers, from northern New Jersey to water-stressed regions worldwide. Addressing these challenges may require a multidisciplinary approach, combining groundwater modeling, contaminant tracking, and targeted remediation strategies to ensure the protection of water resources,

These transport mechanisms are not unique to the Brunswick aquifer. Beyond northern New Jersey, our findings have immediate relevance for urbanized regions reliant on fractured bedrock aquifers, where anthropogenic salinization threatens water security. By demonstrating how fracture-dominated transport interacts with matrix diffusion, this study provides a framework for policymakers to prioritize recharge-zone protection while addressing long-term contaminant retention in fractured systems. The 'salting-out' effects observed in this study suggest similar retention risks for co-migrating contaminants in other fractured systems.

4. Conclusion

In this study, the hypothetical introduction of chlorides into the aquifer at the MSU community, an area with characteristic representation of northern New Jersey, was analyzed using a solute dual-domain transport model. The model revealed varying chloride concentrations, travel distances, and travel times in the regional groundwater flow pattern of a typical hydrogeological environment of northern New Jersey. The transport modeling identifies a defining characteristic of the Brunswick aquifer: a dichotomy between rapid fracture-dominated chloride transport to local discharge zones (threshold exceedance in about 5 years) and matrix-retarded migration to regional discharge points (79 years). This behavioral contrast exemplifies how fracture connectivity governs short-term contamination risks, while matrix storage capacity determines long-term plume persistence—a critical consideration for urbanized fractured aquifer management. Over a 450-year simulation period, chlorides traveled roughly 7,000 meters within a regional flow system, moving northeast from regional recharge area of the First Watchung Mountain toward the regional discharge area (Passaic River ). While groundwater transport processes may take decades for contaminants in northern New Jersey to move from the recharge area to the regional discharge area of the flow system, surface water pathways could significantly accelerate this process in northern New Jersey.

Given that the Third River and Peckman River are tributaries of the Passaic River, chloride-contaminated water exceeding 250 mg/L could reach the Passaic River within 5 years via surface water, not accounting for runoff processes that could further reduce this timeframe. This study scenario describes the complexity of contaminant movement through typical groundwater flow systems of northern New Jersey. Chloride serves as an effective tracer and surrogate for various groundwater contaminants due to its conservative properties. Although other contaminants (or even chlorides) may not follow the exact timelines outlined in this study due to factors like the widespread application of deicing salts and runoff processes, there is a high likelihood that contaminants will follow the modeled groundwater transport pathways. It is important to note that this study focuses primarily on groundwater transport processes, with deicing salts identified as the primary source of chloride, which also participates in other transport mechanisms.

The Brunswick aquifer, characteristic of northern New Jersey’s typical hydrogeological environment, exhibits unique geological features, particularly its fracture systems, that complicate contaminant remediation due to salting-out effects. In such aquifers, contaminants may persist longer, with potential recurrence even after source cessation or remediation efforts. As indicated by Hort et al. [17], contaminant loading into fractured aquifer systems leads to diffusion into the matrix, slowing the contaminant plume. After remediation and the introduction of clean water, adsorbed contaminants in rock pores can diffuse back into the groundwater, complicating cleanup efforts.

This study provides a strong foundation for future water resource investigations in the region. A solute-transport model could be developed to analyze the transport and fate of contaminants under both pump-and-treat and non-pumping remediation methods. Such a model could quantify reductions in contaminant concentrations due to dispersion, dilution, degradation, and diffusion along flow paths originating from areas prone to high contaminant accumulation. Addressing contaminant transport challenges will likely require a multidisciplinary approach, integrating groundwater modeling, contaminant tracking, and targeted remediation strategies to protect water resources and public health. Generally, the findings of this study highlights the urgent need for integrated water management practices to mitigate the impacts of increasing salinity and ensure the long-term sustainability of groundwater systems.