Quantifying the land and population risk of sewage spills overland using a fine-scale, DEM-based GIS model

Model design

The spill model is composed of three components that represent where sewage from a determined spill location may flow across land and accumulate in areas using an elevation raster; see Fig. 2. The first component (purple) is the area surrounding the geocoded spill address location; we call this an accumulation zone. The second component is the flow traceline (green). The third component is an accumulation zone surrounding the end point of the traceline (orange).

Accumulation zones

The accumulation zones (components 1 and 3) are polygons surrounding the start and end points. The start point is the geocoded position of where the spill is indicated as starting. The end point is the point where the traceline (described below) terminates. We create accumulation zones to represent possible pooling areas around these start/end points. The process for creating both these zones is the same for start and end point.

First, we calculate the potential accumulation zone using a dynamically sized buffer based on the spill quantity, see Eq. (1). The size, $s_i$, is dependent on quantity, $q_i$ of the $i^{th}$ spill, Q is the set of all quantities. The formula is a proportional equation to calculate the buffer in relation to the max quantity $(max(Q))$ of all the spills in the dataset. We chose a maximum buffer size of 50 m and the size of the specific spill’s quantity ( $q_i$) is calculated proportionally. We chose to take the $log$ of the spill quantity in order to mitigate outliers in spill quantities that were found in our dataset. An illustration of this calculated possible buffer size can be seen in Fig. 3, as the dark red circle. This buffer size can change depending on the use-case of the SFER model, and the outlier normalization is not necessary and the potential buffer size can be increased or decreased depending on the perceptions of the type of spill (e.g., necessary safety factor) by a decision-maker.

(1)

$$s_i=\frac{log(q_i)}{\frac{log(max(Q))}{50}}$$

After the potential accumulation buffer is created, the areas of elevation that are below the start/end point are aggregated. First, the elevation value of either the start or end point is compared to the elevation value of all the other points in the potential accumulation zone buffer. In Fig. 3, the focus is on the end point accumulation zone, but the process would be the same for the start point. The dark red point in the middle of the figure is the end point and the number 314 is the elevation value in meters above sea level.

To more easily compare the elevation of the end point to surrounding cells in the elevation raster, the elevation raster is converted into points using a raster to point GIS function. A buffer of 4.5 m was created around these points to approximate the area that would be covered by the original 1/3 arc second (approximately 5-m resolution) raster. Simply using the point location instead of buffered points would result in an underestimation of the accumulation zone. See Appendix A.III for further information on building this layer. In Fig. 3, the elevation points are depicted as semi-translucent white points and their corresponding labels are the elevation raster cells values. The elevation points (white points) within the potential accumulation buffer region (dark red polygon) are compared with the start/end point (red point) and aggregated using the PostGIS (2023a) ST_Collect and then turned into a polygon (orange geometry) using the PostGIS (2023b) ST_ConcaveHull functions. In Fig. 3, the orange polygon encapsulates the elevation values that are less than 314 m and represents the accumulation zone.

Assessment of risk

Using the SFER model, exposure risk to populations and to land can be assessed. Two metrics within the domain of sewage spills that could be applicable for decision-makers to assess were developed. First, the risk that historical sewage spills have posed to the environment, and second the public health risk of human community exposure can be examined.

To evaluate public health risk, a metric utilizing an estimated population exposed by the spills within each census tract was generated. Census tracts are regions as defined by the U.S. Census Bureau (U.S. Census Bureau, 2021a). Census tracts were chosen because of their ability to be joined with other demographic information like age or income level. The estimated ratio of people exposed to a spill to the people living in the census tract indicates which census tracts are at higher risk of exposure to its inhabitants. The population estimated by the SFER model is calculated using the LandScan USA Night 90-m raster (Weber et al., 2021). The LandScan raster was converted into points and a 35 m buffer was added; see Appendix A.III for description of how this layer is built and Fig. 4 for a visual representation of the calculation of this ratio. In this figure, the dotted circles with numbers are the LandScan points with the value of the number of persons estimated to live within the corresponding raster cell. The orange polygons are the SFER models. Where these spill models intersect the LandScan points, these points have a bolded outline. The intersecting points with the SFER models are summed, and that value represents the estimated persons exposed to a spill. If there are spills that occur in the same location, the SFER model is first unioned, so the population values impacted are not double counted.

The population risk metric can be seen in Eq. (2). The union of $n$ spill geometries $A_1,A_2,A_3,...A_n$ that intersect with a census tract, $ct$, results in a geometry that is the union of the spills within that $ct$. This unioned geometry represents the area that may be been exposed to sewage within the census tract, and thus with another intersection of Landscan Raster, L, it is possible to sum the population that may be been exposed to the spill. This metric can be applied across all census tracts in the region being analyzed.

(2)

The second metric to calculate risk is a ratio that indicates risk of the land area in a subwatershed region that may be contaminated by a sewage spill originating from a previously reported spill. These subwatershed regions are from the National Hydrology Dataset Plus (U.S. Geological Survey, 2022). These are regions that have been designated as flowing into a single stream segment. In Fig. 5, there are two watershed subregions (blue polygons designated with a and b) and the spills (orange polygons) that occurred within the regions. The total area of subwatershed a is 1,398,358 m², and the one spill that occurred within it is 3,708 m². For subwatershed b, the area is 630,274 m², and of the six spills within this region (five occurred at the same location with quantity differences), their unioned area is 12,186 m².

This land risk metric is shown in Eq. (3). The union of $n$ spill geometries $A_1,A_2,A_3,...A_n$ that intersect with the subwatershed region, $sw$, results in a geometry that is the union of the spills within that $sw$. The resulting geometry represents the area of land that may have been exposed to the spill within the subwatershed region. The land area of the modeled spill footprint, over the entire subwatershed region area represents land risk.

(3)

For both of these metrics, when spill footprints overlap, the union function described in the numerator is used to ensure that neither the population counts (Eq. (2)) nor the land areas (Eq. (3)) of the spill model are counted more than once. This was chosen to highlight areas where the risk is high irrespective of the number of times a spill has occurred at a given location. It is entirely possible to instead use the sum of the spill geometries’ indicating chronic spills, which is a risk on which some decision-makers may wish to focus.

Application of SFER model

To illustrate the utility of the SFER model, it was applied to sewage spill data reported by the Environmental Protection Agency (EPA) of the State of Georgia, USA for the 14-month period between December 2020 and January 2023 (Georgia Environmental Protection Agency, 2023). Georgia reported 1,909 spills, releasing a total of approximately 287 million liters of untreated sewage during this time window. Among the data reported for each spill, the address of spill release and the estimated volume of untreated sewage were provided.

To use these within our model, each spill location was first geocoded using the ArcGIS Online Geocoder service (see Appendix A.II for details). Of the total spills, SFER modeling generated spill footprints for 1,477 of these spills (about 78%). The 432 dropped records were due to geocoding errors and inconclusive quantity of spill amounts.

Of the 1,969 Census Tracts in Georgia, 768 experienced sewage spills in the 14-month window. A total population of 3,099,925 people are reported to live with the affected census tracts. Fine-scale estimates using the LandScan USA Night 90-m rasters of potential population exposed within these census tracts is estimated at 32,760. Of the total subwatershed areas containing a spill zone the overall subwatershed area sums to 161,063 hectares. The total area of the SFER models within this time window was 1,104 hectares. On average, based on the SFER model, 0.53% of the area of an affected subwatershed is impacted by a spill zone. In the maps and figures that follow, we only display census tracts and subwatersheds that have at least one spill incidence within the region.

We have created an interactive map with the layers for the land and population risk across all of Georgia, which can be accessed at elmcdaniel.com/GA_USA_2023_SFERMap. Also, the layers are readily accessible for download within our project repository (https://osf.io/hs6g9/). This map is composed of three layers: subwatershed regions ( $lRisk$), census tracts ( $pRisk$) for night, and census tracts ( $pRisk$) for day. The interactive map is classified using the quantile method with five classes for each of the layers. The lighter colors indicate lower values of the $pRisk$ or $lRisk$ metrics. There are few regions outside of urban areas that have recorded spills. This could be representative of sewage infrastructure differences between urban and non-urban environments or lack of reporting of these types of sewage spills. The difference in perceived impact of sewage spills between the census tracts and subwatershed region map is a result of the relative area covered by spatial units representing each domain. Census tracts must contain a certain population count, and thus in areas with lower population density, census tracts are larger. In areas with less population density, utilizing census tracts as the $pRisk$ region may be less valuable, instead, it would be more beneficial to use finer scale administrative boundaries such as city/place.

Maps of population and land exposure may be used to inform different aspects of public policy making. Decision-makers may use the SFER models to prioritize areas within their jurisdictions where interventions may be needed. We focused on Atlanta as our jurisdiction of interest to illustrate the use of our methodology. We did this for two reasons: (1) sewage spills are of particular concern in the City of Atlanta in Georgia, which historically has been in non-compliance with requirements on regulating Combined Sewer Overflows (CSOs) and Seperate Sewer Overflows (SSOs) (U.S. District Court, Northern District of Georgia, Atlanta Division, 1998, 1999); and (2) the City of Atlanta has a department called Atlanta Watershed Management which is in charge of the sewage system for the city, and thus this could be a planning process they would undertake.

Quantification of risk

Decision-makers that are interested in assessing the risk of exposure to populations may find that it is beneficial to utilize census tracts, estimated population exposed via LandScan data, and our population risk metric as we have described. See Fig. 6A of the population exposure risk score in a choloropleth map that is classified using the quantile method (equal number of tracts per color). Central downtown Atlanta does not have high risk, which can could be attributed to the use of night time population data (in Appendix B, Fig. B1 is the corresponding with the Daytime Raster utilized). Higher risk in Fig. 6A is delegated to the outer rings of the city, where there are more residences and therefore more nighttime population. There are three clusters of adjacent census tracts in the northwest, northeast, and southeast that are within the highest quintile of risk, which may represent areas that decision-makers may choose to prioritize. The high exposure risk that is quantified with the $pRisk$ metric in this case may indicate systemic problems such as aged infrastructure.

Further, it could be beneficial to prioritize regions that contain a higher number of persons vulnerable to gastrointestinal diseases, like persons under the age of 5. This can be done by joining the census tracts to data available in either the Decennial Census or American Community Surveys (ACS), another benefit for utilizing this regional boundary. Figure 6B is a map that takes the top quintile of Fig. 6A, and overlays the counts of persons under 5 years old from ACS 2021 (U.S. Census Bureau, 2021b). In this case, decision-makers could target the western and southern clusters of the city as they have the highest counts of populations under 5 years old. Other demographic parameters of interest in ACS could be used to assess other indicators of vulnerability such as estimated counts of persons above 65 years old, under the poverty line, or medical insurance coverage.

For those interested in assessing environmental risk, they may wish to combine information from the $lRisk$ map with other spatial information such as the percentage of impervious surface area in each subwatershed. In Fig. 6D, the mean imperviousness score was calculated for the top quintile of the subwatersheds in Fig. 6B using the Zonal Statistics as Table tool (Esri, 2023d) from the NLCD 2019 Percent Developed Imperviousness raster (U.S. Geological Survey, 2019). While subwatersheds in the central core of Atlanta contain urban landuse, most are also located in larger subwatersheds. As expected, the highest risk areas are primarily smaller watershed regions with higher percentages of impervious surfaces, and are scattered throughout the city. Impervious surfaces are often used as a surrogate measure of the amount of urbanization, and Melchiorri et al. (2018) have shown that rapidly urbanizing landscapes can adversely affect the natural resources which human populations rely upon. Impervious surfaces can concentrate and funnel overland flow, adversely affecting numerous environmental systems (e.g., nutrient cycling, stream temperature, flood attenuation), especially contaminant filtration and water purification (National Research Council, 2002). Multiple studies have shown that strategies applied to rapidly urbanizing areas can be used protect high quality environmental resources (Atkinson, Hunter & English, 2010) or restore degraded environmental resources (Atkinson & Lake, 2020).

Prioritizing locations

Public health planners and watershed managers working together can simultaneously explore both demographic data that inform public health decisions and watershed management decisions to concentrate limited resources in the most efficient combination. One method to do this is illustrated in the top figure in Fig. 7. This map of Atlanta contains pins indicating where the top quintiles for population risk (Fig. 6A) and land risk (Fig. 6C) intersect. This could be a method for decision-makers to prioritize broad locations of concern. In particular, there are a collection of five pins in the northeast portion of Atlanta that could be prioritized. Additionally, these locations may be qualitatively examined for place context using web mapping tools (e.g., Google Maps, Earth, or Streetview) and mobility data (e.g., SafeGraph). The bottom of Fig. 7 displays the street view (Google Maps, 2022) of a park located adjacent to a community farm where multiple spills were found to have occurred totaling an estimated 662,123 liters of sewage spilled. If decision-makers are aware of certain high-risk areas, they can qualitatively assess these areas in relation to the SFER model.

Limitations

The SFER model is sensitive to the accuracy and resolution of the input datasets used. Inaccuracies in spill locations and elevation data will impact the footprints generated and the derived measures of environmental or population risk. In our study, spill locations were determined by geocoding the reported addresses of where spills were recorded. The exact latitude and longitude of the source of spills as well as the use of higher-resolution DEM can result in a more robust footprint. Further, the use of aggregated population data at the census tract level may introduce spatial biases in population-based assessments of risk. For example, census tracts classified as high risk due to the presence of selected population vulnerabilities may be inaccurate if those populations are clustered in areas with no spills. While high-resolution data on population estimates like LandScan USA may be used to identify overall population exposure at fine-spatial scales, information on specific population characteristics at that scale are not available.

Estimates of land risk using SFER represent relatively simple overland flow models as a function of spill volume and surface elevation. Other environmental parameters may either exacerbate or attenuate the actual environmental risk. For example, if the spill flows through a high-quality wetland, the contaminates may be degraded much more quickly as compared to the spill accumulating along roadways. Likewise, a spill that accumulates in a recreational area such as a community baseball complex or a dog-park may expose more people than census or LandScan USA population estimates would indicate. Including accurate landuse/landcover data into the type of proactive planning that SFER suggests would improve the plans that are developed. The SFER model is limited to use for dry land spills and treats any input location of spills as dry land spills. It also does not account for differences in the surfaces. The incorporation of modeling spills using SFER Model with hydraulic components was out of scope for this model as the focus is on dry land spills.