Rising sea levels (SLR) will cause coastal groundwater to rise in many coastal urban environments. Inundation of contaminated soils by groundwater rise (GWR) will alter the physical, biological, and geochemical conditions that influence the fate and transport of existing contaminants. These transformed products can be more toxic and/or more mobile under future conditions driven by SLR and GWR. We reviewed the vulnerability of contaminated sites to GWR in a US national database and in a case comparison with the San Francisco Bay region to estimate the risk of rising groundwater to human and ecosystem health. The results show that 326 sites in the US Superfund program may be vulnerable to changes in groundwater depth or flow direction as a result of SLR, representing 18.1 million hectares of contaminated land. In the San Francisco Bay Area, we found that GWR is predicted to impact twice as much coastal land area as inundation from SLR alone, and 5,297 state-managed sites of contamination may be vulnerable to inundation from GWR in a 1-meter SLR scenario. Increases of only a few centimeters of elevation can mobilize soil contaminants, alter flow directions in a heterogeneous urban environment with underground pipes and utility trenches, and result in new exposure pathways. Pumping for flood protection will elevate the salt water interface, changing groundwater salinity and mobilizing metals in soil. Socially vulnerable communities are more exposed to this risk at both the national scale and in a regional comparison with the San Francisco Bay Area. Methods Data Dryad This data set includes data from the California State Water Resources Control Board (WRCB), the California Department of Toxic Substances Control (DTSC), the USGS, the US EPA, and the US Census. National Assessment Data Processing: For this portion of the project, ArcGIS Pro and RStudio software applications were used. Data processing for superfund site contaminants in the text and supplementary materials was done in RStudio using R programming language. RStudio and R were also used to clean population data from the American Community Survey. Packages used include: Dplyr, data.table, and tidyverse to clean and organize data from the EPA and ACS. ArcGIS Pro was used to compute spatial data regarding sites in the risk zone and vulnerable populations. DEM data processed for each state removed any elevation data above 10m, keeping anything 10m and below. The Intersection tool was used to identify superfund sites within the 10m sea level rise risk zone. The Calculate Geometry tool was used to calculate the area within each coastal state that was occupied by the 10m SLR zone and used again to calculate the area of each superfund site. Summary Statistics were used to generate the total proportion of superfund site surface area / 10m SLR area for each state. To generate population estimates of socially vulnerable households in proximity to superfund sites, we followed methods similar to that of Carter and Kalman (2020). First, we generated buffers at the 1km, 3km, and 5km distance of superfund sites. Then, using Tabulate Intersection, the estimated population of each census block group within each buffer zone was calculated. Summary Statistics were used to generate total numbers for each state. Bay Area Data Processing: In this regional study, we compared the groundwater elevation projections by Befus et al (2020) to a combined dataset of contaminated sites that we built from two separate databases (Envirostor and GeoTracker) that are maintained by two independent agencies of the State of California (DTSC and WRCB). We used ArcGIS to manage both the groundwater surfaces, as raster files, from Befus et al (2020) and the State’s point datasets of street addresses for contaminated sites. We used SF BCDC (2020) as the source of social vulnerability rankings for census blocks, using block shapefiles from the US Census (ACS) dataset. In addition, we generated isolines that represent the magnitude of change in groundwater elevation in specific sea level rise scenarios. We compared these isolines of change in elevation to the USGS geological map of the San Francisco Bay region and noted that groundwater is predicted to rise farther inland where Holocene paleochannels meet artificial fill near the shoreline. We also used maps of historic baylands (altered by dikes and fill) from the San Francisco Estuary Institute (SFEI) to identify the number of contaminated sites over rising groundwater that are located on former mudflats and tidal marshes. The contaminated sites' data from the California State Water Resources Control Board (WRCB) and the Department of Toxic Substances (DTSC) was clipped to our study area of nine-bay area counties. The study area does not include the ocean shorelines or the north bay delta area because the water system dynamics differ in deltas. The data was cleaned of any duplicates within each dataset using the Find Identical and Delete Identical tools. Then duplicates between the two datasets were removed by running the intersect tool for the DTSC and WRCB point data. We chose this method over searching for duplicates by name because some sites change names when management is transferred from DTSC to WRCB. Lastly, the datasets were sorted into open and closed sites based on the DTSC and WRCB classifications which are shown in a table in the paper's supplemental material. To calculate areas of rising groundwater, we used data from the USGS paper “Projected groundwater head for coastal California using present-day and future sea-level rise scenarios” by Befus, K. M., Barnard, P., Hoover, D. J., & Erikson, L. (2020). We used the hydraulic conductivity of 1 condition (Kh1) to calculate areas of rising groundwater. We used the Raster Calculator to subtract the existing groundwater head from the groundwater head under a 1-meter of sea level rise scenario to find the areas where groundwater is rising. Using the Reclass Raster tool, we reclassified the data to give every cell with a value of 0.1016 meters (4”) or greater a value of 1. We chose 0.1016 because groundwater rise of that little can leach into pipes and infrastructure. We then used the Raster to Poly tool to generate polygons of areas of groundwater rise.

pockmarks are defined as depressions on the seabed and are usually formed by fluid expulsions. recently discovered, pockmarks along the aquitaine slope within the french eez, were manually mapped although two semi-automated methods were tested without convincing results. in order to potentially highlight different groups and possibly discriminate the nature of the fluids involved in their formation and evolution, a morphological study was conducted, mainly based on multibeam data and in particular bathymetry from the marine expedition gazcogne1, 2013. bathymetry and seafloor backscatter data, covering more than 3200 km², were acquired with the kongsberg em302 ship-borne multibeam echosounder of the r/v le suroît at a speed of ~8 knots, operated at a frequency of 30 khz and calibrated with ©sippican shots. precision of seafloor backscatter amplitude is +/- 1 db. multibeam data, processed using caraibes (©ifremer), were gridded at 15x15 m and down to 10x10 m cells, for bathymetry and seafloor backscatter, respectively. the present table includes 11 morphological attributes extracted from a geographical information system project (mercator 44°n conserved latitude in wgs84 datum) and additional parameters related to seafloor backscatter amplitudes. pockmark occurrence with regards to the different morphological domains is derived from a morphological analysis manually performed and based on gazcogne1 and bobgeo2 bathymetric datasets.the pockmark area and its perimeter were calculated with the “calculate geometry” tool of arcmap 10.2 (©esri) (https://desktop.arcgis.com/en/arcmap/10.3/manage-data/tables/calculating-area-length-and-other-geometric-properties.htm). a first method to calculate pockmark internal depth developed by gafeira et al. was tested (gafeira j, long d, diaz-doce d (2012) semi-automated characterisation of seabed pockmarks in the central north sea. near surface geophysics 10 (4):303-315, doi:10.3997/1873-0604.2012018). this method is based on the “fill” function from the hydrology toolset in spatial analyst toolbox arcmap 10.2 (©esri), (https://pro.arcgis.com/en/pro-app/tool-reference/spatial-analyst/fill.htm) which fills the closed depressions. the difference between filled bathymetry and initial bathymetry produces a raster grid only highlighting filled depressions. thus, only the maximum filling values which correspond to the internal depths at the apex of the pockmark were extracted. for the second method, the internal pockmark depth was calculated with the difference between minimum and maximum bathymetry within the pockmark.latitude and longitude of the pockmark centroid, minor and major axis lengths and major axis direction of the pockmarks were calculated inside each depression with the “zonal geometry as table” tool from spatial analyst toolbox in arcgis 10.2 (©esri) (https://pro.arcgis.com/en/pro-app/tool-reference/spatial-analyst/zonal-statistics.htm). pockmark elongation was calculated as the ratio between the major and minor axis length.cell count is the number of cells used inside each pockmark to calculate statistics (https://pro.arcgis.com/en/pro-app/tool-reference/spatial-analyst/zonal-geometry.htm). cell count and minimum, maximum and mean bathymetry, slope and seafloor backscatter values were calculated within each pockmark with “zonal statistics as table” tool from spatial analyst toolbox in arcgis 10.2 (©esri). slope was calculated from bathymetry with “slope” function from spatial analyst toolbox in arcgis 10.2 (©esri) and preserves its 15 m grid size (https://pro.arcgis.com/en/pro-app/tool-reference/spatial-analyst/slope.htm). seafloor backscatter amplitudes (minimum, maximum and mean values) of the surrounding sediments were calculated within a 100 m buffer around the pockmark rim.

Reason for SelectionProtected natural areas in urban environments provide urban residents a nearby place to connect with nature and offer refugia for some species. Because beaches in Puerto Rico and the U.S. Virgin Islands are open to the public, beaches also provide important outdoor recreation opportunities for urban residents, so we include beaches as parks in this indicator.Input DataSoutheast Blueprint 2023 subregions: CaribbeanSoutheast Blueprint 2023 extentNational Oceanic and Atmospheric Administration (NOAA) National Centers for Environmental Information (NCEI) Coastal Relief Model, accessed 11-22-2022Protected Areas Database of the United States (PAD-US) 3.0: VI, PR, and Marine Combined Fee EasementPuerto Rico Protected Natural Areas 2018 (December 2018 update): Terrestrial and marine protected areas (PACAT2018_areas_protegidasPR_TERRESTRES_07052019.shp, PACAT2018_areas_protegidasPR_MARINAS_07052019.shp) 2020 Census Urban Areas from the Census Bureau’s urban-rural classification; download the data, read more about how urban areas were redefined following the 2020 censusOpenStreetMap data “multipolygons” layer, accessed 3-14-2023A polygon from this dataset is considered a park if the “leisure” tag attribute is either “park” or “nature_reserve”, and considered a beach if the value in the “natural” tag attribute is “beach”. OpenStreetMap describes leisure areas as “places people go in their spare time” and natural areas as “a wide variety of physical geography, geological and landcover features”. Data were downloaded in .pbf format and translated ton an ESRI shapefile using R code. OpenStreetMap® is open data, licensed under the Open Data Commons Open Database License (ODbL) by the OpenStreetMap Foundation (OSMF). Additional credit to OSM contributors. Read more on the OSM copyright page. TNC Lands - Public Layer, accessed 3-8-2023U.S. Virgin Islands beaches layer (separate vector layers for St. Croix, St. Thomas, and St. John) provided by Joe Dwyer with Lynker/the NOAA Caribbean Climate Adaptation Program on 3-3-2023 (contact jdwyer@lynker.com for more information)Mapping StepsMost mapping steps were completed using QGIS (v 3.22) Graphical Modeler.Fix geometry errors in the PAD-US PR data using Fix Geometry. This must be done before any analysis is possible.Merge the terrestrial PR and VI PAD-US layers.Use the NOAA coastal relief model to restrict marine parks (marine polygons from PAD-US and Puerto Rico Protected Natural Areas) to areas shallower than 10 m in depth. The deep offshore areas of marine parks do not meet the intent of this indicator to capture nearby opportunities for urban residents to connect with nature.Merge into one layer the resulting shallow marine parks from marine PAD-US and the Puerto Rico Protected Natural Areas along with the combined terrestrial PAD-US parks, OpenStreetMap, TNC Lands, and USVI beaches. Omit from the Puerto Rico Protected Areas layer the “Zona de Conservación del Carso”, which has some policy protections and conservation incentives but is not formally protected.Fix geometry errors in the resulting merged layer using Fix Geometry.Intersect the resulting fixed file with the Caribbean Blueprint subregion.Process all multipart polygons to single parts (referred to in Arc software as an “explode”). This helps the indicator capture, as much as possible, the discrete units of a protected area that serve urban residents.Clip the Census urban area to the Caribbean Blueprint subregion.Select all polygons that intersect the Census urban extent within 1.2 miles (1,931 m). The 1.2 mi threshold is consistent with the average walking trip on a summer day (U.S. DOT 2002) used to define the walking distance threshold used in the greenways and trails indicator. Note: this is further than the 0.5 mi distance used in the continental version of the indicator. We extended it to capture East Bay and Point Udall based on feedback from the local conservation community about the importance of the park for outdoor recreation.Dissolve all the park polygons that were selected in the previous step.Process all multipart polygons to single parts (“explode”) again.Add a unique ID to the selected parks. This value will be used to join the parks to their buffers.Create a 1.2 mi (1,931 m) buffer ring around each park using the multiring buffer plugin in QGIS. Ensure that “dissolve buffers” is disabled so that a single 1.2 mi buffer is created for each park.Assess the amount of overlap between the buffered park and the Census urban area using overlap analysis. This step is necessary to identify parks that do not intersect the urban area, but which lie within an urban matrix. This step creates a table that is joined back to the park polygons using the UniqueID.Remove parks that had ≤2% overlap with the urban areas when buffered. This excludes mostly non-urban parks that do not meet the intent of this indicator to capture parks that provide nearby access for urban residents. Note: In the continental version of this indicator, we used a threshold of 10%. In the Caribbean version, we lowered this to 2% in order to capture small parks that dropped out of the indicator when we extended the buffer distance to 1.2 miles.Calculate the GIS acres of each remaining park unit using the Add Geometry Attributes function.Join the buffer attribute table to the previously selected parks, retaining only the parks that exceeded the 2% urban area overlap threshold while buffered. Buffer the selected parks by 15 m. Buffering prevents very small parks and narrow beaches from being left out of the indicator when the polygons are converted to raster.Reclassify the polygons into 7 classes, seen in the final indicator values below. These thresholds were informed by park classification guidelines from the National Recreation and Park Association, which classify neighborhood parks as 5-10 acres, community parks as 30-50 acres, and large urban parks as optimally 75+ acres (Mertes and Hall 1995).Export the final vector file to a shapefile and import to ArcGIS Pro.Convert the resulting polygons to raster using the ArcPy Polygon to Raster function. Assign values to the pixels in the resulting raster based on the polygon class sizes of the contiguous park areas.Clip to the Caribbean Blueprint 2023 subregion.As a final step, clip to the spatial extent of Southeast Blueprint 2023. Note: For more details on the mapping steps, code used to create this layer is available in the Southeast Blueprint Data Download under > 6_Code. Final indicator valuesIndicator values are assigned as follows:6 = 75+ acre urban park5 = >50 to <75 acre urban park4 = 30 to <50 acre urban park3 = 10 to <30 acre urban park2 = 5 to <10 acre urban park1 = <5 acre urban park0 = Not identified as an urban parkKnown IssuesThis indicator does not include park amenities that influence how well the park serves people and should not be the only tool used for parks and recreation planning. Park standards should be determined at a local level to account for various community issues, values, needs, and available resources. This indicator includes some protected areas that are not open to the public and not typically thought of as “parks”, like mitigation lands, private easements, and private golf courses. While we experimented with excluding them using the public access attribute in PAD, due to numerous inaccuracies, this inadvertently removed protected lands that are known to be publicly accessible. As a result, we erred on the side of including the non-publicly accessible lands.This indicator includes parks and beaches from OpenStreetMap, which is a crowdsourced dataset. While members of the OpenStreetMap community often verify map features to check for accuracy and completeness, there is the potential for spatial errors (e.g., misrepresenting the boundary of a park) or incorrect tags (e.g., labelling an area as a park that is not actually a park). However, using a crowdsourced dataset gives on-the-ground experts, Blueprint users, and community members the power to fix errors and add new parks to improve the accuracy and coverage of this indicator in the future.Other Things to Keep in MindThis indicator calculates the area of each park using the park polygons from the source data. However, simply converting those park polygons to raster results in some small parks and narrow beaches being left out of the indicator. To capture those areas, we buffered parks and beaches by 15 m and applied the original area calculation to the larger buffered polygon, so as not to inflate the area by including the buffer. As a result, when the buffered polygons are rasterized, the final indicator has some areas of adjacent pixels that receive different scores. While these pixels may appear to be part of one contiguous park or suite of parks, they are scored differently because the park polygons themselves are not actually contiguous. The Caribbean version of this indicator uses a slightly different methodology than the continental Southeast version. It includes parks within a 1.2 mi distance from the Census urban area, compared to 0.5 mi in the continental Southeast. We extended it to capture East Bay and Point Udall based on feedback from the local conservation community about the importance of the park for outdoor recreation. Similarly, this indicator uses a 2% threshold of overlap between buffered parks and the Census urban areas, compared to a 10% threshold in the continental Southeast. This helped capture small parks that dropped out of the indicator when we extended the buffer distance to 1.2 miles. Finally, the Caribbean version does not use the impervious surface cutoff applied in the continental Southeast because the landcover data available in the Caribbean does not assess percent impervious in a comparable way.Disclaimer: Comparing with Older Indicator VersionsThere are numerous problems with using Southeast Blueprint