A Rapid Assessment Method to Identify Potential Groundwater Flooding Hotspots as Sea Levels Rise in Coastal Cities

: Sea level rise (SLR) will cause shallow unconfined coastal aquifers to rise. Rising groundwater can emerge as surface flooding and impact buried infrastructure, soil behavior, human health, and nearshore ecosystems. Higher groundwater can also reduce infiltration rates for stormwater, adding to surface flooding problems. Levees and seawalls may not prevent these impacts. Pumping may accelerate land subsidence rates, thereby exacerbating flooding problems associated with SLR. Public agencies at all jurisdiction levels will need information regarding where groundwater impacts are likely to occur for development and infrastructure planning, as extreme precipitation events combine with SLR to drive more frequent flooding. We used empirical depth-to-water data and a digital elevation model of the San Francisco Bay Area to construct an interpolated surface of estimated minimum depth-to-water for 489 square kilometers along the San Francisco Bay shoreline. This rapid assessment approach identified key locations where more rigorous data collection and dynamic modeling is needed to identify risks and prevent impacts to health, buildings, and infrastructure, and develop adaptation strategies for SLR.


Introduction
Sea levels are rising over most of the world's coastlines, and the rate of relative sea level rise (SLR) is projected to accelerate [1]. One of the impacts of SLR will be a rising water table in shallow, unconfined coastal aquifers [2]. Coastal regions that are currently above sea level do not typically manage this shallow coastal groundwater as a resource because it is often contaminated by agricultural chemicals, industry, or urban surface runoff. Maps of depth to this shallow groundwater are rare, although soil contamination is sometimes monitored locally using well samples. As a result, many coastal regions are unprepared to manage the potential impacts of a rising water table.
Shallow saline aquifers and unconfined freshwater aquifers with a direct saltwater interface (i.e., freshwater floating atop higher-density seawater) are affected by tidal fluctuation. These aquifers rise and fall with the tides, and the effects decrease exponentially farther inland [3][4][5]. In the zone where aquifers are affected by tidal flux, they are also affected by SLR. In "flux-controlled" systems, where the rate of groundwater discharge is constant as the sea level rises, SLR causes landward migration of the saltwater toe, otherwise known as saltwater intrusion [6,7]. This saltwater intrusion causes a lift in the level of the overlying freshwater [8]. Therefore, SLR causes an increase in the height as well as the salinity of the water table [4][5][6][8][9][10]. This eventually results in the emergence of groundwater as surface flooding, and also increases surface discharges of streams supplied by groundwater [2]. Before emergence occurs, rising groundwater infiltrates sewer pipes, causing a loss of sewage flow capacity. It also conveys pollutants to nearshore aquatic ecosystems, floods basements, causes heaving of foundations and underground structures, remobilizes soil contaminants, and increases the risk of liquefaction in seismic regions. In coastal areas constructed on former wetland soils, lowering the elevation of groundwater by pumping can accelerate subsidence [11].
Planning for rising and emergent groundwater at a regional scale requires mapping methods that are suitable for large geographic regions and heterogenous subsurface conditions produced by urbanization. Large empirical depth-to-groundwater datasets often exist for urban areas where leaky underground fuel or chemical storage tanks are regulated and monitored. Maps interpolated from these large empirical datasets can be used to support prioritization of limited flood adaptation resources by identifying "hot spots" in the distribution of risk. Similar methods were used to identify gaps in protection of species over large geographic datasets for conservation planning [12], to identify local extremes of the Urban Heat Island effect [13], and for risk assessment [14]. These types of rapid assessments represent conditions over large geographic areas, allowing public resources to be used strategically to gather new data and develop process-based models where future problems are most likely to occur [15]. In the case of groundwater data, the use of an empirical method that interpolates a surface from a dataset of present-day conditions can provide modelers with initial insights into the complex interactions between heterogeneous soil and infrastructure conditions and groundwater elevation and flow in an urban area [4,5]. While empirical methods that use interpolation do not model the dynamics of groundwater, they can be used over very large geographic areas and can reveal localized effects, such as cracked pipe joints, private water pumps, compacted road beds, or faulted sediment. These local anomalies might not appear in a process-based equilibrium model of the water table, yet could present a significant problem for local adaptation. An extensive well dataset allows managers in a coastal region to use interpolation to anticipate the flooding impacts of groundwater as the sea level rises.

Bay Area Sea Level Rise and Groundwater
Sea level has risen 1.1 mm/year on average at the Golden Gate since the historical record began in 1855 [16]. As the rate of rise increases over the next century, flooding is expected to affect a wide range of assets in the San Francisco Bay Area (Bay Area), including low-lying urban areas, two international airports, wetland ecosystems, and essential infrastructure [17,18]. If nothing is done to intervene, ecosystem shifts are likely to occur in San Francisco Bay (the Bay) that could cause all existing inter-tidal marshes to become mudflats with 1.24 m of sea level rise [19]. Extensive impacts on urban development are also anticipated; low-lying coastal homes, businesses, and infrastructure will be in danger of regular flooding as the sea level rises. Planners now have access to maps that predict direct seawater inundation at different sea levels [18,20]. While direct coastal inundation will have a considerable impact in the Bay Area, inundation due to rising groundwater levels is an SLR-induced threat that has received far less attention in coastal adaptation discussions and has been missing from maps of coastal flood risk. Yet, the importance of studying coastal groundwater dynamics in the context of urban and coastal zone management is recognized as an urgent need. In one study, twice as much urban land appeared to be at risk of flooding when rising groundwater was included in coastal flooding predictions [5].
New data on regional rates of land subsidence in the Bay Area indicated that additional flooding could be expected as a result of elevation changes [21]. Subsidence rates of more than 5 mm per year (and up to 10 mm per year in one location, between 2007 and 2010) were identified in areas where urban fill was placed over thick Bay mud deposits. A lower land surface will expose many new areas to flooding from seasonally high groundwater levels, as well as seawater.
The rate of rise in the groundwater surface due to SLR is affected by a number of factors, including tidal forcing, aquifer geology, coastline change, shore slope, surface permeability, precipitation, and groundwater pumping [4][5][6]. Previous studies indicated that the relationship between SLR and the elevation of the water table is unlikely to be consistently linear, especially near streams [2,10,22]. However, several studies used a linear relationship to approximate the effect of SLR on groundwater levels in flux-controlled urban aquifers [4,5]. Since this method is only applicable in zones where the sea level and tidal fluctuations have an influence on the aquifer, one kilometer was used to represent the limit of that zone in these studies [5].
Like many other coastal regions, the Bay Area did not previously have a depth-to-water map for its shallow coastal aquifers, although some local studies were available. In this paper, we present a rapid assessment method to provide this critical planning tool using empirical data.

Study Area
The groundwater basins of the Northern San Francisco Bay Area are largely within valleys formed on alluvial fans. While the deep aquifers are disconnected, the shallow coastal aquifer is continuous in the alluvial deposits [23]. The shallow aquifer in the large Santa Clara Valley basin (which contains five sub-basins) is also unconfined. Due to groundwater withdrawals, saltwater intrusion has historically occurred in this basin. A reduction in pumping and concerted recharge efforts slowed the progression of saltwater inland [23,24]. For the purposes of our analysis, we assume that the shallow coastal aquifer is unconfined and has a direct connection to the Bay. This assumption is reasonable because most of the Bay-front within one kilometer of the Bay is composed of alluvial material and urban fill.

Monitoring Well Data
Our methods follow similar studies for Honolulu, HI [5] and three locations in coastal California (excluding the San Francisco Bay Area) [4] but cover a much larger geographic area than either of these studies. We used a dataset of groundwater monitoring well measurements that contained values for depth to the water table and covered portions of all nine Bay Area counties [25]. The data points were concentrated in heavily developed areas, with fewer wells in the northern Bay Area (Figure 1). We included wells within 1.6 km of the Bay edge to ensure continuity in our interpolated results, although the only results shown were within 1 km of the Bay. We used the San Francisco Estuary Institute's Bay Area Aquatic Resource Inventory delineation of open water and tidal wetland to define the Bay's edge [26].
We selected the minimum depth-to-water value for each well during the years 1996-2016. This represented the seasonal high water table during wetter years, allowing us to estimate the highest elevation of the water table. Where this maximum groundwater elevation occurred, remobilized pollutants and reduced sewer pipe capacity may have been present in an unusually wet year or during an exceptionally high tide event. We selected the minimum depth-to-water value for each well during the years 1996-2016. This represented the seasonal high water table during wetter years, allowing us to estimate the highest elevation of the water table. Where this maximum groundwater elevation occurred, remobilized pollutants and reduced sewer pipe capacity may have been present in an unusually wet year or during an exceptionally high tide event.

Quality Control for Monitoring Well Data
Many urban wells in this dataset were in close proximity to wells used to measure the water quality in deeper aquifers. Therefore, we deleted wells with a minimum depth-to-water value greater than two standard deviations above the mean (i.e., deeper than 6.46 m). By deleting these wells, we had higher confidence that our interpolated surface represented the shallow coastal aquifer that was responsive to SLR and rainfall. A summary of the depth-to-water data is shown in Table 1.

Quality Control for Monitoring Well Data
Many urban wells in this dataset were in close proximity to wells used to measure the water quality in deeper aquifers. Therefore, we deleted wells with a minimum depth-to-water value greater than two standard deviations above the mean (i.e., deeper than 6.46 m). By deleting these wells, we had higher confidence that our interpolated surface represented the shallow coastal aquifer that was responsive to SLR and rainfall. A summary of the depth-to-water data is shown in Table 1.

Tidal Data
We also included tidal data points from a dataset produced for the Federal Emergency Management Agency and regional agencies [27]. To smooth the interpolated surface toward the Bay, we included 603 mean tide line points, and added 0.3 m to the elevation to reflect the expectation that freshwater usually lies above the mean tide line [28]. Since the tidal water levels varied substantially along the shoreline as a result of the hydrodynamics in San Francisco Bay, tide gauge locations alone were insufficient. The tidal dataset we used was calibrated to National Oceanographic and Atmospheric Administration tide gauges and provided extensive spatial coverage along the Bay shore.

Analysis
For each well point in the study area, we extracted a ground elevation from the United States Geological Survey (USGS) Coastal and Marine Geology Program 2 m digital elevation model (DEM). We then calculated the maximum water table elevation at each well point by subtracting the minimum depth-to-water value from this ground elevation, following the method used by Hoover et al. [4]. Next, we applied a set of interpolation algorithms to the groundwater elevations and tidal data points. A flowchart describing our analysis methods is shown in Figure 2.

Tidal Data
We also included tidal data points from a dataset produced for the Federal Emergency Management Agency and regional agencies [27]. To smooth the interpolated surface toward the Bay, we included 603 mean tide line points, and added 0.3 m to the elevation to reflect the expectation that freshwater usually lies above the mean tide line [28]. Since the tidal water levels varied substantially along the shoreline as a result of the hydrodynamics in San Francisco Bay, tide gauge locations alone were insufficient. The tidal dataset we used was calibrated to National Oceanographic and Atmospheric Administration tide gauges and provided extensive spatial coverage along the Bay shore.

Analysis
For each well point in the study area, we extracted a ground elevation from the United States Geological Survey (USGS) Coastal and Marine Geology Program 2 m digital elevation model (DEM). We then calculated the maximum water table elevation at each well point by subtracting the minimum depth-to-water value from this ground elevation, following the method used by Hoover et al. [4]. Next, we applied a set of interpolation algorithms to the groundwater elevations and tidal data points. A flowchart describing our analysis methods is shown in Figure 2.  [29][30][31][32]. The dataset used here was not well-suited to kriging because it did not fulfill the assumption of stationarity necessary for this method. Data variance was not constant across the study area and could not be explained by directional trends. Given these limitations for geostatistical methods, we only compared deterministic interpolation methods.  [29][30][31][32]. The dataset used here was not well-suited to kriging because it did not fulfill the assumption of stationarity necessary for this method. Data variance was not constant across the study area and could not be explained by directional trends. Given these limitations for geostatistical methods, we only compared deterministic interpolation methods.
To maximize the model accuracy, we compared a variety of methods to determine which was most successful at minimizing the prediction error. We compared the root-mean-square error (RMSE) of predicted values from each model using the cross validation function of the ArcGIS Geostatistical Analyst toolbox ( Table 2). Values of RMSE closer to zero indicated a more accurate model, and values of mean error (ME) closer to zero indicated a less biased model. For each interpolation technique, we chose the input parameters (e.g., power, number of neighbors included) that were most successful at minimizing RMSE, rather than those that produced the smoothest output. The method that minimized the RMSE most successfully was the multiquadric radial basis function. A scatterplot of the actual water table elevations (elevation from DEM minus minimum measured depth-to-water) compared to the predicted water table elevations (output from the multiquadric radial basis function interpolation) is shown in Figure 3. Next, we subtracted the interpolated water table surface from the ground surface DEM to produce a depth-to-water map. We excluded areas greater than 1 km from the nearest well due to the increased uncertainty introduced by the lack of well or tidal data points in these areas.
To maximize the model accuracy, we compared a variety of methods to determine which was most successful at minimizing the prediction error. We compared the root-mean-square error (RMSE) of predicted values from each model using the cross validation function of the ArcGIS Geostatistical Analyst toolbox ( Table 2). Values of RMSE closer to zero indicated a more accurate model, and values of mean error (ME) closer to zero indicated a less biased model. For each interpolation technique, we chose the input parameters (e.g., power, number of neighbors included) that were most successful at minimizing RMSE, rather than those that produced the smoothest output. The method that minimized the RMSE most successfully was the multiquadric radial basis function. A scatterplot of the actual water table elevations (elevation from DEM minus minimum measured depth-to-water) compared to the predicted water table elevations (output from the multiquadric radial basis function interpolation) is shown in Figure 3. Next, we subtracted the interpolated water table surface from the ground surface DEM to produce a depth-to-water map. We excluded areas greater than 1 km from the nearest well due to the increased uncertainty introduced by the lack of well or tidal data points in these areas. The "actual" water table elevation was ground elevation minus the minimum measured depth-to-water. The predicted water table elevation was the value extracted from the raster output of the radial basis function interpolation at the well point location. The plot includes points within the final study area only (i.e., within 1 km of the Bay edge). Figure 4 shows the results of our depth-to-water modeling for the coastal Bay Area. A geospatial data file can be downloaded at https://datadryad.org/stash/dataset/doi:10.6078/D1W01Q We found that a shallow groundwater condition exists in many developed areas in the North Bay including Fairfield, Novato, San Rafael, and Petaluma, although fewer well data points were available for the North Bay in general. Many cities in the East Bay also had shallow groundwater along the Bay-front, placing major infrastructure (such as Interstate highways 580 and 880) at risk. Exposure to potential groundwater flooding was perhaps most severe in the Silicon Valley area, where the minimum depth-to-water was already less than one meter in large areas of Mountain View, Redwood City, and San Mateo. Figure 5 shows subset maps of groundwater conditions in selected highly urbanized areas. Figure 3. Actual water table elevation compared to predicted elevation from the interpolation. The "actual" water table elevation was ground elevation minus the minimum measured depth-to-water. The predicted water table elevation was the value extracted from the raster output of the radial basis function interpolation at the well point location. The plot includes points within the final study area only (i.e., within 1 km of the Bay edge). Figure 4 shows the results of our depth-to-water modeling for the coastal Bay Area. A geospatial data file can be downloaded at https://datadryad.org/stash/dataset/doi:10.6078/D1W01Q We found that a shallow groundwater condition exists in many developed areas in the North Bay including Fairfield, Novato, San Rafael, and Petaluma, although fewer well data points were available for the North Bay in general. Many cities in the East Bay also had shallow groundwater along the Bay-front, placing major infrastructure (such as Interstate highways 580 and 880) at risk. Exposure to potential groundwater flooding was perhaps most severe in the Silicon Valley area, where the minimum depth-to-water was already less than one meter in large areas of Mountain View, Redwood City, and San Mateo. Figure 5 shows subset maps of groundwater conditions in selected highly urbanized areas.  . Minimum depth-to-water for the coastal San Francisco Bay Area. Shallow groundwater within one kilometer of the coast is shown in color, with the shallowest areas in red. Our method produced some negative values that suggested groundwater was already emergent, usually where there were no well points in the dataset at the base of a slope or in a valley. These areas (in black) most likely had very shallow groundwater with seasonal surface discharges, but a process-based model would be needed to quantify the volume. likely had very shallow groundwater with seasonal surface discharges, but a process-based model would be needed to quantify the volume.

Projecting Future Conditions
To determine the relationship between 1 m of SLR and a rising water table, we used a simple linear approximation within 1 km of the Bay edge. This replicated the distance used by Rotzoll and Fletcher [5] in Honolulu, HI based on measured tidal efficiencies, and by Hoover et al. [4] in three smaller areas along the California coast. Hoover et al. [4] described this linear approximation of the effect of rising sea levels on groundwater depth as conservative, because additional tidal effects would only increase the impacts on groundwater emergence and shoaling at high tides. Using this linear approximation of sea level rise impacts, areas of the map in Figure 4 where minimum depth-to-groundwater was less than one meter would likely experience groundwater emergence during the wet season of wet years with one meter of SLR. The State of California recommended that public agencies consider one meter of SLR likely at the Golden Gate by 2100 under the RCP 8.5 IPCC emissions scenario, and by 2150 under the RCP 4.5 and RCP 2.6 scenarios [33]. Figure 6 shows the minimum depth-to-groundwater in the highly urbanized Bay Area with 1 m of SLR. Our analysis revealed widespread areas where surface flooding from groundwater emergence is possible. Table 3 reports the extent of flooding from emergent groundwater with 1 m of SLR, compared to the extent of direct flooding from the Bay with the same SLR, based on projections from the Our Coast, Our Future Flood Map (USGS CoSMoS flood model) [20]. To match the groundwater study area, we excluded areas more than 1 km away from a groundwater monitoring well from the direct SLR flooded area calculation.

Projecting Future Conditions
To determine the relationship between 1 m of SLR and a rising water table, we used a simple linear approximation within 1 km of the Bay edge. This replicated the distance used by Rotzoll and Fletcher [5] in Honolulu, HI based on measured tidal efficiencies, and by Hoover et al. [4] in three smaller areas along the California coast. Hoover et al. [4] described this linear approximation of the effect of rising sea levels on groundwater depth as conservative, because additional tidal effects would only increase the impacts on groundwater emergence and shoaling at high tides. Using this linear approximation of sea level rise impacts, areas of the map in Figure 4 where minimum depthto-groundwater was less than one meter would likely experience groundwater emergence during the wet season of wet years with one meter of SLR. The State of California recommended that public agencies consider one meter of SLR likely at the Golden Gate by 2100 under the RCP 8.5 IPCC emissions scenario, and by 2150 under the RCP 4.5 and RCP 2.6 scenarios [33]. Figure 6 shows the minimum depth-to-groundwater in the highly urbanized Bay Area with 1 m of SLR. Our analysis revealed widespread areas where surface flooding from groundwater emergence is possible. Table 3 reports the extent of flooding from emergent groundwater with 1 m of SLR, compared to the extent of direct flooding from the Bay with the same SLR, based on projections from the Our Coast, Our Future Flood Map (USGS CoSMoS flood model) [20]. To match the groundwater study area, we excluded areas more than 1 km away from a groundwater monitoring well from the direct SLR flooded area calculation.   The results of our analysis, based on an interpolation of empirical groundwater well data and a linear relationship between SLR and groundwater levels, can be used to identify hotspots that require a second phase of analysis using higher-resolution elevation and hydrologic data, field measurements of tidal efficiency, and process-based models. Process-based models developed at smaller geographic scales may be able to account for recharge and discharge, the diminishing influence of SLR inland from the coast, wave run-up, and variations in geologic and infrastructure conditions.

Discussion
We created an interpolated surface that estimated the depth of shallow groundwater for 489 square kilometers of San Francisco Bay's coastline using measured depth-to-water and tidal data. This rapid assessment method indicated that many parts of the Bay Area coastline are vulnerable to rising groundwater. Based on these results, many San Francisco Bay Area communities should conduct further modeling studies to prepare for potential flooding from groundwater, in addition to direct flooding from SLR. Our study suggested that there is significant potential for groundwater flooding in important Silicon Valley economic hubs (e.g., Mountain View, East Palo Alto, Redwood City), East Bay cities with fast-growing populations (e.g., Oakland, Hayward, Fremont), and major transportation infrastructure, including freeways (e.g., Interstate 580) and airports (Oakland International Airport, San Francisco International Airport). Our results indicate that flooding from emergent groundwater could impact more land by area than direct SLR flooding, with a SLR scenario of one meter in seven of the nine Bay Area counties, and in the region as a whole (Table 3). However, the calculated area impacted by emergent groundwater does not account for surface discharge to streams and other water bodies.
In addition to groundwater emergence, risks posed to developed areas include increased infiltration and inflow of underground water and wastewater pipes [15], and increased liquefaction risks in active seismic zones. Rising groundwater can also mobilize contaminants from wastewater and legacy soil pollution, producing human and ecosystem health risks. Groundwater emergence is likely to occur even where levees and seawalls are built to serve as barriers to saltwater coastal inundation. These structures alone will be inadequate to prevent flooding and other hazards without additional adaptation measures.

Conclusions
We used a rapid assessment interpolation method to create the first depth-to-groundwater map for the San Francisco Bay shore zone. The empirical data we used reflects existing human impacts from pumping, storm sewer infiltration, and leaky water pipes in a complex urban environment. We maximized the accuracy of the interpolated surface map by testing a variety of methods and selecting the one that minimized errors. The results of the analysis revealed widespread shallow groundwater conditions along most of the shore of San Francisco Bay. Using the conservative assumption of a linear relationship between SLR and shallow, unconfined groundwater depth within one kilometer of the shore, we showed that many densely developed areas are at risk from rising and even emergent groundwater as the sea level rises.
The method presented here is useful as a rapid assessment technique for comparing relative exposure to groundwater hazards and identifying hotspots where localized dynamic modeling is needed [15,34,35]. The minimum depth-to-water surface shown here did not represent a particular point in time, but rather an estimate based on the shallowest measurement taken at each monitoring well in the dataset during the study timeframe. Sampling was not consistent over time in this best-available dataset. Therefore, seasonal changes in precipitation and infiltration were not captured by this minimum depth-to-water method, although they are an important consideration [30]. Since it was empirical rather than modeled, the dataset we used for this interpolation reflected human impacts on coastal groundwater, including current pumping and leaky pipes. In many areas, the results shown here were influenced by local leachate pumping at landfills or other groundwater pumping that was already in place to prevent flooding.
Any interpolation-based method contains errors. More consistent sample point coverage would have reduced the level of error introduced by interpolation. Additionally, the simple linear approximation we used to estimate rising groundwater levels due to SLR did not account for a number of factors that would have been important to consider in a more nuanced modeling effort. Additional factors to consider in future refinements of this technique include the diminishing influence of SLR inland from the coast, the potential effects of tides, waves, and extreme rainfall events, and the need for more accurate local measurements to establish the effects of different geologic conditions and underground pipe and pump systems on the level of groundwater rise. Modeling efforts incorporating measurements of tidal influence and groundwater flow, such as those that Habel et al. conducted in regard to Honolulu, HI [36], are needed in the areas that were identified as potential hotspots by our method.
While previous studies established the existence of rising groundwater due to SLR [2,[4][5][6]9,10,22,31], as well as the potential impact at case study sites [2,4,5,28,31,37], this paper provides a method for building a regional-scale view of the potentially widespread impacts on surface flooding, underground infrastructure, and the health of people and ecosystems. Understanding the full range of SLR impacts is essential for prioritizing adaptation investments, and selecting appropriate strategies in coastal cities [15,35,36]. Other low-lying urban areas around the world with shallow and unconfined coastal aquifers have an urgent need to identify the potential for future groundwater flooding as a result of sea level rise. In eastern and southeastern US, major metropolitan regions around rivers and bays such as Boston, New York, Philadelphia, Baltimore, Washington DC, Norfolk, Charleston, Ft. Lauderdale, Miami, Tampa, and Galveston could benefit from similar assessment methods for groundwater flooding that make use of existing groundwater quality datasets. On the west coast, Seattle, Tacoma, and many smaller cities and towns on bays along the Oregon, Washington, and California coasts are likely to face groundwater flooding. Many low-lying cities along bays and deltas in northwestern Europe, coastal areas of the United Kingdom, coastal Africa, South America, and Southeast Asia face similar threats. The rapid assessment method presented here provides a valuable approach for the identification of hotspots where rising groundwater poses a threat to urban development and human health. Once hotspots are identified, process-based groundwater data collection and modeling efforts will be needed at a local scale to more fully represent the dynamics of rising groundwater in coastal zones and to account for variables such as projected future changes in subsidence, recharge, and discharge rates.