Metrics for marsh migration under sea-level rise in Chesapeake Bay

This process step and subsequent process steps were performed by the same person, Zafer Defne, in ArcGIS Pro (ver.3.3.2) unless otherwise stated. For complex operations, names of specific tools used are given in CAPITAL letters (any critical parameters used are given in parentheses, separated by a semicolon, immediately after the tool name). The input and output file names are provided in [square brackets] when necessary. Units for length and area calculations are meters (m) and square meters (m2) unless otherwise stated.
This processing step defines the extent of migration area.
a) Download the NOAA Sea level rise wetland impacts and migration datasets for the Chesapeake Bay domain [md_marshmigration_2016.zip, va_marshmigration_2016.zip, nc_marshmigration_2016.zip]. These datasets present the potential distribution of each wetland type based on their elevation and frequency of inundation under different SLR rates. Use four scenarios: 0.5, 1.0, 1.5, 2.0 feet of SLR as input.
b) From each scenario select classes of Palustrine Emergent Wetland (15), Brackish/Transition Wetland (17), Estuarine Wetland (18) by EXTRACT_BY_ATTRIBUTES(VALUE = 15 Or VALUE = 17 Or VALUE = 18).
c) Mosaic all four rasters from previous step into to a new raster for Maryland (MD) while setting the coordinate system to NAD1983 (2011) UTM Zone 18N. MOSAIC_TO_NEW_RASTER(Pixel Type="1bit"; Spatial Reference=NAD_1983_2011_UTM_Zone_18N; Mosaic Operator="Maximum"). Export only the values greater than 1 as a new raster and convert to polygon. RASTER_TO_POLYGON(Field = "VALUE").
Repeat the same steps for Virginia (VA) and North Carolina (NC) domains. And merge the results in to a single feature dataset [mm_MD_VA_NC_merge_UTM18_2011.shp].

This processing step delineates the migration area based on the USGS National Hydrography Dataset Plus High Resolution (NHD Plus HR) catchments map. First, the catchments intersecting the marsh units are selected as immediate catchments. Then, if the inland boundary of the migration layer spans across multiple catchments, those catchments are merged to the immediate catchments by proximity. Catchments areas beyond the inland boundary of migration layer are clipped to the boundary. The resulting product is the marsh migration polygons which is the migration zone delineated with catchments.
a) Download Chesapeake Bay marsh units dataset [mu_UVVR_CB.shp] and NHD Plus HR catchments datasets for the domain. [mu_UVVR_CB.shp]; NHDPlusCatchment feature datasets [NHDPlusCatchment_0206.shp], [NHDPlusCatchment_0207.shp], [NHDPlusCatchment_0208.shp] and [NHDPlusCatchment_0301.shp].
b) Project all source feature datasets to the same coordinate system with the previous processing step, NAD1983 (2011) UTM Zone 18N. Resulting projected data sets are: Marsh units [mu_UVVR_CB_UTM18_2011.shp]; NHDPlusCatchment feature datasets [NHDPlusCatchment_0206_UTM18_2011.shp], [NHDPlusCatchment_0207_UTM18_2011.shp], [NHDPlusCatchment_0208_UTM18_2011.shp] and [NHDPlusCatchment_0301_UTM18_2011.shp]
c) Merge all catchment feature datasets to create the total catchment area [HDPlusCatchment_merge_0206_0207_0208_0301.shp].
d) Select parts of migration area that intersect with marsh units. SELECT_LAYER_BY_LOCATION(Input features=mm_MD_VA_NC_merge_UTM18_2011.shp; Selecting features=mu_UVVR_CB_UTM18_2011.shp; Relationship=Intersect; Search distance =1 meter). Export features as [mmp_by_mu.shp].
e) SELECT_LAYER_BY_LOCATION(Input features=HDPlusCatchment_merge_0206_0207_0208_0301.shp; Selecting features= mmp_by_mu.shp, Relationship=Intersect). Export features as [cathments_selected.shp].
f) Intersect the two feature datasets to identify the immediate catchment polygons. PAIRWISE_INTERSECT(Input features=catchments_selected.shp, mmp_by_mu.shp ). Export features as [mmp_iniside_cat.shp].
g) Remove the immediate catchments from the overall marsh migration area so that the catchments in the migration area but beyond the immediate ones can be processed. PAIRWISE_ERASE(Input features=mm_MD_VA_NC_merge_UTM18_2011.shp; Erase features=mmp_inside_cat.shp; Output=mm_outside_cat.shp). MULTIPART_TO_SINGLEPART(Input features=mm_outside_cat.shp; Output features=mm_outside_cat_single.shp) and SELECT_LAYER_BY_LOCATION(Input features=mm_outside_cat_single.shp; Selecting features=mu_UVVR_CB_UTM18_2011.shp; Relationship=Intersect; Search distance =1 meter). Export features as [mmp_outside_cat.shp].
h) Union the immediate catchments with the remote catchments to create a union set of catchments. UNION(Input features=mmp_iniside_cat.shp, mmp_outside_cat.shp; Gaps allowed=Yes, Output features=mmp_all.shp).
i) Prepare and split the union of catchments to immediate catchments and remote catchments so that the remote catchments can be merged by proximity. Specifically, ERASE marsh units from the union of catchments and MULTIPART_TO_SINGLEPART to obtain [mmp_all_single.shp]. Create two feature datasets from this dataset: For immediate catchments SELECT_LAYER_BY_LOCATION(Input features=mmp_all_single.shp; Selecting features=mu_UVVR_CB_UTM18_2011.shp; Relationship=Intersect; Search distance=0.001 meter) and export as [input_mmu.shp]. For remote catchments SELECT_LAYER_BY_LOCATION(Input features=mmp_all_single.shp; Selecting features=mu_UVVR_CB_UTM18_2011.shp; Relationship=Intersect; Search distance=0.001 meter, Invert Spatial Relationship=Yes) and export as [near_mmu.shp].
j) Run the HYDUNITLOOP tool to merge catchments. At each iteration step, the script uses NEAR tool to find remote catchments from the [input_mmu.shp] within 1 meter of an immediate catchment [near_mmu.shp] and merges them to it using UNION and DISSOLVE tools. The catchments merged with an immediate catchment are removed from the remote catchments dataset. The script iterates until there is no change in the number of remote and immediate units. The remaining catchments are those that are more than 1 meter away from any immediate ones and they are discarded. The resulting final dataset is [final_mmu.shp]. HYDUNITLOOP(Input features=input_mmu.shp; Near features=near_mmu.shp; Output features=final_mmu.shp].
k) Because of the resolution difference between the datasets, processing sometimes creates sliver polygons around the edges. To clean up the final output, remove detached marsh migration polygons with area smaller than a threshold value. Specifically, first PAIRWISE_DISSOLVE(Input features=final_mmu.shp; Output feature class=final_mmu_diss.shp) and CALCULATE_GEOMETRY_ATTRIBUTES(Input features=final_mmu_diss.shp; Geometry attributes field=ASQM; Property=Area (geodesic); Area unit=Square meters; Coordinate system=NAD_1983_2011_Contiguous_USA_Albers). Then, select polygons with areas greater than 300 square meter to export [fmmu_gt_300.shp] and use them to select from [final_mmu.shp] dataset by SELECT_BY_LOCATION(Input feature=final_mmu.shp; Relationship=Within; Selecting features=fmmu_gt_300.shp; Search distance=0 meters). Finally, export the selected features as [marsh_migration_polygons.shp].

Calculate total marsh migration area, ratio of total migration area to marsh unit area and migration rate.
Total migration area for a marsh unit is calculated as the sum of all catchment areas within in the migration zone that the marsh unit borders, where catchments are clipped to the SLR extent landward. If same migration area is bounded by multiple marsh units then the marsh area proportioned by the length of the shared boundary between the migration area and each unit.
Marsh migration rate for each unit is defined as the migration area covered per year by that marsh unit under a SLR scenario. Similar to total area calculation, if same migration area is shared by multiple marsh units a migration rate is calculated per each migration area proportional to the shared boundary. However, if multiple catchments are within the total marsh migration area, the maximum of the migration rates is assigned as the final migration rate.
a) First, calculate the migration rate in in terms of area per time by dividing the migration area for each unit by the time it takes for a SLR of 2 feet. This step is done in MATLAB (version 2021b). For sea-level rise projections Sweet and others (2022) SLR data within the region were used. A total of 40 points including 20 stations and 20 grid points were within the bounding box with lower left corner of (-78.0000N, 35.7950W) and upper right corner of (-74.9600N, 40.0000W). SLR rate is calculated as SLR_RATE=(RSL2100+RSL_OFFSET)/100-RSL_VLM, where RSL2100 is the relative sea level by 2100 under a global mean sea level rise scenario, RSL_OFFSET is the offset to initiate the projection at year 2000, and RSL_VLM is the relative sea level contribution from vertical land motion. Interpolate the SLR_RATE from 40 points over a regular grid using GRIDDATA function with natural neighbor interpolation (grid size of 30 longitudinal by 42 latitudinal points). Assign to each marsh migration polygon the interpolated SLR_RATE value from the nearest grid point. For each polygon, calculate the time it takes for a SLR of 2 feet by dividing the migration area by the SLR_RATE. Do this for the three scenarios considered: GMSL rise of 0.3 meters, 0.5 meters and 1.0 meters by year 2100. Join these values to the marsh migration polygons table as YRS03, YRS05, YRS10, respectively.
Sweet, W.V., Hamlington, B.D., Kopp, R.E., Weaver, C.P., Barnard, P.L., Bekaert, D., Brooks, W., Craghan, M., Dusek, G., Frederikse, T., Garner, G., Genz, A.S., Krasting, J.P., Larour, E., Marcy, D., Marra, J.J., Obeysekera, J., Osler, M., Pendleton, M., Roman, D., Schmied, L., Veatch, W., White, K.D., and Zuzak, C., 2022, Global and Regional Sea Level Rise Scenarios for the United States: Updated Mean Projections and Extreme Water Level Probabilities Along U.S. Coastlines. NOAA Technical Report NOS 01. National Oceanic and Atmospheric Administration, National Ocean Service, Silver Spring, MD, 111 pp.
b) Next, add a new field to marsh migration polygons dataset to transfer their FID info to the shared boundary segments. CALCULATE_FIELD(Input table=marsh_migration_polgons.shp; Field name=FID_MM; Field Type=Long; Expression="FID_MM=!FID!"). Then intersect the marsh units and the marsh polygons to calculate the length of each shared boundary. INTERSECT(Input features=mu_UVVR_CB.shp, marsh_migration_polygons.shp; Output feature=MM_MU_boundary.shp; Output type=Line). Calculate the length of each boundary segment with CALCULATE_GEOMETRY_ATTRIBUTES(Input features=MM_MU_boundary.shp; Geometry attributes field=M_L; Property=Length(geodesic); Length units=Meters; Coordinate system=NAD_1983_2011_UTM_Zone_18N).
c) Delete field FID_MM from marsh migration polygons and spatial join with the shared boundary dataset to create marsh polygons with boundary segment lengths. SPATIAL_JOIN(Target features=marsh_migration_polygons.shp; Join features=MM_MU_boundary.shp; Output feature class=MM_MU_boundary_len.shp; Join operation=One to many; Match option=Intersect; Fields=FID_MM,FID_CMU,M_L,YRS03,YRS05,YRS10). SELECT(Input features=MM_MU_boundary_len.shp; Output feature class=MM_MU_boundary_len_Select.shp; Where=TARGET_FID is equal to FID_MM).
d) Calculate marsh migration area associated with each boundary segment by CALCULATE_GEOMETRY_ATTRIBUTES(Input features=MM_MU_boundary_len_Select.shp; Geometry attributes field=MM_A; Property=Area (geodesic); Area unit=Square meters; Coordinate system=NAD_1983_2011_Contiguous_USA_Albers) and SUMMARY_STATISTICS(Input table=MM_MU_boundary_len_Select.shp; Output table=M_L_sum_by_MM.dbf; Statistics field property=M_L; Statistics type=Sum; Case field=TARGET_FID) and JOIN_FIELD(Input table=MM_MU_boundary_len_Select.shp; Input field=TARGET_FID; Join table=M_L_sum_by_MM.dbf; Join field=TARGET_FID; Transfer method=Use field mapping; Field map=TARGET_FID, SUM_M_L; Index join fields=Do not add indexes).
e) Calculate potential marsh migration area associated with each marsh unit by CALCULATE_FIELD(Input table=MM_MU_boundary_len_Select.shp; Field name=PMA; Expression="!M_L! / !SUM_M_L! * !MM_A!") and then PAIRWISE_DISSOLVE(Input features=MM_MU_boundary_len_Select.shp; Output feature class=MM_MU_boundary_area_summary.shp; Dissolve fields=FID_CMU; Statistics fields=PMA; Statistics type=Sum; Create multipart features=Yes).
f) Finally, create the marsh migration dataset by joining the migration summary dataset to the UVVR dataset. Also calculate a ratio of migration area to marsh unit are at this step. JOIN_FIELD(Input table=mu_UVVR_CB.shp; Input field=FID_CMU; Join table=MM_MU_boundary_area_summary.shp; Join field=FID_CMU; Transfer method=Use field mapping; FieldMap:AMIG_M2,MIGR03,MIGR05,MIGR10). CALCULATE_FIELD(Input table=mu_UVVR_CB.shp; Field name=AMIGRAT; Expression="!AMIG_M2!/!ATOT_M2!"). Keep only FID_CMU, ATOT_M2, AVEG_M2, AMIG_M2, FLG fields and the transferred fields, and export features as mu_migration_CB.shp.