δ34S results versus δ18O for all analyzed samples are shown in Figure 8a where they can be compared to typical domains of known sulphate sources reported by Mayer [48]. Results are also shown with regards to sulphate concentrations in Figure 8b (where the point size represents the relative sulphate concentration). The isotopic signature of leachate sulphate from the backfill eluates is close to that observed in groundwater. Compared to Mayer's sulphate sources [48], groundwater samples match with the 'soil SO4 ' zone, which includes organic (ester and CS-mineralization) and inorganic sulphides. The eluate samples from the backfill are within the anthropogenic source zone of atmospheric deposition near the boundary of the 'soil SO4 ' zone.

As described previously, the studied backfill soil may be highly heterogeneous: a mixture of sandy loamy soil with recycled construction and mining wastes that contain shale and coal in small quantities. The isotopic signature of sulphate in this backfill therefore represents the signature of different anthropogenic sources.

In an urban context, Bottrell et al. [17] studied, using isotopes, sulphate sources around the city of Birmingham. Sulphate isotopic signature, for a group of wells in the city center and some industrial sites, was quite similar to our study results with δ34S between 0.2‰ and 3.2‰ and δ18O between 2.8‰

and 7.73‰. This isotopic signature was interpreted as from SO4 2− urban pollution resulting from artificial (made ground) soils and/or from wastewater.

**Figure 8.** δ34S versus δ18O values of sulphate from groundwater and eluates samples compared to typical values depending on: (**a**) Sulphate sources and (**b**) with regards to SO4 2− concentrations. The isotopic composition for the sulphate sources are taken from Mayer [48]. Concentration of SO4 2− in groundwater versus: (**c**) δ34S and (**d**) δ18O.

The isotopic signature obtained in our study reflects a mixture between background sulphate in the aquifer and additional sulphate released from the backfill layer. In the study by Jurado et al. [49], the isotopic signature of SO4 2− in the same aquifer for points with low sulphate concentrations (between 34.54 and 71.49 mg/L) was between −2.08‰ and −0.13‰ for δ34S, and between 2.1‰ and 4.3‰ for δ18O. These results are similar to our sampling points with low SO4 2− concentrations, and therefore we can consider these values as the natural background signature of sulphate in the chalk aquifer.

Locally, in the study area, the higher the sulphate concentrations increase, the more we have an enrichment of SO4 2− in groundwater at δ34S and δ18O, and a trend towards the isotopic signature of backfill eluates Figure 8c,d. Sampling wells showing the highest SO4 2− concentrations are located below the backfill of the industrial site.

The increase in sulphate concentration is clearly accompanied by an increase of δ34S-SO4. For δ18O-SO4, this trend is not so evident because of the background isotopic signature of sulphate in groundwater comes mainly from the mineralization of organic matter which causes a depletion in δ18O-SO4 but not of δ34S-SO4 which does not show significant fractionation through this process [50].

#### **6. Synthesis of Process Leading to Groundwater Mineralization Changes**

The main processes explaining groundwater mineralization changes that occurs at the industrial site are conceptually summarized in Figure 9:

**Figure 9.** Conceptual scheme summarizing the obtained results with the main process controlling the groundwater mineralization under the industrial site (Scale not respected).

The investigated backfill layer had been in place for more than 40 years prior the present study. It is made of a mixture of materials with considerable calcium sulphate release capacity. The detected traces of oxidized pyrite with hematite indicate that the backfill soil contained pyrite in the past. Gypsum was also present, either with recycled building materials or formed as a result of calcite reaction to acidity. Some sulphate in the backfill soil were in amorphous form. In addition, the improper handling of products containing CAHs has produced a contamination of soil and groundwater. The mixture of pollutants had enough time to percolate across the unsaturated zone with the infiltrated water to reach the saturated zone.

The current groundwater hydrochemistry changes have resulted from the combination of an increase of calcite dissolution along with the migration of calcium and sulphate from the backfill soil, as verified by the correlation of sulphate concentrations with non-carbonated calcium (Ca-HCO3). Furthermore, isotopic analyses results confirm the mixing between the sulphate released from the backfill soil with the background sulphate in the aquifer.

Considering the leaching test results for the first 4 fractions, where easily leachable components are leached quite extensively, calcium and sulphate are the dominant compounds in the backfill eluates. Results extrapolation at the field scale shows that for a period between 0.86 to 8.59 years, the infiltrated water passing through the backfill layer produces leachate with sulphate concentrations between 9.2 and 3.3 times the background concentrations of sulphate in the aquifer. While for calcium, for the same duration, concentrations in infiltrated water correspond to values between 1.9 and 0.7 times the background concentration of calcium in the aquifer.

On the other hand, geochemical simulations provided an assessment of calcite dissolution in groundwater in the presence 1,1,1-TCA degradation reactions. Results showed that for the maximum concentration of 1,1,1-TCA (1100 μg/L) observed in groundwater, HY/DH degradation reactions decrease the pH by a value of 1.01 units for a 1-year degradation time. Consequently, the calcite dissolution produces an additional amount of 0.64 times the background calcium concentration, and 0.93 times the background bicarbonates concentration.

#### **7. Conclusions**

Through this study, the understanding of a pollution problem in the considered chalky aquifer was improved. Using a combination of tests on the backfill material extracted from the site and groundwater quality analyses, the mechanisms that lead to changes in groundwater mineralization at this site was clarified. These changes are mainly due to two sources. The first is the increase in calcite dissolution as a buffer reaction to the acidity released by the degradation of 1,1,1-TCA by HY/DH. Geochemical simulation results showed that calcite can release up to 95.28 mg/L of calcium and 289.93 mg/L of bicarbonate during one year of 1,1,1-TCA degradation.

The backfill layer at the industrial site represents a second source that affected the hydrochemistry by releasing calcium sulphate that migrated to the saturated zone with recharge water. Leaching test results showed that for a period between 0.86 and 8.59 years, the average concentration in backfill leachate is 302.52 mg/L for sulphate and 104.39 mg/L calcium.

Thus, in the present case study, the improper backfill materials quality caused an additional pollution source influencing groundwater mineralization. A better control of backfilling materials would have prevented a part of the groundwater contamination. Furthermore, the advanced characterization of backfill soil and groundwater hydrochemistry provided an appropriate estimation of the effect of 1,1,1-TCA degradation reactions on the hydrochemistry compared to backfill leachates effect, leading to a more correct conceptual site model that will improve the remediation plan.

The current research demonstrated also the benefits of analyzing major inorganic chemical elements in cases of local pollution with organic pollutants such as CAHs. These analyses not only make it possible to identify other potential pollutions, but also to improve the understanding of CAHs degradation reactions in groundwater.

**Author Contributions:** Conceptualization, Y.B., S.B. and A.D.; methodology, Y.B., S.B., P.O. and A.D.; investigation, Y.B. and D.G.; resources, P.J., D.G. and A.D.; writing—original draft preparation, Y.B.; writing—review and editing, S.B. and A.D.; supervision, S.B. and A.D.

**Funding:** This research was partly funded by the research unit Hydrogeology and Environmental Geology at the University of Liège.

**Acknowledgments:** Many thanks to Joël Otten for XRD analysis on soil samples and for the analysis of inorganic chemical elements performed on water samples at the University of Liège. The authors are grateful to Dr. Kay Knöller from UFZ for the isotopic analysis. Many thanks to Dr. Mohammed Tayeb Sadani and to Youcef Hakimi for their enriching discussions.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

#### *Article*

# **Groundwater Contribution to Sewer Network Baseflow in an Urban Catchment-Case Study of Pin Sec Catchment, Nantes, France**

### **Fabrice Rodriguez 1,2,\*, Amélie-Laure Le Delliou 1,3, Hervé Andrieu 1,2 and Jorge Gironás 4,5,6,7**

Received: 21 December 2019; Accepted: 21 February 2020; Published: 3 March 2020

**Abstract:** Sewer systems affect urban soil characteristics and subsoil water flow. The direct connection observed between baseflow in sewer systems under drainage infiltrations and piezometric levels influences the hydrological behavior of urban catchments, and must consequently be considered in the hydrologic modeling of urban areas. This research studies the groundwater contribution to sewer networks by first characterizing the phenomenon using experimental data recorded on a small urban catchment in Nantes (France). Then, the model MODFLOW was used to simulate the infiltration of groundwater into a sewer network and model dry weather flows at an urban catchment scale. This application of MODFLOW requires representing, in a simplified way, the interactions between the soil and the sewer trench, which acts as a drain. Observed average groundwater levels were satisfactorily simulated by the model while the baseflow dynamics is well reproduced. Nonetheless, soil parameters resulted to be very sensitive, and achieving good results for joint groundwater levels and baseflow was not possible.

**Keywords:** groundwater; urban hydrology; drainage; modeling; sewer; baseflow

#### **1. Introduction**

Urbanization modifies land use and affects soil, sub-soil and subsurface processes in different ways. Urban features as well as surface and underground infrastructure can have a strong impact on groundwater levels. Some causes explaining a reduction in groundwater include an infiltration decrease due to additional imperviousness, groundwater pumping for various urban water uses, and groundwater flow into drainage trenches [1]. On the other hand, leakage from water supply and waste water networks becomes a source of recharge for urban groundwater [2–5]. Furthermore, leakage from waste water systems is also a possible source of groundwater contamination [6]. Despite aquifer levels and stream flow having proved to be related in various rural contexts [7,8], few studies have focused on examining this relationship in urban catchments [5,9–11]. Soil water infiltration in sewers, however, is a phenomenon affecting urban hydrology, sewerage and waste water treatment plant management [12,13]. Initially, the presence of soil water in wastewater sewerage was revealed

by [14,15], when flow variations in separate wastewater sewerage during wet weather periods were commonly attributed to inappropriate or irregular stormwater runoff connections. These authors showed that groundwater seep into the sewers through defective cracks once the water table level reaches the depth at which sewers are buried. This phenomenon has been recently reinforced by [13] for several Flemish catchments; some interesting stable isotopes methods have been used to detect and estimate groundwater intrusions in sewers in Nancy, France [16] and in Brussels, Belgium [17] and temperature-based methods in Trondheim, Norway [18]. Similarly, [9] showed that soil water may explain the variations of runoff coefficients observed on a small urban catchment in the city of Nantes, France, as groundwater drainage flow becomes significant once the water table exceeds a given threshold. The flowrate in sewers is the sum of runoff during rain events, wastewaters from housing (for separative wastewater or combined sewers), and a groundwater infiltration component, representing the baseflow in sewers; this component is often part of "extraneous" water in sewers. Interestingly, [19] suggest the appropriate terminology of urban "karst" to denote the multiple water soil and pipe interactions below the surface. Despite the relevance of improving our understanding of the impact of groundwater drained by urban sewer systems, only a few efforts are currently devoted to the problem. Groundwater drainage by sewers is often noted in literature [20], although is rarely quantified when assessing the urban water budget. Only a few papers deal with the determination of the urban water budget, likely due to the lack of reliable data of all the water budget components (i.e., rainfall, evapotranspiration, surface and groundwater flow rates). This lack of knowledge must be addressed as urban stormwater practices increasingly consider innovative infiltration-based technologies and approaches to mitigate the hydrological impacts of urbanization [21]. Thus, a better understanding and quantification of the fate of urban soil water becomes essential [22,23].

The first attempts to represent the interaction between surface water and the aquifer in urban integrated hydrological models were simple and conceptual. Aquacycle [24] and the Storm Water Management Model (SWMM, [25]) are examples of such models. In a comprehensive review of ten stormwater models, [26] identified the models able to simulate the groundwater baseflow component, mostly through a conceptual linear reservoir. Some of the more recent physically-based modeling efforts consider soil drainage by sewers. An integrated sewer-aquifer model was developed [27] to test the effects of future buried drains on both groundwater flow and sewer infiltration phenomena. This integration proved to be relevant to describe the groundwater-sewer interactions. The Network Exfiltration and Infiltration Model (NEIMO) model [28] was developed in close connection with already available groundwater models such as MODFLOW [29] or FEFLOW [30], to simulate infiltration and exfiltration processes. This modeling approach requires developing an integrated modeling framework like Urban Volume and Quality (UVQ) developed by CSIRO, Australia [31,32]. MODFLOW was used [33] to analyze the relevant infiltration parameters at a local scale, and compared a 1D-infiltration approach and a MODFLOW modeling approach to model groundwater infiltration at a larger scale. Overall, coupling groundwater models with the simulation of infiltration/exfiltration processes is a major challenge, especially due to the high uncertainty associated with these models [34]. The coupling of a hydraulic model was realized by [35] with both a groundwater model and a sewer failure estimation approach to identify the sections potentially affected by infiltration in coastal urban areas. Satisfying results were obtained by [11] with regard to the groundwater infiltration by coupling two commercial models MIKE URBAN for simulation of sewer flow and MIKE SHE for simulation of groundwater transport. The Urban Runoff Branching Structure Model (URBS-MO) [36] was developed to represent surface and subsurface water flows while focusing on the impact of sewer networks. Although the model simulates groundwater drainage by sewer, real observations at a local scale to test this component of the model were not available. Finally, [5] coupled the hydrological/water management WEAP model and MODFLOW to verify a strong stream-aquifer interaction in areas with shallow groundwater, as well as quantify the local recharge associated with pipe leaks and inefficient urban irrigation. Overall, the literature review shows that much more experimental data and modeling using comprehensive 3D modeling tools such as MODFLOW, are needed to better

understand and quantify the impact of groundwater infiltration and groundwater levels on low flows in urban environments.

This study investigates the interactions between urban groundwater and baseflow in wastewater and stormwater sewer systems, and its first goal is to better understand their role on the hydrological budget of an urban catchment. These interactions are analyzed and quantified using a field study conducted in the "Pin Sec" catchment, located in the city of Nantes (France) [37]. The second objective consists of testing the ability of a groundwater modeling tool to mimic these interactions. The results are used to build a numerical model focused on simulating the hydrological behavior of the urban soil and low flows in the sewer system. The outline of this paper is as follows: the next section presents the case study and the methods adopted in this work, which focus on using MODFLOW to model baseflow and interactions between the sewer systems and groundwater at catchment scale. The results section presents and compares observed and simulated groundwater levels and baseflow discharges. Finally, the last section summarizes our results and highlights the main conclusions.

#### **2. Materials and Methods**

#### *2.1. In-Situ Interactions Between Baseflow and Groundwater Levels*

#### 2.1.1. Case Study

The experimental area is part of the Observatoire Nantais des Environnements Urbains (ONEVU) initiative, devoted to long-term monitoring of water and pollutant fluxes and soil-atmosphere exchanges in urban settings [37]. The Pin Sec catchment is located in the East of Nantes (France) between the Loire and the Erdre rivers. The area is under an oceanic climate, with frequent but not very intense rainfalls. The Pin Sec neighborhood developed between 1930 and 1970 is mainly composed of single and multi-family housing (Figure 1). This 31 ha catchment (called urban catchment) belongs to a larger hydrogeological catchment of 120 ha, both of them are represented on Figure 1. This catchment has a gentle slope, the highest altitude of the catchment is 28 m asl and the lowest is 13 m. The Pin sec catchment has an imperviousness of 45% and is equipped with a 50 years old separate sewer network. The wastewater and stormwater sewer systems have total lengths of 7.3 and 4 km, and mean depths of 2 and 2.9 m respectively. Although the wastewater sewer network is denser than the stormwater system, both are mainly superimposed, and both present cracks or faulty sealing joints. Both system outlets, while collecting water from the same geographic area, are not located in the same place. The stormwater system is connected to the Gohards river, whereas the wastewater sewer system drains into a downstream combined sewer system. Gohards river stream is an old river which was buried between 1945 and 2012 and which has been re-opened through a renovation urban planning of this area in 2012. It is a perennial stream.

Recorded data include (1) 5 min rainfall records from 1999 to 2010, deduced from the average of 3 rain gauges covering the study area, (2) 5-min flow discharge records at the outlet of both the stormwater and the wastewater systems from 1 September 2006 to 31 August 2008, (3) 20-min groundwater levels measured with pressure sensors in eleven piezometers located throughout the catchment from 1 September 2006 to 31 August 2008 (one piezometer, PZGN being only used to map the groundwater level contours), and (4) other meteorological records as well as evapotranspiration estimated with the Penman method, measured 12 km away from the catchment outlet, from 1999 to 2010. We focus our analysis on the main period between 1 September 2006 and 31 August 2008 (i.e., two hydrological years), but we also used year 2002 for some modeling tests as rainfall characteristics were representative of the average conditions in the area.

**Figure 1.** Pin Sec catchment within the Observatoire Nantais des Environnements Urbains (ONEVU) territory and location of the piezometers (green circles). The dotted brown and continuous black lines represent the boundary of the urban catchment and the hydrogeological catchment respectively. The large arrows represent the main sewer flow directions. By convention, north is placed at the top of the following maps. PZPS is located 47◦14 43.0 N 1◦31 09.1 W.

#### 2.1.2. Geology and Hydrogeology

The city of Nantes has developed on the Armorican massif, and especially in the south Armorican shear zone. The geology is mainly composed of a mica-schists bedrock and a stack of alluvial deposits along the Loire River. On the Pin sec catchment, the mica-schists bedrock is covered with altered mica-schists, silty eolian deposits and alluvial deposits along the Gohards stream (Figure 2a) [38,39]. The piezometers are mainly located in the top layers and especially within altered mica-schists, with a drilling depth varying from 4 to 8 m (Figure 2b); for two of them, the depth was limited due to hard non altered micaschists. The hydraulic conductivities estimated from in-situ measurements by water bail tests [40] using the Hvorslev method [41] varies from 8.5 10−8 to 1.810−5 m s−1 [42]. Piezometer PZGN is out of the hydrogeological catchment and was not considered for further results. From now on groundwater level will designate the level of the saturated zone, which is equivalent to the local groundwater level in our case study.

The piezometers were used for mapping the behavior of the water table on the catchment. Groundwater level contours were interpolated using inverse distance weighting (IDW) techniques with Vertical Mapper© in Mapinfo platform, and used to determine groundwater main flow directions. The groundwater levels of the high water table periods are presented in Figure 3. Note that the upstream condition of the water table is given by two piezometers (PZCRI and PZGO) located at the upper limit

of the catchment, and that groundwater moves in the northwest-southeast direction towards the stream outlet following the surface slope. The low water table condition is characterized by a similar profile.

**Figure 2.** (**a**) Main geological layers on the Pin sec catchment: (non-altered micaschists are not visible because below the other layers; (**b**) Soil profiles of the set of monitored piezometers, ordered from upstream to downstream. Colors representing the geological layers on (**a**) refer to the Legend in (**b**).

**Figure 3.** Groundwater level contours (hydraulic heads in m) deduced from the piezometric data interpolation on 7th March, 2007.

#### 2.1.3. Baseflow in Sewer Networks

Both waste- and stormwater systems outlets are equipped with continuous flow rate measurement devices that use triangular weirs. Nonetheless, measuring the lowest flow rates within sewers has proved difficult as dry-weather flows do not necessarily reach the minimum water depth (~3 cm) needed for proper velocity measurement. Diurnal patterns of the wastewater flow rates show a typical domestic double-bump flow during dry weather time (Figure 4), and an increase of this flow during rain weather, due to parasitic waters (due to both wrong sewer connexions and groundwater infiltration). The stormwater flow rate is typical for urban catchments, with a fast hydrological response related to the runoff produced by impervious areas. In this study, baseflow is considered on a daily basis. Stormwater baseflow is basically the minimum flowrate of each day. Because of the small size of the catchment, minimum wastewater flows corresponding to water draining from the water table are assumed to occur in the middle of the night [43–45]. Thus, we assume that the flow rate at 4:00 A.M. is representative of this drainage flow. Those days in which rain takes place are removed from the wastewater and stormwater discharge records to assess baseflow at the outlet of both systems during the study period (i.e., 65% of the time).

**Figure 4.** Dynamics of flowrates in both sewers (wastewater in red and stormwater in blue) during 10 days in May, 2007. The yellow dots show the minimum nighttime wastewater flow, used for the daily baseflow estimation.

#### *2.2. Groundwater Modeling at the Urban Catchment Scale*

Our aim in proposing a model to represent soil-sewer interactions is to improve our understanding of water fluxes in urban catchments. We use the groundwater model MODFLOW [29] to simulate these interactions because of (1) its widely-used extensive modeling capabilities [46], (2) its drainage modeling options to account for soil discontinuities using fine discretizations, and (3) the possibility of easily linking MODFLOW with other models. Few studies, however, have used the model in urban groundwater modeling [5,33,47–49], and in doing so, quite large grid resolutions (100–500 m) are adopted. In our study MODFLOW is used at a much finer resolution to represent sewers or buried trenches and their discontinuities in detail. Simple application conditions have been retained, because the use of MODFLOW in an urban context is not frequent and a cautious approach has been adopted. Moreover, the urban soil is both vertical and horizontal heterogeneous, and therefore not well documented.

#### 2.2.1. Main Modeling Principles and Assumptions

MODFLOW uses a finite-difference approximation and an implicit scheme [29] to solve the 3D groundwater flow diffusion equation. In our implementation, the geological features of the catchment are deduced from the geological characteristics found in the case study (i.e., silt, altered and non-altered mica-schist). Two soil layers are defined from base to top: (1) a basement layer consisting of non-altered mica-schists; (2) a topsoil layer consisting of silty eolian deposits and altered mica-schists on the main part of the area, and alluvial deposits on the downstream part; this topsoil layer contains the sewer trenches.

Three MODFLOW modeling options were tested to represent the sewer trench, as this element cannot be explicitly implemented in the model: (1) a basic modification of the hydrodynamic properties of the cells containing the trench; (2) a field trench parallel to the surface slope with a small drainage capacity, using the drain package (DRAIN); and (3) a representation of the sewer network using a river with the river package (RIV).

A preliminary sensitivity analysis was carried out at a street scale for a 100 × 100 m unit element, both in steady state and transient conditions (not presented here). This analysis helped to understand the impact of the modeling configurations (grid resolution and the sewer trench representation option) on both groundwater levels and flows simulation. Moreover, the influence of the main parameters used in the model (i.e., sewer trench conductance, hydraulic conductivities, specific yield and specific storage) was analyzed. The choices about the parameters or configurations described below were not adopted to replicate experimental results, but to retain realistic characteristics of the subsurface processes to be represented.

The grid resolution in our analysis varied between 1 and 10 m and affected mainly the water flows. Finally, a 4 m grid cell resolution was chosen. This size is quite large to represent the sewer trenches encompassing sewer pipes, but is a good compromise given the area of the study catchment (120 ha). Along the three tested options, the option using the DRAIN package to represent the sewer trench proved best to represent the saturated zone drainage by a sewer, because it can better reproduce the decreasing groundwater near the trench. Specific yields and conductivities are the main sensitive parameters, influencing both the groundwater levels and the flows. Due to the low impact of anisotropy in our sensitivity analysis, conductivities are considered homogenous in the three directions. The water flux increases with the hydraulic conductivities, while the groundwater level decreases. The initial hydraulic conductivities are estimated from in-situ measurements described above (§2.1.2). The influence of the specific yield of the topsoil layer is significant for the sewer trench outflow: the downstream flow reduces by 7% when specific yield decreases from 20% to 5%, which are typical values assessed by [50].

The field drainage system representing the sewer trench is characterized by a drain hydraulic conductance (L2/T) (conductance thereafter), which depends a priori on the material and characteristics of the drain. Since the conductance is usually unknown, it must be estimated by calibration [49,51]. We first assessed the range of variation of the conductance from the sewer baseflow observations. As discussed by [51], the discharge rate to drain cells may be calculated as the product of the conductance and the head gradient (h-d), where h is the hydraulic head and d is the elevation of the drain. We integrated this relation at the whole sewer network scale, which allows estimating a global conductance associated with the total groundwater flowrate in sewers. Assuming that the baseflow in sewers is only due to groundwater infiltrations, the application of this relation for the high water level period during year 2007 and for the total length of flooded sewer pipes generates a conductance of ~0.9 m2/day. This parameter only changes the water flows and its effect on groundwater levels is less significant. In the end, the specific storage slightly affects the modeling results, as it ranges theoretically from 3.3 × 10−6 m−1 (rock) to 2 × 102 m−1 (plastic clay) [51]. A medium value of 1 × 10−5 m−1 was adopted initially.

#### 2.2.2. Urban Groundwater Modeling at the Catchment Scale and Application to the Case Study

The MODFLOW model is used on the Pin Sec catchment to assess soil water-sewer interactions at the urban catchment scale. Furthermore, certain simplifications not considered in the sensitivity analysis for the unit element were adopted. The two soil layers are represented in this way: (1) the basement layer is located in the entire catchment at a mean depth of 18 m, and (2) the topsoil layer is located in the valley area around the old stream bed, with a mean depth of 2 m. Because both the stormwater and wastewater systems are often placed within the same trench below the street surface, both networks have been combined to represent the draining trench in the soil. Thus, the field drain represents both stormwater and wastewater networks.

Boundary conditions are stipulated as follows and summarized on Figure 5: (1) a zero-flux condition is defined in the upstream boundary of the hydrogeological catchment by considering no-flow cells and (2) a river boundary condition representing the Gohards stream is defined as the downstream condition by using the River package (RIV) on a 200 m long stream; it is a head-dependant boundary condition. The combined 'field drain', representing the sewer systems, drains into this stream. A one hydrological year simulation in transient conditions is performed on the catchment using a daily time step. The transient simulation is initialized using the steady state case calibrated to the field data observed the first day of simulation. The recharge is calculated on a daily basis using the urban surface hydrological model URBS-MO, already successfully applied to urban catchments in the area to simulate the various components of the water budget [36]. The recharge is assumed to be uniform throughout the catchment, and equal to the infiltration calculated by URBS-MO (i.e., rainfall intensity minus both the actual evapotranspiration and surface runoff). Rainfall and potential evapotranspiration data available for the Pin Sec catchment were used with URBS-MO. The initial set of parameters for the simulation is deduced from the observations and a first steady state simulation conducted on the catchment (Table 1). Specific yield and specific storage values are the same as those used for the sensitivity analysis; they vary in space along with the two geological layers defined in this catchment. Due to the lack of information about the spatial distribution of cracks and faulty sealing joints in sewers and to an homogeneous age of the pipes, the drain conductance is assumed to be uniform and the initial value is deduced from the groundwater baseflow observed in the catchment sewers presented in next section.

**Figure 5.** Main boundary conditions applied to the catchment scale model. The 4 × 4 m grid is not perceptible at this scale; the relief is represented with gray elevation contours (altitude asl in m).

The model simulates (1) the hydraulic head and groundwater flux at any grid cell of the domain, and (2) the drain flux at any point of the field drain system. Special attention is paid to the hydraulic head at field piezometers and the flow at the outlet of the catchment, where a flow gauge is located

(Figure 1). The period going from September 2006 to September 2007 is used for calibration, while the period between September 2007 and September 2008 is used for validation.

**Table 1.** Parameters used within the catchment modelling for the different layers.

#### **3. Results**

#### *3.1. Groundwater Dynamics and Experimental Relationship Between Groundwater Level and Baseflow*

The temporal dynamics of groundwater levels are presented in Figure 6 for two of the piezometers; more details may be found in [39]. The mean difference between the high and low water table was approximately 1.7 m during the period 2006–2007, and 1.5 m during the period 2007–2008. This difference varied from ~3 m upstream to 1 m downstream (Table 2).

**Figure 6.** Daily dynamics of the groundwater levels in piezometer PZCS (black) and PZD (grey) from 2006 to 2008 on the Pin Sec catchment. The corresponding surface elevations are plotted with black and grey dotted lines. Bars correspond to monthly rainfall.

Groundwater infiltration in sewer networks depended on several factors such as the location and shape of the water table, the depth of the sewers, and their conditions (e.g., presence of defects). [9] collected information in a small urban catchment in the metropolitan area where the Pin Sec catchment is located, and noticed that groundwater drainage flow became significant once the water table reached a depth of 1.5 m below the surface, and increased as the water table rose. In fact, the mean sewer depth in the area is ~1.2 m below the surface.

**Table 2.** Main features of the groundwater levels during the studied period (September, 2006 to September 2008). The piezometers are presented from upstream (left) to downstream (right). *Zsoil* is the altitude of the corresponding piezometer; 'gw depth' is the groundwater depth from the surface level; R2(SW) and R2(WW) are the determination coefficient of the polynomial regression relationship for stormwater and wastewater baseflow, respectively (2006–2007 period) (See Figure 7 for piezometer PZPS).

**Figure 7.** Relationship between daily groundwater drainage flow *Q* in (**a**) stormwater and (**b**) wastewater sewers, and the groundwater level *H* for the 2006–2007 period. The groundwater level is observed on piezometer PZPS, located on the main sewer intersection. Polynomial regression relationships are y = 0.0112x2 − 0.3693x + 3.0460 (R2(SW) = 0.817) for (**a**) stormwater and y = 0.0083x2 − 0.2711x + 2.2056 (R2(WW) = 0.822) for (**b**) wastewater.

The connection between the average groundwater level of the catchment H and daily baseflow Q in the sewer systems is evident from Figure 7. A specific piezometer (PZPS) located near the main sewer intersection and representative of the average behaviour of the groundwater level dynamics has been chosen; it is considered to be in the catchment center according to groundwater main flow directions. Figure 7a shows a Q~H2 relationship for the baseflow in the storm sewer, which is typical of ideal field drains used in rural hydrology [52]. Groundwater drainage is observed at the outlet of the stormwater

network during the entire observation period, which demonstrates that sewer infiltrations take place all year round. The behavior is similar for the wastewater sewers (Figure 7b). The wastewater baseflow barely varies during summer when the water table is low (between July and October). Baseflow becomes substantial again when the groundwater level for this piezometer PZPS exceeds a threshold value of ~16.75 m (i.e., 3.2 m below the surface). Waste- and storm-baseflows differ because the drainage density of the wastewater system is higher than that of the stormwater system. This is the reason why the Q~H2 relationships are more consistent for the wastewater baseflow that the stormwater baseflow, as revealed by the determination coefficient variation R2(SW) or R2(WW) (Table 2). Additionally, these relationships are generally more pronounced for the downstream piezometers, closer to the baseflow monitoring locations.

This observation is confirmed by the map showing the sewer receiving groundwater contributions during the high water table period of 2007 (Figure 8). This map was realized through the superimposition of groundwater level contours (Figure 3) and the sewer depth levels deduced from GIS data. The total sewer length soaked during winter is 2.6 km for the wastewater sewers and 1.1 km for the stormwater sewers. Finally, the annual subsurface water volume drained into both sewer systems deduced from the sum of daily base-flows between September 2006 and August 2007 is 28% (wastewater sewer) and 14% (stormwater sewer) of the total rainfall respectively. This result proves that the base flow discharges in artificial sewer systems can be a major component of the urban water budget, as 42% of the total annual rainfall is drained to the sewer systems.

**Figure 8.** Sewer pipe network affected by groundwater during winter (2007/03/07): (**a**) stormwater network and (**b**) wastewater network. Black and light gray sewer pipes are located 1 m below and above the groundwater table respectively. Dark gray sewer pipes are between them.

#### *3.2. Comparison between Modeling Results and Observed Data at the Catchment Scale*

Three metrics were adopted for the evaluation of the model: the root mean square error (RMSE), a Bias error, and the determination coefficient R2:

$$RMSE = \sqrt{\frac{1}{n} \sum\_{i=1}^{n} [V\_2(t) - V\_1(t)]^2} \tag{1}$$

$$Bias\,\,error = \frac{\sum\_{i=1}^{n} V\_2(t) - V\_1(t)}{\sum\_{i=1}^{n} V\_1(t)}\tag{2}$$

$$R^2 = \frac{\left(\sum\_{i=1}^n \left[V\_2(t) - \overline{V\_2}\right] \left[V\_2(t) - \overline{V\_1}\right]\right)^2}{\sum\_{i=1}^n \left(V\_1(t) - \overline{V\_2}\right)^2 \sum\_{i=1}^n \left(V\_1(t) - \overline{V\_1}\right)^2} \tag{3}$$

where *n* is the number of time steps involved, *V*1 and *V*2 are the observed and simulated values, *V*1 and *V*2 the temporal average of *V*1 and *V*2.

Absolute error (m) was used as well to estimate the difference between simulated and observed groundwater levels in specific piezometers.

#### 3.2.1. Implementation and Calibration of MODFLOW at the Catchment Scale

In the MODFLOW model of the catchment, we first calibrated the hydraulic conductivity, specific yield and specific storage of each soil type by minimizing the global bias error between observed and simulated groundwater levels. Then we calibrated the drain conductance by maximizing the coefficient of determination R2 between observed and simulated combined baseflows (i.e., stormwater + wastewater sewers baseflows) at the outlet of both sewer systems. The final calibration presented in Table 1 is achieved through a trial and error method. The obtained hydraulic conductivities are two orders of magnitude higher than the observed hydraulic conductivities using the piezometers. This result is consistent with previous statements following the application of other models to ditch-groundwater interactions [53], and confirms the difficulty to simulate groundwater flows and groundwater levels using a few specific point observations of the hydraulic conductivities. High soil heterogeneity, particularly in urban areas, cannot be well characterized using soil properties measurements sampled in a few locations, which leads to a poor simulation of the general groundwater behavior. Despite the differences between observed and calibrated conductivities, the observed differences between the different types of soil are consistent with the theory. Hence, the highest hydraulic conductivity is obtained for alluvial deposits, and the smallest one for mica-schist. Finally, the calibrated drain conductance is higher than initially estimated, which allows a better simulation of high values and the dynamics of the baseflow.

In order to assess the influence of the sewer system on the evolution of the groundwater levels, a complementary simulation of the catchment was done on the same period and with the same set of parameters by removing the sewer system and the corresponding field drain. This simulation exercise shows that the groundwater reaches the soil surface during the winter 2007–2008 on a significant part of the basin; that result clearly confirmed the importance of the drainage ensured by the sewer system [39].

#### 3.2.2. Groundwater Level Distribution Assessment

The temporal and spatial simulation of groundwater levels is assessed using available piezometric records. Figure 9 illustrates the simulation of the isopiezometric contours for the high water table period (7th March, 2007), and can be compared to the observed isopiezometric contours at the same date shown in Figure 3. At this period, the absolute errors between simulated and observed values may be larger than 2 m in some piezometers, especially on the south catchment boundary. They are however better in the middle of the catchment and near the outlet. In addition, absolute errors are higher around the piezometers distant from the sewer trenches. Note that the simulated contours in Figure 9 are plotted on a surface area larger than the experimental one, as the simulated catchment exceeds the piezometric set boundary. Consequently, the simulated flow directions differ significantly from observed ones near the urban catchment boundary, but are rather similar within the urban catchment. Nonetheless, at this time the simulated saturation level is often higher than the observed one. The simulated vs. observed groundwater level dynamics is assessed by computing the RMSE, bias error and R2 for all the piezometer locations, both for the calibration and validation periods (Table 3). R2 values are often smaller than 0.5 for the calibration stage, but increase for the validation stage. In general, the model underestimates the groundwater level for half of the piezometers (PZCPS, PZCS, PZGO, PZUV, PZJV) and overestimates it for the others. The boxplot analysis of the bias error (Figure 10) shows that the underestimation is mainly concentrated on the upstream part of the catchment (PZCRI, PZGO, PZCS, PZJV, PZUV), and might be explained by either the upstream boundary condition or the uniform recharge assumption. Indeed, the spatial analysis of the imperviousness coefficient of the Pin Sec catchment reveals that underestimated piezometric observations are located in more pervious areas within the catchment (i.e., green public parks). Hence, using a spatially distributed recharge could allow a more accurate modeling of the groundwater levels dynamics in these areas. Overall, the simulated groundwater levels compared properly with observed ones both for the calibration and validation periods (Figure 11), despite an underestimation of the model compared to the observed data at the beginning of the validation period, when a particularly dry autumn took place. In the absence of rain, the groundwater model simulates a groundwater level decrease higher than observed, which to some extent questions both the upstream boundary condition assumption and the recharge estimation. Moreover, the quite stable groundwater level observed during this period could be explained by a groundwater influx coming from upstream, which would disprove the assumed upstream boundary condition. Nonetheless, no observations outside the catchment are available to clarify this issue. Finally, drinking water leakage in the soil can occur, which would increase the recharge throughout the year and prevent the high decrease of groundwater levels observed in the simulation; this phenomenon has been discussed in several cities like Bucharest, Romania [4,49], St Louis, Missouri(US) [54] or Santiago, Chile [5].

**Table 3.** Metrics of simulated groundwater levels simulated in each piezometer for the calibration and validation stages. Data availability indicates the percentage of valid data available during the simulation period. Columns *in italics* denote piezometers with poor data availability.

**Figure 9.** Groundwater level contours (hydraulic heads in m) simulated in transient state during the low groundwater table period (7th March, 2007) in the Pin Sec catchment. Absolute errors are represented by red and blue bars at each piezometer. Bar sizes show a variation from −2.66 m (PZUV) to 2.11 m (PZCRI).

**Figure 10.** Box-plot characterizing the absolute error distribution between simulated and observed groundwater levels for the validation period (2007–2008). The piezometers are presented from upstream (left) to downstream (right). Small circles are minimum and maximum values red mark are the median, whereas the ends of the whiskers represent the 25th and 75th percentiles.

Figure 9 shows a small curvature of the simulated groundwater level contours near the sewer trenches caused by the local groundwater drawdown. Groundwater decrease is quite significant in places where sewer density is high (i.e., central and north portions of the catchment). This phenomenon is enhanced when the water table rises above the sewer system, and becomes less significant in summer. For simulated data samples, groundwater drawdown is higher during the wet year 2006–2007 with a decrease of about 0.3 m near the sewer trench, and a zone of influence ranging up to 120 m away from the trench axis.

**Figure 11.** Comparison between simulated (dotted line) and observed (bold line) average hydraulic heads in the catchment from September 2006 to August 2008.

#### 3.2.3. Baseflow Rate Assessment

The simulated soil water flux combines two components simulated in MODFLOW: groundwater flux within the soil and groundwater drainage inside the sewer trench. Both components flow in the same general direction forced by the catchment topography and the downstream condition imposed by the Gohards stream. Groundwater drainage may be compared with the baseflow discharge observed in both waste- and storm- sewer systems during the 2006–2007 simulation (note that the period 2007–2008 is not presented due to the lack of data). Figure 12 shows that the daily groundwater drainage component varies significantly, ranging from 50 m3/day during very dry periods (i.e., the end of summer 2007) to 2600 m3/day during the particularly humid 2006–2007 winter. Indeed, ~630 mm of rain (i.e., ~80% of the catchment mean annual rainfall) fell between November 2006 and April 2007 in the area. Unfortunately, a direct comparison between the simulated groundwater drainage flow and observed baseflow is difficult because of the lack of flow rate data. Considering the data actually available, the comparison focuses on the dynamic evolution of the baseflow. Simulated groundwater drainage and in-situ observations over the simulated period vary in a quite similar manner. Nonetheless the overestimation of the observed flows could not be reduced by further calibration, as drain conductance variations do not affect groundwater drainage significantly. As seen in the sensitivity analysis, soil hydraulic conductivities have more impact on the simulation than drain conductance. Because soil hydraulic conductivities affect not only groundwater drainage but also groundwater level, a simultaneous enhancement of both is not possible.

Special attention must be paid to the relationship between baseflow rates and groundwater levels. As discussed above and shown in Figure 7, both are strongly linked through a power-law function. Figure 13 overlaps simulated and observed values when considering the combined baseflow, and shows that the model can satisfactorily reproduce this groundwater level-baseflow relationship despite of the baseflow overestimation during high water table periods. Although only 32% of the experimental results during the 2006–2007 period are available simultaneously for both sewers, and apart from the low groundwater levels occurring in summer and discussed above, the validity of the relationship holds quite well.

**Figure 12.** Comparison between simulated (dotted line) and observed (bold line) daily combined baseflow at the catchment outlet from September 2006 to September 2007. The recharge evolution shows the infiltration and evaporation periods during this simulation.

**Figure 13.** Combined daily baseflow in sewers vs. piezometric levels between September 2006 and September 2007. Both observed values (black points) and simulated values (grey points) are shown in the (0–0.02 m3/s) range of flowrates, along with their respective polynomial regression relationships: y = 0.0225x2 − 0.7404x + 6.1016 (R2 = 0.803) for the observed values and y = 0.0276x2 − 0.9060x + 7.4436 (R2 = 0.908) for the simulated values.

#### **4. Conclusions**

This study presents an experimental analysis and modeling of the interactions between groundwater and urban sewer systems in a 31 ha urban catchment located near Nantes, France. Groundwater contributions to flow in stormwater and wastewater systems due to sewers' defects were deduced from observed flow rates, whereas piezometric records were used to characterize the dynamics of groundwater levels in the catchment.

Overall, the experimental analysis reveals a strong co-fluctuation of groundwater levels and sewer baseflow, which is more pronounced in downstream piezometers. This study highlights that two or three piezometers located near the catchment outlet are likely adequate to estimate the groundwater-baseflow relationship in an urban catchment. The total baseflow volume in wastewater sewers is larger than in stormwater sewers because of its larger coverage. On an annual basis, the total volume of soil water drained by both sewers is 42% of the total rainfall.

A MODFLOW model was implemented in the study catchment to simulate the groundwater dynamics and contribution to the sewer systems on a daily basis. This application proves the model ability to represent satisfactorily the spatial and temporal evolutions of groundwater levels, although low groundwater levels were underestimated. Spatial differences between modeled and observed groundwater levels in some locations could be explained by the spatial variability of the land uses and the urban soil configuration around each piezometer. The simulation of combined waste and storm water sewer baseflows is not entirely relevant because the model fails to reproduce the high values of baseflow. The calibration of the model lead to hydraulic conductivities values higher than those observed in the piezometers. In fact, the calibrated values are similar to those reported in the literature, and they are attributed to possible secondary porosities and hydraulic conductivity distributions comparable to those of a karstic system [55]. Hydraulic conductivity varies significantly in urban soils due to differences in their characteristics, the existence of different land-uses, and the presence of previous buried constructions.

More work is needed to improve our knowledge of urban soil water flow paths. Baseflow must be better measured, as flow meter sensors are not always reliable to be used with low flows. In addition, more observations at a local scale are needed to better understand groundwater-sewer interactions and improve large scale modeling efforts. Such data will come from either (1) physical models relying on better estimations of soil parameters, which can represent the influence area of a sewer trench, or (2) sensors located near the sewers. Groundwater models applied to urban catchments will require using a graded grid size with mesh refinements near sewers and more urbanized areas. Furthermore, better spatially distributed recharge estimations can also improve the performance of these models. This could be done by coupling the groundwater and surface water models, adopting better modeling frameworks like the modeling chain proposed by [56], or implementing other specific coupling modeling approaches. Finally, a better understanding of the occurrence and spatial distribution of leaks from drinking water pipes is also needed to improve the characterization of the recharge.

**Author Contributions:** This study was designed by F.R., A.-L.L.D., and H.A. The manuscript was prepared by F.R. and revised by H.A. and J.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study received financial support from the French research programs ANR and INSU-EC2CO within the project Hy2ville, and from the Région Pays de la Loire. Funding from CONICYT/FONDAP grants 15110017 and 15110020 were also available.

**Acknowledgments:** We thank Nantes Métropole (Direction de la Géographie et de l'Observation) for providing the GIS data.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

# *Article* **Nitrogen Mass Balance and Pressure Impact Model Applied to an Urban Aquifer**

#### **Mitja Janža 1,\*, Joerg Prestor 1, Simona Pestotnik 1 and Brigita Jamnik 2**

Received: 13 February 2020; Accepted: 16 April 2020; Published: 19 April 2020

**Abstract:** The assurance of drinking water supply is one of the biggest emerging global challenges, especially in urban areas. In this respect, groundwater and its management in the urban environment are gaining importance. This paper presents the modeling of nitrogen load from the leaky sewer system and from agriculture and the impact of this pressure on the groundwater quality (nitrate concentration) in the urban aquifer located beneath the City of Ljubljana. The estimated total nitrogen load in the model area of 58 km2 is 334 ton/year, 38% arising from the leaky sewer system and 62% from agriculture. This load was used as input into the groundwater solute transport model to simulate the distribution of nitrate concentration in the aquifer. The modeled nitrate concentrations at the observation locations were found to be on average slightly lower (2.7 mg/L) than observed, and in general reflected the observed contamination pattern. The ability of the presented model to relate and quantify the impact of pressures from different contamination sources on groundwater quality can be beneficially used for the planning and optimization of groundwater management measures for the improvement of groundwater quality.

**Keywords:**urban hydrogeology; groundwater quality; sewer system; agriculture; groundwatermodeling

#### **1. Introduction**

In 2018, 55% of the world's population lived in urban areas, and this proportion is expected to increase to 68% by 2050 [1]. Europe, with 74% of its population living in urban areas in 2018, is expected to reach an urban population of 80% in 2040 and approaching 85% by 2050 [1]. Considering this population growth in light of climate change projections, the estimates show that with current practices, the world will face a 40% shortfall between the forecasted demand and available supply of water by 2030. Water security is thus emerging as one of the biggest global risks in terms of developmental impact [2].

Groundwater will have an important role, as it represents over 95% of the world's available freshwater reserves [3]. Compared to surface water, it is better protected and a more stable water resource. Better quality, proximity of the resource, and the relatively low cost of pumping wells make groundwater an increasingly important source for city populations in the future.

Aside from socio-economic benefits, urbanization has also created environmental problems, manifested in the deterioration of air and water quality as well as a significant change in water balance [4–7]. Urban subsurface infrastructures strongly influence the chemical as well as the quantitative status of groundwater. The subsurface infrastructure of Bucharest (sewer system and water supply network) contributes more than 71% of the groundwater recharge. In addition, more than 77% of the groundwater discharge is associated with non-natural sinks (mainly the sewer system and the drain of the main sewer collector) [4]. In the Baltimore metropolitan area, the infiltration

of groundwater into wastewater pipelines was found to be the most important factor of urban development affecting the groundwater storage, which can cause its decrease by 11.1% [8].

The assurance of good quality groundwater in the urban environment is a complex task, due to numerous human activities and the number of potential sources of contamination [9,10]. Additionally, the limited ability for observation and the delay between the contamination and its detection make it complicated to define the causal relationship between human activities that exert pressures and their impacts on groundwater quality.

Wastewater exfiltration from the sewer system is an important factor influencing groundwater quality, with its impact being strongly variable both in space and time [11–14]. A review of exfiltration modeling approaches and measurement techniques in terms of both volume and pollutant load showed a wide range of exfiltration rates estimated on the catchment scale, ranging from 1% to 56% of dry-weather flow (DWF) [15].

Nitrate is one of the key groundwater contaminants, and measurements of its concentration are regularly included in groundwater quality monitoring. In Slovenia, about 37% of the groundwater in alluvial aquifers has a poor chemical status, mainly due to high concentrations of nitrate (above 50 mg/L) [16]. Background concentrations of nitrate largely depend on the geological properties of aquifers; in Slovenia, they are below 10 mg/L [17]. Increased concentrations are related to human activities, in particular to agriculture. The main source of diffuse nitrogen input on agricultural land in Slovenia are mineral fertilizers and livestock manure, which account for 84% of the total input of nitrogen [18]. Among the important non-agricultural sources of nitrate in groundwater, wastewater has been identified as one of the major sources of nitrogen in urban aquifers [19].

The impact of nitrogen sources on groundwater have been investigated in a number of studies. Different sources of nitrogen and their contribution to the nitrate content in the groundwater have been analyzed [20–24]. The research on the role of land use [25], vulnerability of groundwater [26,27], processes of leaching [28], and the migration and transformation of nitrogen in groundwater [29–31] have furthered the knowledge on factors controlling groundwater nitrate contamination and its mitigation. The spatial and temporal variation of nitrate in groundwater has been commonly simulated with numerical groundwater flow and transport models [24,29,32,33]. A fundamental step in the development of such models is the estimation of nitrogen loading. The broad range of nitrogen sources in urban areas makes reliable estimation of nitrogen load a difficult task in these areas [19].

The drinking water supply of the City of Ljubljana relies on groundwater, primarily on its abstraction from the Ljubljansko polje aquifer, situated beneath the city. The well fields Kleˇce (55%), Hrastje (12%), and Šentvid (9%) (Figure 1) assure around three-quarters of all the water needs in Ljubljana's drinking water supply system [34]. To determine the urban contamination in the Ljubljansko polje alluvial aquifer, an urban water balance modeling approach was applied [35]. The study focused on an area (approximately 0.76 km2) located within the City of Ljubljana, and indicated that residential land use in an urban area with a high unsaturated thickness may have a significantly smaller impact on groundwater quality than agriculture or industry.

The first catchment scale pressure impact model for the Ljubljansko polje aquifer was developed in the frame of the INCOME project [36]. The aim of the model was to simulate the concentrations in groundwater for a number of pollutants. The model was able to explain practically the entire mass of nitrogen in the sewer system by the known sources. On the basis of these results and an encouraging decrease of observed nitrate concentration in groundwater, more optimistic groundwater management objectives for the year 2017 were set. The targeted average nitrate concentration in the aquifer was set below 15 mg/L and in all observation wells below 25 mg/L [37]. Nitrogen input from agriculture was identified as the most important source of uncertainty in the model. The explanation of high nitrate concentration at the local level, that is, individual observation wells, was pointed out as an important groundwater management activity.

In the present study, a refined and updated groundwater flow and transport model [38] was applied in the catchment scale approach. The aim of the modeling was to analyze the causal relationship between pressure (nitrogen load) arising from the leaky sewer system and from agriculture, and its impact on the quality of groundwater (nitrate concentration) in the Ljubljansko polje aquifer [39]. Nitrogen load from the leaky sewer system was estimated with the empirical model of spatially distributed sewer exfiltration [40] and together with estimated nitrogen surpluses from agriculture [41] used as input into the groundwater transport model. Modeled and observed nitrate concentrations in groundwater were analyzed to assess the uncertainty of the model and to identify priority areas for remediation.

**Figure 1.** Map of the study area showing observation locations, groundwater level isolines, model boundaries, artificial surfaces, and well fields.

#### **2. Materials and Methods**

#### *2.1. Study Area*

The study area (57.91 km2) is located in the flat part of the catchment of the alluvial Ljubljansko polje aquifer forming the subsurface of the City of Ljubljana, the capital of Slovenia. It corresponds to the area of the local groundwater model (Figure 1). The mean altitude of the area is 292 m above sea level (m a.s.l.) and prevailing land cover types are artificial surfaces (59%) and agricultural areas (39%) [42].

The Ljubljansko polje alluvial aquifer is developed in a tectonic basin, filled with highly permeable Quaternary gravel and sand beds, which are partly conglomerated (cemented with calcite cement) [43]. The thickness of the sedimentary fill is up to 100 m down to its very low permeable pre-Quaternary basement, composed of Carboniferous–Permian siliciclastic rocks. The groundwater table is on average 25 m below the surface. The unconfined intergranular aquifer with high transmissivity (Figure 2) is mainly recharged from the Sava River in the northwestern part of the study area. It drains into the river in the northwestern part, which largely predisposes the groundwater flow direction (Figure 1). The main water balance components of the regional model area (Figure 1) for the 3 year modeling period (2013–2015) are presented in Table 1.

**Figure 2.** The calibrated transmissivity of the aquifer.

**Table 1.** Water balance components (modeling period 2013–2015).

#### *2.2. Nitrogen Mass Balance and Pressure Impact Model*

The applied methodology [36,37,40] comprises the following main steps: (1) estimation of the spatial distribution of nitrogen load from the leaky sewer system and agriculture (nitrogen mass balance model), (2) use of the nitrogen load as input into the groundwater transport model, and (3) analysis of the modeling results with observed nitrate concentrations in the observation wells.

#### 2.2.1. Nitrogen Load from the Leaky Sewer System

Assessment of nitrogen load from the leaky sewer system was based on the model of spatially distributed sewer exfiltration developed by Prestor et al. (2011) [36] for the sewer system of the City of Ljubljana. The main steps of workflow for the estimation of nitrogen load from the sewer system are presented in Figure 3. Firstly, the exfiltration rate (losses) from the sewer system was determined by the mass balance of wastewater released to sewer system and outflow from the sewer system. Then, the total mass of pollutant in the sewer system was calculated. In the next step, the exfiltration rate from the tested section of the sewer system was assessed and then upscaled to the spatial distribution of exfiltration of the entire sewer system.

**Figure 3.** The flowchart of the workflow for the estimation of nitrogen (N) load from the sewer system (adapted after Pestotnik et al., 2019 [40]).

Exfiltration Rate (Losses) from the Sewer System

Daily flow rate measurements at the wastewater treatment plant (WWTP) for the period 2000–2010 were compared with in-house water consumption for that period, provided by the drinking water supply and the wastewater system operator. In the analysis, 13 DWF measurements were used, mainly in the winter period, when rainfall infiltration and outdoor water use have a low influence on the sewer system water balance. Nevertheless, the measured flow rate at the WWTP exceeded the quantity of in-house water consumption, which indicated the presence of extraneous water. Estimation of the extraneous water was based on the assumption that 0.01 km2 of surface in the catchment of WWTP contributes 0.05 L/s of inflow of water [44], which accounted for the entire catchment 276.6 L/s. On the basis of this, the exfiltration rate from the entire sewer system in the catchment of WWTP was estimated to be 168 L/s (on average 0.21 L/s/km) [40].

#### Mass of Nitrogen in Wastewater

The estimation of the mass of nitrogen in water is based on the mass balance model [36,37]. The aim of this model was to investigate the shares of known and unknown sources of different pollutants (e.g., nitrogen, phosphorus, chloride, sulphate, copper, chromium, lead, organic carbon) in the sewer system. The mass of pollutants released into sewer system was calculated in the model as a sum of inputs from households, from roads and other traffic surfaces, and from industry. The model results were compared with two extensive analyses of chemical composition of wastewater at WWTP, made in 2008. More frequent measurements were performed only for nitrogen. The concentration of total nitrogen was correlated with discharge and both varied over time depending on the amount of rainfall (Figure 4). The yearly mass of other pollutants discharged at WWTP was a product of discharge and concentrations defined in the two mentioned analyses. Data from the first analysis was used as representative for the dry-weather periods (121 days), and from the second measurement as representative for the wet periods (121 days).

Good agreement between the modeled and observed content of phosphorous (−2.1%), nitrogen (+5.9%), and lead (13.1%) was found [36,37,40]. The content of other pollutants could not be well explained by the known sources.

**Figure 4.** Concentration of total nitrogen at wastewater treatment plant (WWTP) presented with brown dots and 5 days moving average (brown line), discharge presented with blue dots and 5 days moving average (blue line), and rainfall high (gray line) measured at Ljubljana Bežigrad weather station (year 2008).

Spatial Distribution of Sewer Exfiltration

The model of spatial distribution of sewer exfiltration [36] was based on a detailed inspection of a 14.82 km section of the sewer system, which represents 2% of the entire sewer system. Closed-circuit television (CCTV) provided by the sewer system operator was used to detect the pattern of sewer damages in this section. Among 1246 detected damages, 104 were selected as significant and characterized in more detail to quantify exfiltration. This pattern and information on type of material, age, and size of pipes [45] were used to estimate the spatial distribution of the exfiltration and nitrogen load from the entire sewer system (Table 2).

**Table 2.** Estimated exfiltration and nitrogen load from pipelines of the entire sewer system [40].

For the calculation of the exfiltration flow rate from the leaky areas, an approach based on Darcy's law was used. It assumes linear dependency between the hydraulic gradient (I), the hydraulic conductivity of the clogging layer (kf), the leakage area of the pipe (Aleak), and the exfiltration flow rate (Qex) [46]:

$$\mathbf{Q\_{c\alpha}} = \mathbf{k\_f} \times \mathbf{A\_{leak}} \times \mathbf{I} \tag{1}$$

where the hydraulic gradient (I) is the ratio between the sum of the water level inside the sewer (hs) and the pipe wall thickness (hw), and the pipe wall thickness (hw):

$$\mathbf{H} = (\mathbf{h}\_{\rm s} + \mathbf{h}\_{\rm w}) / \mathbf{h}\_{\rm w} \tag{2}$$

In the calculation, three different hydraulic conditions in the pipelines were considered: (1) dry-weather conditions, days with mainly waste water discharge from water consumers (DWF); (2) days with significant rainfall inflow; and (3) extreme storms events with very high rainfall inflow. An analysis of rainfall data of a 10 year period (2000–2010) was performed. It showed that on average,

in 1 year, DWF conditions last for 146 days, significant rainfall inflow occurs in 218 days, and once a year extreme storms are expected. It was assumed that sewer pipes are 10% full in DWF period, 50% in significant inflow period, and 100% in storm events (Figure 5). These data were used together with data on the diameter of pipes and a representative sample of damages for different materials and ages of pipes to calculate the spatial distribution of yearly exfiltration (Figure 6).

**Figure 5.** Scheme showing the different degrees (10%, 50%, 100%) of fullness in sewer pipes.

**Figure 6.** The spatial distribution of the estimated nitrogen load from the leaky sewer system (in kg/day per 50 m grid cell).

#### 2.2.2. Nitrogen Load from Agriculture

Nitrogen load from agriculture used in the model was the nitrogen surplus calculated as the difference between nitrogen added to an agricultural land (mineral fertilizers, livestock manure, atmospheric deposition, and mineralization of organic matter) and nitrogen removed from the agricultural land (removed field crops and denitrification in the soil) [41]. On the basis of literature and statistical data, the nitrogen mass balance for different agricultural land use areas (arable fields, vineyards, orchards, and meadows and pastures) was calculated, and the results were provided in the polygon information layer. These data were used to generate a model input grid (Figure 7).

**Figure 7.** The spatial distribution of the estimated nitrogen load from agriculture (in kg/day per 50 m grid cell).

#### 2.2.3. Groundwater Transport Model

Groundwater flow and transport simulations were carried out using the MIKE SHE/MIKE 11 modeling framework [47], which enabled transient simulation of the most important processes of the hydrological cycle in the study area. The solute transport model used in the study [38] is a refined and local version of the calibrated regional hydrological model [48]. The regional model (Figure 1) is a spatially distributed (200 m grid cells) and transient model with implemented daily values of time-dependent variables (temperature, precipitation, river levels and discharges at boundaries, groundwater level observations, abstractions, and leakage from the water supply system). The calibration and validation of the model was based on a comparison of the observed and simulated groundwater heads at the locations of 12 observation wells distributed over the entire study area, and the Sava River levels at one gauging station [48]. Additionally, the solute transport model was validated against the observed concentrations of the trichloroethylene pollution plume discovered in 2004 [49].

The local model (Figure 1) has one calculation layer and at boundaries uses hydraulic conditions from the regional model. Its model domain is discretized into horizontal grid cells of the size 50 × 50 m. The estimated nitrogen load was used as a spatially distributed, subsurface and constant mass input into the saturated zone of the aquifer. Neither retardation nor degradation processes were considered in the simulations. For the nitrate background concentration, 6 mg/L was used. This concentration is comparable to the lowest observed nitrate concentration (7.2 mg/L in observation location 18, Table S1 and is slightly higher than the observed concentration (4.9 mg/L) in the Sava River [50].

Scenarios for a 6 year simulation period (2010–2015) were performed. The last 3 years (2013–2015) of simulations were considered for the analysis. The first half of the simulation period was used for the model warm-up, which assured the avoidance of biased results affected by initial conditions.

#### **3. Results**

The estimated sources of nitrogen from the sewer system and from agriculture are distributed over an area of 19.84 km2 (34%) and 27.81 km2 (48% of model area), respectively. The mean nitrogen load in these model cells from the sewer system is 45 g/day and from agriculture is 52 g/day (Table 3). In the local model area of 58 km2 , this accounts for a total nitrogen of 334.22 ton/year, 38% arising from the leaky sewer system (Figure 6) and 62% from agriculture (Figure 7).

**Table 3.** Statistics for the estimated nitrogen load expressed in g/day per 50 m cell.

Nitrate concentration in groundwater was available from 58 observation wells and 2 springs. The observed locations were numbered according to their geographical positions, in accelerating order from the west towards the east (Figure 1). The modeled total nitrate concentration (Figure 8) is the sum of the modeled impacts of nitrogen load from the sewer system and agriculture (Figure S1, Figure S2) with an assumed natural background concentration of 6 mg/L. Median values of observed and modeled concentrations in a 3 year period (2013–2015) at observation locations (obs) are presented in Table S1 and summary statistics in Table 4.

**Figure 8.** The spatial distribution of modeled compared to observed nitrate concentrations in groundwater. **Table 4.** Summary statistics for modeled and observed nitrate concentrations at observation locations (Table S1).

\* Sum of the modeled impacts from the sewer system and agriculture with an assumed natural background of 6 mg/L. \*\* Difference between observed and modeled total concentrations.

The average difference between the observed and modeled concentrations at the observation locations was found to be 2.7 mg/L (Table 4, Figure 9). The distribution of differences was left-skewed (Sk = −0.73, Figure 10). The correlation coefficient between the observed and modeled concentrations was 0.51 (Figure 11).

**Figure 9.** Spatial distribution of differences between the observed and the modeled nitrate concentrations in the groundwater in mg/L with red showing and overestimation by the model.

**Figure 10.** Histogram of differences between the observed and modeled nitrate concentrations at the observation locations.

**Figure 11.** Scatterplot of observed and modeled nitrate concentrations at the observation locations, with fitted line (black line), 95% confidence band (light-brown dotted line), and 95% prediction band (blue dashed line).

#### **4. Discussion**

The modeled spatial distribution of nitrate followed the groundwater flow pattern. It reflected the general characteristics of the observed concentrations, which are higher in the southeastern and central parts of the model area (Figure 8). In spite of higher total nitrogen load from agriculture, the highest nitrate concentrations (>25 mg/L) were observed in areas (obs 2, 3, 4, 6 and 8) with predominant impact from sewer system (Figure 8).

The model gives insight into the impacts of two contamination sources in the capture zones of abstraction wells that are used for drinking water supply. Nitrate concentrations in observation wells in the Šentvid (obs 5 and 7) and Kleˇce well fields (obs 11–17, 19, 20, 21, 23, 25, 28 and 29) were found to be relatively low (on average below 15 mg/L) and, according to the model, arise predominantly from agriculture (72% of anthropogenic impact). Nitrate concentrations in observation wells in the Hrastje well field (obs 50–56) were higher (on average above 20 mg/L), and the shares of nitrogen load from agriculture (47%) and the sewer system (53%) are comparable. These findings give a bigger picture of the impact of different nitrogen sources on groundwater quality in comparison to the local study, focused on a small part of the urban area [35].

On average, the modeled concentrations were found to be 2.7 mg/L lower than the observed concentrations (Table 4, Figure 9). The discrepancy is relatively low and gives a good estimate of the nitrogen mass balance in the study area. Overestimation of observed nitrate concentrations by a factor of 4.5 was reported in the case study solely by using the nitrogen mass balance estimate for infiltration water [22]. This emphasizes the important contribution of the integration of groundwater dynamics in modeling in order to reduce the uncertainty in the estimate of the nitrate concentration in groundwater.

The modeled concentrations were most evidently underestimated in the area of the highest transmissivity and groundwater discharge (Figure 2). This zone stretches from the northeastern part of the model area, in which the most intensive recharge from the Sava River occurs and further follows the groundwater flow path. This discrepancy could be attributed either to the overestimated groundwater discharge rates or to underestimated nitrogen load, which in this area arises mainly from agriculture. In the urbanized part of the study area, underestimated concentrations could be related to unknown non-agricultural sources of nitrogen [19]. On the contrary, the overestimated nitrate concentrations are mainly in the southern part of the model area. This zone stretches from the southeastern part of the model area, in which the highest concentrations were observed (obs 2, 3 and 4) and follows the groundwater flow path towards the southeast.

The differences between the modeled and observed nitrate concentrations are related to the uncertainties arising from different sources. The level of complexity of used models was adjusted to available data and knowledge. The selection of the models followed the principle as simple as possible, as complex as necessary [51,52]. The used estimates of nitrogen load were based on the limited available data and a simplified nitrogen mass balance model. The wastewater exfiltration model used characteristics of 2% of the sewer system and scales up these findings to the entire sewer system. Due to this generalization and the lack of observations and the complexity of exfiltration process [14], the quantification of impact of sewer leaking on groundwater quality is uncertain. The available quantification of nitrogen load (surpluses) from agriculture was based on literature and statistical data [41], and thus it did not consider site-specific agricultural practices and their changes over time.

The groundwater model is a mathematical representation of an aquifer that is highly productive, characterized by a saturated zone up to 70 m deep and a fast groundwater flow with up to 20 m/day. The simulated pollution plumes are narrow and fast progressing, thus difficult to observe [49]. The modeling of the nitrate concentration is based on simulated groundwater discharge and flow direction. These two processes in the aquifer are difficult to measure and the lack of available data limits the possibility of validating the groundwater transport model. Especially in the areas close to the boundaries of the model, in which the transmissivity of the aquifer (Figure 2) is low, modeled concentrations were found to be very high (Figure 8) and sensitive to the interpreted shape of the aquifer basement. It is assumed that in these areas, the uncertainty of the model is the highest.

Data on nitrate concentration in groundwater was acquired from three monitoring networks managed by the Slovenian Environmental Agency, the City of Ljubljana, and the water supply company JP VOKA SNAGA. The samples were taken from observation wells using withdrawal method. No multilevel sampling technique for collecting depth-specific groundwater samples was applied. Representative depth-specific groundwater sampling in unconsolidated and coarse-grained aquifers is difficult to obtain. In such aquifers, even the application of packers within a screened well does not assure representative depth-specific sampling, due to the high risk of vertical circulations in the annular space around the well screen or casing [53]. The available data on groundwater sampling were insufficient for characterization of vertical distribution of nitrate in the aquifer. Due to this, the median value of measurements (2013–2015) in each well was considered as a representative concentration of the whole groundwater column at the observation location. This fact led to the simplification of the model and the use of a single calculation layer, which limited the simulation of vertical stratification and interpretation of the site-specific vertical concentration distribution.

#### **5. Conclusions**

The presented modeling approach enabled relation and quantification of the impact of pressure from contamination sources on groundwater quality. Modeled spatial distribution of nitrate content in groundwater pointed out the leaky sewer system as the prevailing source of nitrogen load inducing locally high nitrate concentration (above 25 mg/L). Reconstruction of the part of the sewer system in the catchments of observation wells with high nitrate content was recognized as a priority measure for the improvement of groundwater quality. The developed model can be used as a supportive tool in the optimization of mitigation measures, providing the simulation of effects of different remediation scenarios in a cost–benefit analysis. This approach could be applied also to assess the efficiency of groundwater management plans and to analyze the contribution of the stakeholders' activities to the improvement of groundwater quality.

Nitrate concentration in groundwater could be well explained by the known sources of nitrogen that approved the applied integrated mass balance and groundwater modeling concept. Application of this concept to other pollutants revealed a high share of unknown sources for several pollutants; among them, chromium (CrIV) was pointed out as the most threatening for the water supply [40]. Identification and remediation of unknown sources of pollutants whose content in groundwater is high, considering drinking water quality standards, or increases in time, is a priority for groundwater management and the assurance of safe drinking water supply for the City of Ljubljana.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4441/12/4/1171/s1, Table S1: Median values of modeled and observed nitrate concentrations (2013–2015), Figure S1: The spatial distribution of the modeled nitrate concentration arising from the sewer system; Figure S2: The spatial distribution of the modeled nitrate concentration arising from agriculture.

**Author Contributions:** Conceptualization, M.J. and J.P.; methodology, J.P., S.P. and M.J.; investigation, J.P., S.P. and M.J.; data curation S.P., J.P. and M.J.; writing—original draft preparation, M.J.; writing—review and editing, M.J., J.P. and B.J.; project administration, J.P., B.J. and S.P.; funding acquisition, J.P., B.J. and M.J. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by projects AMIIGA (Interreg Central Europe, No. CE32), INCOME (LIFE07, ENV/SLO/000725), IGCP 684 (UNESCO-IUGS), and by the Slovenian Research Agency (research core funding Groundwaters and Geochemistry (P1-0020)).

**Acknowledgments:** The authors wish to thank Petra Megliˇc and Dejan Šram for the help in cartography. The authors are very grateful to reviewers for very constructive reviews that significantly improved the manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Managing Stormwater by Accident: A Conceptual Study**

**Carly M. Maas 1,2,\*, William P. Anderson, Jr. 3 and Kristan Cockerill 4**

**Abstract:** Stormwater-driven road salt is a chronic and acute issue for streams in cold, urban environments. One promising approach for reducing the impact of road salt contamination in streams and adjacent aquifers is to allow "accidental wetlands" to flourish in urban areas. These wetlands form naturally as a byproduct of human activities. In this study, we quantified the ability of an accidental wetland in northwestern North Carolina, USA, to reduce the timing and peak concentration of road salt in a stream. Monitoring suggests that flow and transport processes through the wetland reduce peak concentrations and delay their arrival at the adjacent stream. We expand these findings with numerical simulations that model multiple meltwater and summer storm event scenarios. The model output demonstrates that small accidental wetland systems can reduce peak salinities by 94% and delay the arrival of saltwater pulses by 45 days. Our findings indicate that accidental wetlands improve stream water quality and they may also reduce peak temperatures during temperature surges in urban streams. Furthermore, because they find their own niche, accidental wetlands may be more effective than some intentionally constructed wetlands, and provide opportunities to explore managing stormwater by letting nature take its course.

**Keywords:** accidental wetland; road salt; urban hydrogeology; headwater stream

#### **1. Introduction**

The idea of "accidental" wetlands is gaining traction as a potentially valuable component of various ecosystems. These are water bodies that "result from human activities, but are not designed or managed for any specific outcome" [1]. As the moniker implies, these "accidental" systems arise on their own, often in vacant lots or various low spots where stormwater or irrigation runoff collect. In a study of the South Platte River basin, for example, 89% of the extant wetlands exist because of various irrigation conveyances [2]. Accidental wetlands typically form in low-lying abandoned or underutilized landscapes with poor drainage [3]. Flooding often occurs in these areas due to a lack of stormwater management and sediment transport processes, which contribute to repeated sediment deposition, allowing wetland plants to grow and establish a habitat [3]. Brooks et al. [4] documented that in arid, urban areas, effluent from human sources can account for up to 90% of streamflow during dry periods. Any wetland-like areas on those rivers exist because of human activity. Although they are not yet heavily studied, there is evidence that these accidental wetlands offer benefits to rival constructed versions [1,5,6]. In fact, Palta et al. [3] argue that "accidental wetlands may provide more services than designed environments because the latter are commonly over-designed for a limited set of specific functions." There is also evidence that small wetlands, which many accidental wetlands are, can be as or more effective than larger ones [7,8]. Although they are just beginning to be intentionally documented and studied, Palta et al. [3] offer evidence highlighting the ubiquity of accidental wetlands and their potential to offer diverse benefits. In the

**Citation:** Maas, C.M.; Anderson, W.P., Jr.; Cockerill, K. Managing Stormwater by Accident: A Conceptual Study. *Water* **2021**, *13*, 1492. https://doi.org/10.3390/ w13111492

Academic Editor: C. Radu Gogu

Received: 26 April 2021 Accepted: 24 May 2021 Published: 26 May 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

limited literature available, researchers have found that accidental wetlands can remove pollutants, mitigate heat, store groundwater, restore surface-mined lands, and provide the social benefits of additional green space in urban settings [3,9–11].

In these few known studies, accidental wetlands were shown to help with pollutant remediation, and a key pollutant in cold urban areas is road salt. Like accidental wetlands, in recent years researchers have paid increasing attention to de-icing salt in stormwater runoff as a key issue for stream quality, especially in urban areas. Although there are numerous salt sources in urban environments, including water-softening, septic field discharge, natural rock weathering, and wastewater treatment plants [12], the dominant source is road salt applied to the road surface in order to melt snow and ice [13,14]. The annual road salt application in the United States in 2014 was about 56.5 million metric tons, costing about \$1.18 billion [15]. Elevated salt levels in waterways contribute to negative consequences, including contaminating drinking wells [16] and various ecosystem impacts. In their review of the literature, Hintz and Relyea [17] found that salt negatively affects all freshwater species, but that the level at which salt becomes a concern is highly variable across species.

There are complex relationships between surface and groundwater systems, such that salt remains available for long periods of time as it travels between surface and groundwater systems [18–20]. Furthermore, numerous studies have focused on how effective various "green infrastructure" (e.g., rain gardens, bioretention systems) are in cold climates (see Kratky et al. [21] for a review). Retention systems capture salt as it runs off and can slow its path to waterways [18,22]. Retention systems do not, however, eliminate salt or retain it for long periods. Studies have found that salinity remains a key stream stressor even in restored streams and in those with stormwater management efforts in place [23,24]. In addition to its own negative consequences, salt can mobilize various metals and lower plants' uptake of metals, thereby reducing the effectiveness of bioretention [21,25–27]. Cockerill et al. [20], however, did find that retention had the potential to reduce salinity spikes, which can reduce the long-term salt levels in a given hydrologic system. Flattening the curve, as it were, and alleviating those salt spikes, may over time lower the summer pulses of salt into that system. Given Hintz and Relyea's [17] findings on the variability in how much salt triggers negative consequences, reducing salinity spikes and reducing summer inputs may reduce the effects on some aquatic species.

Although researchers are studying accidental wetlands and the negative effects of road salt on urban waterways, as far as we know, these two subjects have not been assessed in concert. We offer a conceptual study addressing the potential for accidental wetlands to ameliorate salinity levels. Because constructed wetlands offer highly variable results depending on scale, lifespan, and local hydrologic conditions [28–31], paying more focused attention to sites where wetlands form themselves may offer a cheaper and equally or more effective alternative means of lessening road salt impacts. Though no wetland will completely remove the salt from an urban stream system, there is still potential to improve urban water quality by documenting and encouraging accidental wetlands.

The project documented here is a testimony to the value of paying attention to local environments and the role of serendipity in research. Copeland [32] describes serendipity in science as "the intersection of chance and wisdom" and notes that the value of the serendipitous process or event can only be assessed in hindsight. In our case, authors Anderson and Cockerill initially noticed a "wetland" forming in a concrete culvert on the Appalachian State University campus and joked about its potential value. As we continued to observe the accumulation of sediment and increased vegetation, we began thinking that the site actually warranted more focused attention. At that point, author Maas was seeking a project and began to intentionally explore the effect of that wetland on road salt contamination for a headwater stream in an urban setting. In searching for background information to guide this more intentional study, we encountered the premise of "accidental wetlands" and recognized that we had one.

Once intentionality was established, this study focused on these research questions:

#### *Site Description*

Boone Creek and the accidental wetland investigated in this study are located in the small yet heavily urbanized town of Boone, North Carolina (Figure 1). Boone Creek is a headwater tributary of the South Fork New River, located within the Blue Ridge Mountains of northwestern North Carolina [20]. The town of Boone has high runoff ratios that result in flashy, flood-prone streams due to heavy urbanization near the stream. Although the total catchment area has 24.3% impervious surface cover, much of the more natural land use occurs on the mountain slopes near the catchment boundaries (Figure 2). Within 25 m of Boone Creek, the impervious surface cover ranges from 1% to 75% [33]. The Town of Boone has about 19,500 residents (U.S. Census Bureau, Quick Facts, https://www.census.gov/ quickfacts/fact/table/boonetownnorthcarolina,US/PST045219, accessed 25 May 2021) and Appalachian State University (ASU) adds an additional 20,000 students to the local population (Appalachian State University Facts, https://www.appstate.edu/about/ accessed on 25 May 2021). The stream is in a mountainous area, with a total catchment relief of nearly 500 m, although the main channel gradient is fairly low, at 2%. The stream has previously had trout populations and retains a state classification as a "trout stream" [34]. Other tributaries to Boone Creek have gradients of greater than 10%. The catchment of the basin considered in this study has an area of 5.2 km2 [20].

**Figure 1.** The Accidental Wetland (AW) site is located in Boone, North Carolina, USA (left top insert map) at 36◦12 54.6 N 81◦40 57.0 W. The solid blue line denotes the surface water of Boone Creek and the dashed blue line signifies culverted sections of Boone Creek. The circles indicate data collection sites, Accidental Wetland (AW) and Varsity Gym (VG). The bottom left insert map shows the Boone Creek watershed (heavy purple line) with the black triangles indicating peak elevations within the watershed. The arrow represents the direction of flow. ASL denotes above sea level. Adapted from ref. [20].

**Figure 2.** Land use map of the Boone Creek watershed. The Accidental Wetland site is located in the light blue and pink regions, which designate medium and high intensity development, respectively. The stream includes culverted and nonculverted segments. The map was created using ArcGIS® software by Esri, Redlands, CA, USA. Accessed on 10 February 2021. ArcGIS® and ArcMap™ are the intellectual property of Esri and are used herein under license. Copyright © Esri. All rights reserved. For more information about Esri® software, please visit www.esri.com, accessed on 10 February 2021.

Boone Creek (Figure 1) runs through several culverts, with the longest being about 600 m long [20]. The area experiences 89 cm of snow per year, requiring 6 × 105 kg (600 metric tons) of deicing salt applied to the impervious surfaces on average annually [12]. Electrical conductivity data, converted to equivalent chloride, indicate that the stream regularly violates the US Environmental Protection Agency's (EPA) chronic threshold of a four-day average chloride concentration, exceeding 230 mg/L, and also exceeds the acute threshold of a one-hour average chloride concentration, exceeding 860 mg/L on a regular basis in the winter salting season [20,35]. Between 1 July 2014 and 1 July 2015, Boone Creek exceeded the EPA's chronic threshold for 10% of the time, or 36 days, and exceeded the acute threshold for the equivalent of 8.5 days [20].

Our accidental wetland site is unique because it was created in a constrained concrete channel and did not form on an existing substrate. Subsequently, there is no direct connection between the wetland and the natural groundwater reservoir. We are treating the accidental wetland as an isolated "aquifer" that is disconnected from the natural groundwater flow system. Throughout the rest of this paper, our references to groundwater flow and solute transport are solely related to this process in the small accidental wetland that formed in the concrete culvert that ultimately feeds water to Boone Creek through seepage at its downgradient end. The only sources of water to the accidental wetland are from runoff from impervious surfaces or direct precipitation. Figure 3 shows the total area of the drainage basin to the accidental wetland, delineating pervious areas and impervious areas, such as buildings and sidewalks. The 30-year normal precipitation in Boone as measured at the BOONE 1 SE, NC US weather station between 1981 and 2010 is 1338 mm, with roughly equal precipitation during each season (https://www.ncdc.noaa.gov/cdo-web/datatools/normals, accessed on 21 April 2021). As Figure 4a–f show, the 600 m culvert containing Boone Creek is open for five meters to

allow stormwater to enter from a concrete conveyance that drains part of campus, as well as part of the Town of Boone.

**Figure 3.** Total area of the drainage basin for the Accidental Wetland (AW) site. Gray areas indicate impervious surfaces such as sidewalks, roads, and rooftops and the white areas represent green spaces.

**Figure 4.** Evolution of the Accidental Wetland (AW) from 2005–2019. The photos were taken on the bridge over the wetland, facing downstream toward Boone Creek. The wetland was removed in September of 2018.

The focus of this study was not on how accidental wetlands form, because the exact timing for this specific accidental wetland was unclear and did not factor in our findings. It may have originated in 2011–2012 from large storm events, carrying significant sediment into the culvert. By 2015, authors Anderson and Cockerill observed that the channel had become a heavily vegetated wetland (Figure 4b), and we began observing it regularly. Up to that point, we had paid little attention to the site and the few photos available before 2015 show spotty vegetation in 2005 (Figure 4a) and we observed some sediment accumulation along one edge by 2012. Since 2015 the system evolved to include diverse vegetation throughout the entire channel and a meandering "stream" through the vegetation.

A storm event in October 2017 deposited a large amount of sand/soil at the culvert end of the wetland. At that time, authors Maas and Anderson decided to assess what was happening more formally. They installed three wells and an electrical conductivity logger (HOBO Conductivity Data Logger, U24-001) along the length of the wetland in 2017. The monitoring wells showed large fluctuations in water levels throughout the year, with high levels during intense storm events during warm months and during the wet and cold winter. In September 2018 the university removed all of the vegetation and accumulated sediment in the accidental wetland as part of flood control measures on campus (Figure 4e,f).

#### **2. Materials and Methods**

#### *2.1. Field Data Collection*

We collected handheld probe measurements in the accidental wetland from 7 September 2017 through 6 February 2020. We used a Yellow Springs Instrument (YSI) 556 MPS (YSI, Yellow Springs, OH, USA) to measure salinity, dissolved oxygen, temperature, and other parameters. These data allowed us to gain an understanding of the salinity distribution within the wetland both spatially and temporally, and guided our conceptual model of transport through the sediment. For examples of studies using the YSI 556 MPS, see [36,37]. We collected handheld measurements on a regular basis at four locations along the wetland and at one location in Boone Creek (Figure 5). The YSI probe measured salinity in PSU, and we converted these data to equivalent chloride concentrations by multiplying by 606.6 mg/L Cl− per PSU [38]. PSU is approximately equivalent to parts per thousand (ppt) of NaCl; therefore, for every 1 g/L of NaCl, 606.6 mg/L of the solution is Cl− and 393.4 mg/L is Na+ [38]. All subsequent data in this study are presented in terms of equivalent chloride concentrations.

Datalogging at the Varsity Gym (VG) monitoring site has collected electrical conductivity (EC) data and stage-discharge data at 15-min intervals since July 2014, approximately 200 m downstream of the wetland (Figure 1). EC is a proxy for salinity [20]. The stream's salinity is less than that of the wetland due to the mixing of runoff after winter storm and precipitation events and the constant input of baseflow, which most often has a lower salinity than winter runoff. These data are scaled to provide a continuous dataset of boundary conditions for the groundwater flow and transport simulations, which are described later in the manuscript.

We used the Hazen method [39] on five soil samples collected from throughout the wetland to analyze grain size distributions. We used these data to estimate the accidental wetland's hydraulic conductivity using

$$\mathbf{K} = \mathbf{C} \left(\mathbf{d}\_{10}\right)^2,\tag{1}$$

where d10 is the diameter of the 10% finer grains (cm) and C is a coefficient based on grain size and sorting. We averaged the samples to estimate a bulk hydraulic conductivity for the wetland, and we further calibrated this value with subsequent numerical modeling experiments.

**Figure 5.** A map of the Accidental Wetland (AW). The stars indicate the locations from which the pictures were taken. Red circles indicate runoff surfaces at the inflow area of the wetland. The sites, denoted by white circles, were named Inflow, Accidental Wetland Site 1 (AW 1), Accidental Wetland Site 2 (AW 2), Seepage face, and Boone Creek. (**a**) Facing upstream towards the stormwater culverts, or the inflow sources. (**b**) Standing on bridge facing downstream towards seepage face and Boone Creek. (**c**) On grass facing Boone Creek stormwater culvert and the seepage face. The diameter of the culvert is 1.5 m.

#### *2.2. Numerical Modeling of Groundwater Flow and Solute Transport*

In this study, we used FEFLOW (DHI-WASY, DHI Group, Horsholm, Denmark) to simulate groundwater flow and solute transport through the wetland, utilizing salinity and runoff data from 2018 to formulate boundary conditions. FEFLOW is a three-dimensional finite-element groundwater flow and solute transport model that is well-documented in the literature and has been used to study subsurface solute transport (e.g., [40,41]). It has also been used in previous studies of the Boone Creek watershed [20]. The numerical model is not an attempt to recreate the empirical values collected from the accidental wetland through our sampling efforts; rather, the goal is to demonstrate with the model the groundwater flow and transport processes that are causing the decrease in salinity along the length of the wetland. The model does not account for any uptake of salt by the wetland plants or short-term retention by the wetland sediment, so the model output should be considered a conservative estimate of the effectiveness of the accidental wetland in delaying the arrival of salt-laden runoff.

Figure 6 shows the model domain that was used in groundwater flow and solute transport simulations for this study. The two-dimensional simulations represent a 21.2-meter-long by 0.2 m cross section of the accidental wetland. The base of the model is a no-flow boundary and represents the interface between the wetland and the concrete culvert. The stormwater inflow boundary at the right of the model domain, including two meters of the upper surface, is a time-dependent head and salinity boundary, the values of which are scaled to measurements taken at a long-term monitoring site on Boone Creek (see Figure 1 for the location of the Varsity Gym stream gauge relative to the accidental wetland). The seepage boundary on the left, representing the flow of water out of the wetland, is a time-dependent head boundary, with high water levels occurring only when a high stream stage in Boone Creek extends to the accidental wetland. A fluid flux boundary covers most of the top boundary of the model domain and represents the recharge to the mini-aquifer. Five simulated monitoring wells along the base of the model were used to collect output data for this paper. It was at these locations that we performed a semi-quantitative model calibration of our simulations to the field data.

**Figure 6.** Model domain used in the flow and solute transport models. The flow enters on the right side of the figure, represented by black circles, and moves through the wetland to the seepage face in the left bottom corner, indicated by blue circles. The no-flow boundary (red line) represents the concrete culvert.

#### **3. Results**

#### *3.1. Field Data*

As shown in Figure 7, chloride values were highest at the upstream end of the wetland (at the inflow and AW1 and AW2 sites) and decreased with increasing distance into the accidental wetland. The highest salinities occurred at the first three sampling locations, and the lowest salinities left as discharge through the seepage face at the end of the accidental wetland (Table 1). As expected, the salinities decreased as the runoff moved downstream through the wetland, indicating that the wetland mitigated salinity by reducing peak

values and delaying their arrival to the stream. About 47% of the time, the salinity at the seepage face was lower than the salinity at the inflow. This depended on the timing of the sampling because there may have been times when a fresh rainfall/runoff event occurred as a lagged salt event was still working its way through the wetland. Overall, however, as shown in Figure 7 and Table 1, the average and median chloride concentrations were highest at the upstream end of the wetland and were lowest at the downstream end of the wetland, near the seepage face.

**Figure 7.** Boxplots of field chloride concentrations at each site throughout the length of the wetland. The middle red line marks the average chloride concentration. Chloride concentrations decrease towards the seepage face.

**Table 1.** Summary of chloride concentration statistics of the field measurements.

#### *3.2. Synthetic Boundary Conditions*

Salinity data, collected at the inflow point to the AW, were collected with a handheld probe. Salinity data were also collected on a continual basis at the VG stream gauging site, just downstream of the AW. A comparison of these data suggests that salinity values at the AW inflow during melt events were 4.5 times greater than those measured at the VG site (Figure 8). Because we did not have continuous data from the AW inflow, we created a synthetic dataset of accidental wetland inflow salinities by taking the VG salinity data and multiplying those values by 4.5 (Figure 9). In doing so, we are acknowledging that we are not studying the exact conditions arriving at the AW, but are instead using a synthetic time series that is consistent with the observed salinities at the inflow point. These scaled salinities were used in the numerical model as the inflow boundary conditions.

**Figure 8.** Relationship between the inflow concentrations of the Accidental Wetland and Boone Creek sites from hand samples. The black line represents the factor of 4.5 between the Boone Creek and inflow data. The open blue circles are the Boone Creek chloride concentration compared to inflow concentration collected during sampling.

**Figure 9.** Stream stage–time plot (top panel) and concentration–time plot (bottom panel) from the Varsity Gym (VG) site from 1 January 2018 to 1 January 2019. The salinity data (bottom panel) were scaled by 4.5 to match concentrations observed in the Accidental Wetland and were used in the numerical model as inflow boundary conditions.

#### *3.3. Numerical Modeling of Groundwater Flow and Solute Transport*

The model output using the synthetic boundary condition data shows that the accidental wetland at the end of the stormwater culvert reduced peak salinity concentrations and delayed their arrival to Boone Creek. Dispersive transport properties in the wetland reduced peak salinity concentrations as the runoff water flowed as groundwater through the wetland. Here, groundwater refers to the flow of shallow subsurface water in the

wetland substrate. There was no connection between the flow in the wetland and the natural water table due to the concrete culvert. The wetland groundwater velocities varied with the influx of water during storm and meltwater events. The gradients became high and in turn, increased the velocities, thus lowering the lag times; the storage of salt was of longer duration during dryer periods.

Initial simulations were carried out to examine the potential for a multi-year build-up of salt in the accidental wetland. To assess this, we ran a five-year simulation repeating our years' worth of scaled boundary data five times. Figure 10 shows the results of this simulation. The upper panel shows the scaled inflow data and the output seepage data. The lower panel shows the other simulated monitoring wells at a different vertical scale. The salinity concentration decreases through the wetland and approaches a concentration of zero during the summer months. The model output demonstrates that there is little, if any, salt accumulation in the wetland from season to season over the five years. This result is not surprising, given the small size of the wetland and the relatively fast circulation of water through the system.

**Figure 10.** Concentration–time plots of output from the numerical model over five years of simulation time. The top panel shows the differences in salinity between the inflow to the wetland (black line) and the seepage face (dark green). The bottom panel shows the output from the simulated monitoring wells. Salinities decrease with distance travelled through the wetland, with the lowest salinities occurring at the seepage face.

Figure 11 shows simulation output, demonstrating the reduction in peak salinities as dispersion spreads out the salt plume with a transport distance downgradient. Figure 11 focuses on one year of the simulation output from Figure 10. It also demonstrates that there is a notable delay in the arrival of the center of mass of the plume to Boone Creek. The model output in the Figure compares the raw input data, represented by monitoring well 1, to the output data of the seepage face. Peak salinity values at the inflow lag up to approximately 45 days and have a peak salinity reduction of up to 94% by the time the solute reaches the seepage face. These reduced salinities and delayed arrival times, caused by typical solute transport processes in the accidental wetland, prevent the stormwater culvert from being a source of the high-salinity surge that is typical of meltwater events during 'salt season.'

We also performed model sensitivity analyses to understand the effects of hydraulic conductivity, which we had only estimated from Hazen Method calculations and adjusted through model calibration, on modeled salinity concentrations at the seepage face. We varied hydraulic conductivity in the modeled wetland while keeping other parameters, such as boundary conditions and longitudinal and transverse dispersivities, equivalent

to the base case scenarios. As expected, higher hydraulic conductivity values increase the groundwater velocity, which in turn results in a quicker pulse of salinity through the wetland system and higher seepage salinities (Figure 12). Solute transport processes suggest that macrodispersion increases with groundwater velocity; however, the small scale of the aquifer limits the time available for dispersion to take place, and the result is an advection-dominant transport process. For example, at a simulated hydraulic conductivity of 19.8 m/d, peak concentrations rise above 730 mg/L Cl− and have a similar lag time to the base case simulation. Lower hydraulic conductivity values result in lower groundwater velocities, which increase the residence time of the salt in the aquifer. The low-permeability simulation, utilizing a hydraulic conductivity of 5 m/d, lagged the other sensitivity simulations by three months and had a peak concentration of roughly half of the base and high-permeability simulations.

**Figure 11.** Concentration–time plot from the start of the typical salt season to show the mitigation of salt from the inflow to the seepage face. The top panel plots the inflow (black line) and seepage face (dark green line) salinities. The bottom panel plots the salinities of the different monitoring wells from the numerical model. Closer to the seepage face, lower salinities can be observed.

**Figure 12.** Concentration–time plot for the sensitivity analysis of the numerical model of salinities at the seepage face. K = 15 m/d (dark green) was the hydraulic conductivity used in the above simulations. The model is sensitive to changes in hydraulic conductivity.

#### **4. Discussion**

Our study demonstrates that accidental wetlands help remediate salt-laden stormwater runoff by increasing the residence time of the salt and delaying its arrival to a receiving stream, thereby decreasing the peak concentration of the salinity pulse but not the overall quantity of salt that is delivered to the stream.

Our field-sampled data support this conclusion. The median and mean values at the seepage face at the end of the wetland were lower in concentration than those at the upper end of the wetland, where salt-rich runoff entered the wetland during salt events. Our groundwater modeling and solute transport simulations also support this conclusion in a comprehensive way. The accidental wetland slows the delivery of salt-laden runoff, thereby lowering simulated salinities by up to 94 percent and lagging the timing of the peak concentrations by up to 45 days. This type of retention system does not remove the salt, but by slowing down its migration to the adjacent stream, it can lower the chance of chronic chloride violations and greatly reduce the potential for acute salinity spikes. This distinction between acute and chronic is important to stream health in the same way that dosage is relevant to any toxin. As noted previously, in their review of salinity impacts, Hintz and Relyea [17] found significant variability across species and across salinity levels. For example, they found that some species adapt to chronic salinity levels, although acute spikes are lethal. Because road salt will continue to enter streams, finding ways to at least lower the intense spikes of salt should be beneficial. If enough accidental wetlands are encouraged within urban watersheds, rather than removed from them, this could have a positive influence on the overall water quality in streams.

Furthermore, sensitivity analyses with the model showed that hydraulic conductivity values typical of sands are ideal in order for this to take place. This was the primary material in the accidental wetland at our study site and served as the base case for our simulations. It is safe to assume, however, that materials with properties too far outside of this range, such as gravels or clays, would not be effective. The hydraulic conductivity of highly-permeable materials such as gravels would allow rapid migration through the accidental wetland, preventing much of a reduction in peak concentrations or an increase in lag times. Conversely, silts and clays would have so little permeability that surface flow would dominate and would allow the quick delivery of the salt-laden urban runoff to the stream.

We think that there are likely other accidental wetlands similar to the one studied here in other environments; that is, a small scale, concrete substrate wetland forming in an area "designed" for another purpose. The formation of multiple accidental wetlands along a watershed would retain a higher percentage of the road salt, thus decreasing the concentration of the salinity pulses and preventing higher salinities in streams themselves. We can potentially improve the water quality of our urban streams by watching for and then encouraging accidental wetlands in multiple areas throughout a watershed.

#### *4.1. Other Accidental Benefits*

In addition to salinity concerns, accidental wetlands also influence other dynamic hydraulic issues such as stream flashiness and temperature surges in urban environments. Stream flashiness is caused by high percentages of impervious surfaces in the watershed, especially in the lower reaches in the Town of Boone and on the ASU campus, and results in rapid stream-stage increases and reversed-gradient events, which drive road salt into riparian aquifers along the length of urban stream reaches. These processes lead to both chronic and acute stream contamination over much longer time frames than what occurs in natural, unurbanized streams [42]. In our study location, urban, flashy runoff provides a significant source of contamination to the riparian aquifer along Boone Creek, a focus of years of urban hydrology research, and these groundwater–surface water interaction processes are especially important in relation to the chronic and acute chloride contamination of Boone Creek.

When stormwater runoff encounters an accidental wetland, the amount of runoff entering the stream decreases and lags according to the travel time that the wetland requires. Temperature surges occur when heated runoff enters a stream and quickly increases the stream temperature [20]. Depending on the formation of the accidental wetland, any detained water has time to cool, thereby reducing the temperature of hot runoff in temperature surges. If enough of these wetlands form, temperature surges will perhaps be greatly reduced or eliminated within a watershed. Long-term average stream temperatures are also improved by accidental wetlands via natural shading [3]; in fact, several small trees, along with other knee-high vegetation, did grow in the accidental wetland we studied. In addition to reducing impacts from salt, our case suggests that accidental wetlands in urban areas can reduce the effects of urban stream syndrome by mitigating flashiness and heat.

#### *4.2. Management Implications*

As already noted, campus facilities management removed the accidental wetland in 2018. Removing this wetland aligns with historic management approaches as well as continued perceptions about stormwater, flooding, and vegetated waterways. Traditional stormwater management has focused on shunting water as quickly as possible to the nearest stream to avoid flooding. Although new approaches, including constructing wetlands, bioswales, and retention systems, have been implemented, there remains a lack of understanding and misperceptions about these kinds of management techniques. In his overview of the literature, Everett [43] concluded that while there is some public awareness of alternative management approaches, it is not yet mainstream knowledge. People tend to focus on the amenity value or the aesthetics of entities like wetlands or bioswales rather than on their function. In some cases, this creates a negative feedback loop, in which people do not understand what the infrastructure is designed to do—or in this case what it accidentally does—and they therefore do not manage or maintain it appropriately, which leads to further negative views as aesthetics and/or function decline.

In our Boone Creek case, a misperception that the wetland was reducing the channel's ability to handle runoff and thereby increasing flooding risk was the rationale to remove it. This mistaken perception potentially reduced, rather than improved, both flood management and water quality. Because all of the vegetation and sediment was removed, runoff entering the culvert now goes straight into the stream, carrying any salt with it (Figure 2f). Because accidental wetlands arise without intention and are not part of explicit management plans, they can be ephemeral and can be quickly removed. If paying attention to when and where accidental wetlands form became part of overall urban stormwater management and planning, it could offer reduced flooding, improved water quality, and economic benefits, as the wetland costs nothing to establish.

#### **5. Conclusions**

The accidental formation of a wetland in a cold, urban environment inspired us to study the impacts of this system on mitigating road salt. Returning to our research questions for this conceptual study, we conclude that the accidental wetland decreased the peak concentrations of salt and chloride contamination and delayed the arrival of salt to the stream. Simple numerical simulations demonstrated that the saline pulse from urban runoff was delayed by up to 45 days and peak concentrations declined by up to 94%. Although the salt is not removed from the system, its arrival was delayed and the peak concentrations were lowered, ultimately improving water quality and potentially reducing negative consequences for some aquatic species.

We also think that accidental wetlands are likely to serve several functions and potentially offer multiple benefits to stream quality. For example, our model provided a conservative estimate regarding peak salinity reduction because we did not include uptake by plants. Future research assessing the role of vegetation and/or soil in contributing to reducing salinity impacts in accidental wetlands is warranted. Additionally, our results

suggest that accidental wetlands may be effective in lowering runoff-induced thermal pollution due to the delay in the runoff entering the stream. Although this delay does not actually remove the salt, it most likely does reduce the temperature, offering a potentially significant improvement in cold-water streams. Again, more detailed analyses of this potential contribution from accidental wetlands would be valuable. Our work highlights the potential for accidental wetlands to improve water quality issues in urban streams. The only way to decrease salinization issues in cold regions is to discontinue the use of road salt during winter storm events. Because of safety concerns, totally eliminating salt is not practicable and therefore, we need to better understand and explore reactive measures to retain salt and delay its delivery into streams. Our study builds on previous work, suggesting that accidental wetlands cost nothing but offer multiple benefits, which are potentially as valuable as some stream restoration or intentional stormwater retention measures.

**Author Contributions:** Conceptualization, C.M.M. and W.P.A.J., methodology, C.M.M. and W.P.A.J.; software, C.M.M. and W.P.A.J.; validation, C.M.M., W.P.A.J., and K.C.; formal analysis, C.M.M. and W.P.A.J.; investigation, C.M.M. and W.P.A.J.; resources, C.M.M. and W.P.A.J.; data curation, C.M.M. and W.P.A.J.; writing—original draft preparation, C.M.M. and W.P.A.J.; writing—review and editing, C.M.M., W.P.A.J., and K.C.; visualization, C.M.M. and W.P.A.J.; supervision, W.P.A.J. and K.C.; project administration, W.P.A.J. and K.C.; funding acquisition, C.M.M. and W.P.A.J. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data are contained within the article and more detailed data may be requested from the authors.

**Acknowledgments:** The authors thank the Department of Geological and Environmental Sciences, the College of Arts and Science, and the Office of Student Research at Appalachian State University for support of the work. We thank undergraduate researchers Emily Fedders and Andrew Barringer for their help in maintaining the stream monitoring network. We also thank the comments and suggestions of four reviewers of a previous version of this paper.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

# *Article* **Modeling Shallow Urban Groundwater at Regional and Local Scales: A Case Study in Detroit, MI**

**Sadaf Teimoori \*, Brendan F. O'Leary and Carol J. Miller**

Department of Civil and Environmental Engineering, College of Engineering, Wayne State University, 5050 Anthony Wayne Drive, Detroit, MI 48202, USA<EMAIL_ADDRESS>(B.F.O<EMAIL_ADDRESS>(C.J.M.) **\*** Correspondence<EMAIL_ADDRESS>