Relationship between Groundwater Nitrate Concentration and Density of Onsite Wastewater Treatment Systems: Role of Soil Parent Material and Impact on Pollution Risk

Our study was conducted in the town of Charlestown, situated on the south shore of the state of Rhode Island (USA) (Fig. 1). The town has an area of ${107.02\mathrm{km}}^2$ and a population of 8,072 (United States Census 2022). It has three tidally influenced coastal lagoons—locally known as salt ponds—that connect with Block Island Sound through engineered tidal inlets (Fig. 1). A large portion of Charlestown (${41.1\mathrm{km}}^2$) is within the ${82\mathrm{km}}^2$ Salt Ponds watershed, which includes portions of the towns of Narragansett, South Kingstown, Charlestown, and Westerly. The portion of the Salt Ponds watershed within the Town of Charlestown is referred to here as the Coastal watershed.

Charlestown’s coastal zone contains the highest density of septic systems in the town. The delineated Coastal watershed encompasses one-third of Charlestown’s land area but contains 58% of the Town’s developed parcels. Charlestown’s economy is primarily based on its coastal zone where tourism, recreation, and coastal businesses thrive. Recent data indicate that nearly 70% of municipal revenue is generated from within the Coastal watershed (Town of Charlestown, personal communication, 2023).

The geologic setting of Rhode Island’s Salt Ponds watershed has been characterized (Boothroyd et al. 1998; Masterson et al. 2007; Stone and Borns 1986). Briefly, the geology of the area is dominated by surficial, glacially derived materials deposited by the retreat of the Laurentide ice sheet which overlie the irregular, eastward sloping surface of Paleozoic bedrock. The northern boundary of the watershed is delineated by the ridge of the Charlestown moraine, a recessional moraine characterized by hummocky, loose ablation/debris flow till. The moraine is the surface water drainage divide between the Pawcatuck River watershed to the north and the Salt Ponds watershed to the south. South and seaward of the Charlestown moraine, surficial geology throughout the watershed consists of predominantly thick ($>15\mathrm{m}$), stratified glacial fluvial sand and gravel deposited as outwash plains, braided river deltas, and ice contact deposits. Glacial till, typically an unconsolidated, relatively compact mixture of sand, silt, and gravel overlies bedrock and dominates the Quonochontaug headland area and a small headland between Green Hill Pond and Eastern Ninigret Pond.

We analyzed groundwater quality data from 367 individual private potable water wells serving homes within the Coastal watershed (Fig. 1). Well water samples were collected from October 2008 to June 2022 and analyzed for ${\mathrm{NO}}_3^{-}\text{-}\mathrm{N}$, abbreviated here as ${\mathrm{NO}}_3^{-}$. Of these, 209 wells were sampled as part of a voluntary town public health outreach effort during six distinct sampling events from 2013 to 2022. Applications to participate in the program were primarily submitted by owners of dwellings situated within the Coastal watershed in proximity to coastal surface water bodies and in densely developed areas, where health risk of impacted groundwater was expected by personnel from the Office of Wastewater Management to be highest, and in less densely developed areas, where risk was presumed to be lower. Six sampling events were undertaken, in the neighborhoods of Charlestown Beach Road (2013; $n=28$), Tockwotten Cove (2015; $n=33$), Quonochontaug (2016; $n=11$), Green Hill Pond and Eastern Ninigret Pond (2018; $n=41$), Allen’s Cove and Eastern Ninigret Pond (2021; $n=65$), and Green Hill Pond and Central and Eastern Ninigret Pond (2022; $n=31$) (Fig. 1). Samples were collected without regard for the time of year at both seasonal and year-round occupied dwellings. Additional private well samples from different parts of the watershed, referred to as Town Samples ($n=109$), were collected from 2008 to 2021 and laboratory analytical reports were submitted to the town, either as part of community engagement activities through the Charlestown Conservation Commission or to meet regulatory requirements of the Rhode Island Department of Health’s (RIDOH) Rules and Regulations Pertaining to Private Drinking Water Wells (RIDOH 2008). This regulation requires that newly installed private wells be sampled and analyzed for potability prior to the issuance of a Certificate of Occupancy from a municipal building official, and at real estate transactions. A further 49 private well samples were collected in 2010 in the Quonochontaug headland area by University of Rhode Island researchers to assess nutrient loading risk (Donohue 2013).

All wells were actively in use at the time of sampling. Water samples were collected by privately contracted well samplers licensed by the State of Rhode Island from interior taps prior to any domestic water treatment mechanisms. Samples were transported to a RIDOH certified drinking water laboratory under a chain of custody protocol and analyzed for ${\mathrm{NO}}_3^{-}$ using EPA Method SM 4500 ${\mathrm{NO}}_3^{-}\mathrm{F}$ (APHA 1998).

Because they were privately owned, we could not access individual wellheads, and well depth is not known. However, given the location and geology of the study area, we inferred that wells in this study are sited in the shallow unconfined aquifer based on the following:

To approximate the mean thickness of freshwater above sea level in the nearshore wells of the study area—those within 100 m of the coastal feature—we used spatial analysis. Georeferenced seasonal high-water table reported on approvals for 392 septic system design permits by the Rhode Island Department of Environmental Management (RIDEM) sited within the study area were spatially linked in ArcGIS Pro v. 2.5 (Esri Inc., Redlands, CA) to a digital elevation model of ground surface lidar and mean sea level (MSL). We then calculated the mean elevation difference between MSL and the site water table elevation. The resulting mean ($n=392$) thickness of freshwater above sea level within 100 m of the coast in the study area was 1.31 m. Based on Eq. (1), the thickness of the freshwater layer before intersecting a saltwater lens in the near coast portion of the study area is approximately 52 m.

We used spatial analysis to determine OWTS density relative to sampled potable well locations and ${\mathrm{NO}}_3^{-}$ concentration in well water. Sample sites were spatially oriented using ArcGIS Pro. Sampled well locations were mapped by joining addresses of dwelling units where private well data were collected to the Town Assessor’s property parcel data set. The Rhode Island Geographic Information System (RIGIS) (2020) E-911 data set, which identifies locations of individual structures on properties, was used to geolocate each dwelling (RIGIS 2021). The E-911 dataset was also used to approximate OWTS and well locations in the study area, because precise locations for these are not georeferenced. Therefore, each E-911 point represents a dwelling or structure with an OWTS and drinking water well. This is a reasonable assumption given that in the study area, the average parcel size is 0.12 ha, and the locations of wells and OWTS are typically within 15 m from any structure. The representative E-911 points are referred to here as parcel points.

To calculate the number of OWTS in the vicinity of a private well (development density), we first removed undeveloped parcels from the parcel point dataset. For the remaining parcel points, a 121.9-m radius, representing the RIDOH public wellhead protection buffer area, was placed around each point, creating a buffer layer. This radius represents a land area of 4.67 ha. Parcel points were then spatially joined with the buffer layer and the number of parcel points that fell within each buffer area were used as the input to the Point Density tool, with density calculated based on parcel points per hectare. The resulting density layer was used to predict groundwater ${\mathrm{NO}}_3^{-}$ concentration by extracting raster values of well water concentrations and locations to analyze measured ${\mathrm{NO}}_3^{-}$ concentrations at varying development (OWTS) densities.

Soil parent material was classified as till or glacial fluvial deposits, as defined by the USDA Soil Survey Geographic (SSURGO) Database Web Soil Survey (Soil Survey Staff 2023) and delineated using ArcGIS Pro using the 2018 Rhode Island Soil Survey layer (RIGIS 2018) and the Quaternary Geologic Map of the Carolina and Quonochontaug Quadrangles (Boothroyd et al. 2001). Till includes ablation till, loess over ablation till, recessional moraine deposits, and lodgment till, whereas glacial fluvial deposits include eolian sand and/or outwash deposits, fluvial deposits, and loess over fluvial deposits (Soil Survey Staff 2023).

A list of soils series associated with sampling sites is included in Table S1. Saturated hydraulic conductivity of soils with till parent material ranged from 4 to $42\mu \mathrm{m}/\mathrm{s}$, and from 42 to $703\mu \mathrm{m}/\mathrm{s}$ for glacial fluvial parent material. The permeability ratings for till soils are Moderately Rapid to Rapid and glacial fluvial soils are Rapid to Very Rapid. Each well sample location was assigned a designation of till or glacial fluvial deposit when data were analyzed by soil parent material.

We used spatial analysis to geolocate sampling locations by joining well nitrate concentration data to the town’s parcel data in ArcGIS and exporting the resulting coordinates as a new data layer. The output raster was then used to extract density values to the sampled well locations, ${\mathrm{NO}}_3^{-}$ concentration for each well and associated OWTS density, and the mapped soil parent material.

Statistical analyses and linear regression were carried out using SigmaPlot (v. 14.5; Inpixion, Palo Alto, CA). A Mann-Whitney U-test was used to examine differences in well water ${\mathrm{NO}}_3^{-}$ concentration and in OWTS density between soil parent materials. Levene’s test was used to compare the variance of ${\mathrm{NO}}_3^{-}$ concentration among water samples from different soil parent materials. Student’s $t\text{-}\mathrm{test}$ was used to compare slopes and intercepts of linear regressions. A $P$ value less than or equal to 0.05 was used as a threshold for statistical significance for all tests.

We used linear regression to predict ${\mathrm{NO}}_3^{-}$ concentration as a function of OWTS density in both till and glacial fluvial soil parent material aquifers combined and separately. Values of slope and y-intercept were then used to map predicted ${\mathrm{NO}}_3^{-}$ concentrations as a model in ArcGIS, with the OWTS density raster ($x$) input to the line equation using a GIS raster calculator tool, producing a raster layer of predicted ${\mathrm{NO}}_3^{-}$ concentrations.

Cumulative probability distribution was estimated using the equation:where P = probability of groundwater ${\mathrm{NO}}_3^{-}$ concentration occurring in a sample density class; n = number of samples in each density class; and X = nitrate concentration detected at individual OWTS densities in each density class.

The cumulative probability was estimated for ${\mathrm{NO}}_3^{-}$ concentrations in wells from 10 different OWTS density classes ($0.4<$ to $>3.6\mathrm{OWTS}/\mathrm{ha}$) to determine the proportion of samples with ${\mathrm{NO}}_3^{-}$ values in the range of individual risk categories based on criteria established in the Guide for Completing Source Water Assessments in Rhode Island (Table 1) (RIDOH et al. 2010; URI and RIDOH 2010). In terms of human health risk, the USEPA threshold for nitrate-N in drinking water is $5\mathrm{mg}/\mathrm{L}$, half the MCL of $10\mathrm{mg}/\mathrm{L}$.

Estimates of yearly N load from individual OWTS within the Coastal watershed were based on the following assumptions:

The concentration of ${\mathrm{NO}}_3^{-}$ in wells from combined glacial fluvial and till soil parent material data sets ($n=367$) ranged from 0.0 to $15.0\mathrm{mg}/\mathrm{L}$, with a median value of $2.9\mathrm{mg}/\mathrm{L}$ (Fig. 2). The density of OWTS ($\text{systems}/\mathrm{ha}$) from combined glacial fluvial and till soil parent material data sets ($n=367$) ranged from 0.01 to 4.81, with a median value of 1.2 (Fig. 2). Density values were not normally distributed.

For wells in glacial fluvial parent material ($n=207$), ${\mathrm{NO}}_3^{-}$ concentration ranged from a 0.0 to $9.7\mathrm{mg}/\mathrm{L}$, with a median value of $2.7\mathrm{mg}/\mathrm{L}$ (Fig. 2). The samples from wells sited in till parent material ($n=160$) had ${\mathrm{NO}}_3^{-}$ concentrations from 0.0 to $15.0\mathrm{mg}/\mathrm{L}$, with a median value of $3.7\mathrm{mg}/\mathrm{L}$. Values of ${\mathrm{NO}}_3^{-}$ concentration for wells from both types of soil parent material were not normally distributed. The difference in ${\mathrm{NO}}_3^{-}$ concentration between soil parent materials was statistically significant.

The density of OWTS in glacial fluvial parent material ($n=207$) ranged from 0.04 to 4.8, with a median of 1.0 (Fig. 2). For till parent material ($n=160$), density ranged from 0.012 to 3.81, with a median of 1.3. Values of OWTS density were not normally distributed for either parent material, and there were no statistically significant differences in density between parent materials.

There was a statistically significant linear relationship between OWTS density and groundwater ${\mathrm{NO}}_3^{-}$ concentration for wells located in glacial fluvial ($r^2=0.167$) and in till ($r^2=0.189$) parent material (Fig. 3). Furthermore, the slope value ($\mathrm{mgN}$ · $\mathrm{ha}/\mathrm{L}$ · OWTS) for till (0.99) was significantly higher than for glacial fluvial (0.63) parent material. In contrast, the y-intercept values of the regression lines, which correspond to nitrate concentration when OWTS density is $0\mathrm{OWTS}/\mathrm{ha}$, were $2.3\mathrm{mg}/\mathrm{L}$ for till and $1.9\mathrm{mg}/\mathrm{L}$ for glacial fluvial parent material and were not significantly different.

The spatial distribution of groundwater ${\mathrm{NO}}_3^{-}$ concentrations in soils with till and glacial fluvial parent material was modeled using the linear relationship between ${\mathrm{NO}}_3^{-}$ concentration and OWTS density. For comparison, measured ${\mathrm{NO}}_3^{-}$ values were overlain on the model as points with symbols representing different concentration ranges. The general concordance of modeled and measured values suggests that the linear relationship between OWTS density and nitrate levels can reasonably represent the spatial distribution of nitrate values. The model predicts groundwater ${\mathrm{NO}}_3^{-}$ concentrations of $5\mathrm{mg}/\mathrm{L}$, corresponding to Extreme pollution risk (Table 1) in the Eastern Ninigret Pond and Green Hill Pond areas (Fig. 4) and in the Quonochontaug headland (Fig. 4), both of which have a $\text{density}>2.0\mathrm{OWTS}/\mathrm{ha}$.

Analysis of the data based on cumulative probability distribution shows that fewer than 5% of the 207 wells located in glacial fluvial parent material had a probability of having ${\mathrm{NO}}_3^{-}$ concentrations of $5\mathrm{mg}/\mathrm{L}$ or more, which corresponds to Extreme pollution risk (Table 1), in areas where OWTS density was $<2.0\mathrm{OWTS}/\mathrm{ha}$ (Fig. 5). In contrast, at densities above $2.0\mathrm{OWTS}/\mathrm{ha}$ nearly 10% of the sites had a probability of ${\mathrm{NO}}_3\text{-}\mathrm{N}$ concentrations $\ge 5.0\mathrm{mg}/\mathrm{L}$.

Nitrate concentration values from wells located in till parent material ($n=160$) were more variable than those in glacial fluvial material and had a higher proportion of values above the Extreme risk threshold of $5\mathrm{mg}/\mathrm{L}$ (Fig. 5). At low densities ($0–0.8\mathrm{OWTS}/\mathrm{ha}$), 10% or fewer wells had levels of ${\mathrm{NO}}_3^{-}$ concentrations above $5\mathrm{mg}/\mathrm{L}$. However, at densities of $0.8–1.2\mathrm{OWTS}/\mathrm{ha}$, 50% of wells had a ${\mathrm{NO}}_3^{-}$ concentration above $5\mathrm{mg}/\mathrm{L}$, and 8% exceeded the EPA drinking water standard of $10\mathrm{mg}/\mathrm{L}$. This trend was observed with increasing OWTS density across all sites sampled. These results indicate that, in areas of soils with till parent material, Extreme pollution risk is expected at OWTS densities of $0.8–1.2\mathrm{OWTS}/\mathrm{ha}$. (Fig. 5).

When both soil parent materials and all OWTS density classes ($0-3.6\mathrm{OWTS}/\mathrm{ha}$) are considered, more than 75% of all sites had ${\mathrm{NO}}_3^{-}$ levels above the High risk threshold of $2.0\mathrm{mg}/\mathrm{L}$.