EPA/600/R-15/003 I September 2015 I www2.epa.gov/research
United States
Environmental Protection
Agency
Predictive Seagrass
Habitat Model
Office of Research and Development
National Health and Environmental Effects Research Laboratory
-------
-------
EPA/600/R-15/0031 September 2015
Predictive Seagrass Habitat Model
Prepared by:
Naomi E. Detenbeck and Steven Rego
Atlantic Ecology Division
National Health and Environmental Effects Research Laboratory
Narragansett, RI 02882
National Health and Environmental Effects Research Laboratory
Office of Research and Development
U.S. Environmental Protection Agency
Washington, DC 20460
-------
Notice
The development of the model described in this document has been funded by the U.S.
Environmental Protection Agency (EPA], with CIS support for input data preprocessing
received under SES3 BPA GS-35F-4594G, Task Order 1520 Remote Sensing and
Geographical Information Systems Support (RSGIS) TDD 2-1 Narragansett. This report has
been subjected to the Agency's peer and administrative review and has been approved for
publication. Mention of trade names or commercial products does not constitute
endorsement or recommendation for use.
Acknowledgements
James Hagy, Giancarlo Cicchetti, and Michael McManus provided comments on a preliminary
version of this report
-------
Table of Contents
Figures vi
Tables viii
Acronyms ix
Executive Summary for Coastal Managers xi
Abstract xii
1. Introduction 1
1.1 Purpose 1
1.2 Conceptual Model 2
1.3 Existing Approaches 2
1.3.1 Preliminary Transplant Suitability Index: 2
1.3.2 Predictive modeling approaches: 3
1.4 Study Site (Narragansett Bay) 4
2. Methods 6
2.1 Data Sources 6
2.2 Data Pre-processing 6
2.2.1 Salinity 6
2.2.2 Temperature 6
2.2.3 Sediment Type (Grain Size and Percent Total Organic Carbon) 6
2.2.4 Seagrass (Current/Historical) 8
2.2.5 Transparency 8
2.2.6 Wave Exposure Data - WEMo and WAVES 10
2.2.7 Distance to physical disturbance source: hardened shorelines and marinas 10
2.2.8 Unsewered residential development on high infiltration soils 11
2.2.9 Canada goose grazer density 11
2.3 Statistical Model Development 11
2.3.1 Strategies to limit memory requirements 11
2.3.2 Modeling approaches at different scales 12
2.3.3 Selection of best-fitting models 13
2.3.4 Model diagnostic tests 17
2.3.5 Development of coordinate framework to assess spatial autocorrelation 17
2.3.6 Assessment of spatial autocorrelation ranges 17
2.3.7 Eliminating or incorporating SAC into models 18
2.3.8 Strategies for eliminating multicollinearity 19
2.3.9 Alternative approaches 20
3. Results 21
3.1 Final Models 21
3.1.1 Seagrass Grid Presence/Absence 21
3.1.2 Predictive models for shoreline segment P/A and minimum or maximum depth
for seagrass 26
4. Discussion 46
4.1 Data and modeling challenges and constraints for statistical predictive seagrass models 46
4.2 Comparative performance of alternative modeling approaches 47
4.3 Relative importance of different factors in predicting seagrass P/A 49
4.3.1 Seagrass sensitivity to different environmental variables in Narragansett Bay 49
4.3.2 Seagrass sensitivity to different environmental variables in other regions 50
4.3.3 Limitations to existing regression models to predict seagrass in U.S. estuaries 51
4.3.4 Comparison of estimates for light compensation depth and optimum light levels 52
4.3.5 Seagrass sensitivity to factors other than transparency 53
-------
4.4 Model application for assessing management scenarios 58
4.5 Model predictions of seagrass increase with decreased nitrogen loading 61
4.6 Future improvements 62
5. Conclusions 65
References 68
R-package references 77
Appendices
A. National and Regional Data Sources for Seagrass Habitat Model Development A-l
B. Tutorial on Finding and Downloading Data for Seagrass Habitat Prediction Models
Using EPA's Estuary Data Mapper B-l
B.I. Overview of Estuary Data Mapper B-l
B.2. Installation of Estuary Data Mapper Tool B-2
B.3. Zooming to Area of Interest B-3
B.4. Identifying Data Sources B-7
B.4.1. System Boundaries B-7
B.4.2. Seagrass B-9
B.4.3. Depth B-10
B.4.4. Transparency B-ll
B.4.5. Energy Environment B-13
B.4.5.1. Wave Energy Model Inputs B-13
B.4.5.2. Relative Exposure B-13
B.4.5.3. Current velocity B-14
B.4.6. Sediment Characteristics B-14
B.4.7. Water Quality B-16
B.4.7.1. Temperature B-16
B.4.7.2. Salinity B-16
B.4.8. Nitrogen Concentration and Loading B-17
B.4.8.1. Nutrient Concentrations B-17
B.4.8.2. Atmospheric Loads B-17
B.4.8.3. Watershed Sources B-20
B.5. Downloading Data B-23
B.6. Terminating an EDM Session and EDM Updates B-24
C. R Packages and Commands Used in development of Predictive Seagrass Habitat Models C-l
C.I. Exploratory analysis (graphics package) C-3
C.2. General linear mixed models with diagnostics plots (MASS, rms packages) C-3
C.3. General linear mixed model with interaction plots (effects package) C-9
C.4. Generalized additive mixed models (mcgv package) C-9
C.5. Spline correlogram evaluations (ncf package) C-10
C.6. Community correlograms with anisotropy (CommunityCorrelogram package) C-ll
C.7. Cross-validation and ROC construction (pROC package) C-12
D. Exploratory Analyses and Diagnostic Tests D-l
D.I. Models for Seagrass Grid Cell Presence/Absence D-2
D.1.1. Exploratory Analyses D-2
D.I.2. Evaluation of Model Assumptions in Preliminary Seagrass Grid P/A Models D-8
D.2. Models for Shoreline Segment Presence/Absence D-13
D.2.1. Exploratory Analyses of Shoreline P/A Relationships D-13
D.2.2. Test of Model Assumptions during Initial Model Development for Shoreline P/A D-16
D.2.3. Test of Model Assumptions during Initial Model Development for Shoreline
Relative Abundance D-l7
IV
-------
D.2.4. Test of Model Assumptions during Initial Model Development for
Maximum Seagrass Depth D-20
E. Quality of the Data and Limitations on Use of the Data E-l
Appendices Figures and Tables
Figure B-l. Estuary Data Mapper home page B-2
Figure B-2. First tab of Estuary Data Mapper with default display extent B-3
Figure B-3. Selection for zooming in to a state on first tab of Estuary Data Mapper B-4
Figure B-4. Process of zooming into an individual estuary in Estuary Data Mapper B-5
Figure B-5. Estuary Data Mapper zoomed in to Narragansett Bay B-5
Figure B-6. Zooming in to an estuarine watershed in Estuary Data Mapper B-6
Figure B-7.Get Data tab of Estuary Data Mapper B-7
Figure B-8. Save Data tab on Estuary Data Mapper B-8
Figure B-9. Retrieving seagrass data in Estuary Data Mapper B-9
Figure B-10. Retrieving topobathymetry data in Estuary Data Mapper B-10
Figure B-ll. Selecting STORET variables for data retrieval in Estuary Data Mapper B-ll
Figure B-12. Retrieving remotely sensed light attentuation data for Chesapeake Bay
in Estuary Data Mapper B-12
Figure B-13. Retrieving Coastal Vulnerability Index data in Estuary Data Mapper B-13
Figure B-14. Retrieving current velocity data in Estuary Data Mapper B-14
Figure B-15. Drop-down selections for selection of sediment parameters from EPA National
Coastal Assessment dataset B-15
Figure B-16. Retrieval of remotely sensed temperature data with Estuary Data Mapper B-16
Figure B-17. Selection of remotely sensed salinity data for download in Estuary Data Mapper B-17
Figure B-18. Selection of nitrogen loading data sets in Estuary Data Mapper B-18
Figure B-19. Retrieval of CMAQ monthly nitrogen deposition data with Estuary Data Mapper B-19
Figure B-20. Retrieval of NADP nitrogen deposition data with Estuary Data Mapper B-19
Figure B-21. Nitrogen submenus for CMAQ and NADP nitrogen deposition B-20
Figure B-22. Submenus for nitrogen source data at gridded or estuarine scales B-21
Figure B-23. Selection of N loading data from SPARROW models B-22
Figure B-24. Save Data tab in Estuary Data Mapper, illustrating saving shapefiles B-23
Figure B-25. Final tab in Estuary Data Mapper to close out program and check for updates B-24
Figure C-l. General sequence of CIS and R analyses to create generalized linear mixed models
to predict seagrass presence/absence C-2
Figure D-l. Distribution of independent continuous variables entered into general linear models
and general additive models D-3
Figure D-2. Spine plots showing conditional probabilities of seagrass presence as a function
of potential independent predictor variables D-6
Figure D-3. Heterogeneity of variance for residuals of initial GLMM model to predict seagrass
presence at the scale of 10 meter grid cells D-10
Figure D-4. Spline correlogram with 95% confidence intervals from bootstrapping based on
Pearson residuals from working version of general linear mixed model D-ll
Figure D-5. Correlogram showing spatial autocorrelation for offshore distance based on
residuals from model predicting seagrass presence/absence for Shoreline 14 D-12
Figure D-6. Spine plots showing conditional probabilities of seagrass shoreline presence
(1) versus absence (0) as a function of potential independent predictor variables D-14
Figure D-7. Diagnostic plots for initial model 3 predicting shoreline seagrass presence/absence....D-18
Figure D-8. Diagnostics for model 3j D-19
Table A-l. National and regional data sources for seagrass habitat model development A-2
Table D-l. Potential independent variables included in models predicting seagrass presence D-9
-------
Figures
Figure 1. Map of Narragansett Bay and Rhode Island with color-coded bathymetry.
WWTF = Wastewater treatment facility 5
Figure 2 a. Map of zone delineating boundaries of seagrass presence/absence grid
(0 to 5 meter depth) in Narragansett Bay and b) close-up of 10-meter grid cells with seagrass
presence (red) and absence (black). In Figure 2b, the shoreline is to the south of the grids 7
Figure 3. Shoreline code assignment for Narragansett Bay. A value of -99 was assigned to
unused segments and isolated shorelines 9
Figure 4. Seagrass endpoints defined at different scales: presence/absence (P/A) at individual
grid cell centroids, P/A within transect perpendicular to shoreline (assessed after points were
snapped to shoreline), relative frequency of points along perpendicular transect, minimum depth
(z) of occurrence along transect, and maximum depth of occurrence along transect 12
Figure 5. Process for assigning x-, and y-coordinates to seagrass grid cell centroids 18
Figure 6. Diagnostic plots for model 1 predicting seagrass grid occupancy. Spatial autocorrelation of
residuals is virtually eliminated. Heterogeneity of variance is greatly reduced, a) Spline correlogram
of Pearson residuals for Shoreline 14. b) Plot of Pearson residuals versus predicted value 24
Figure 7. ROC curve for model 1 based on a) initial fit (area under curve = 0.9767) and
b) 10-fold cross-validation (area under the curve is 0.7144) 24
Figure 8. Nonlinear effects of a) centered Salinity (PSU), b) centered average Secchi depth 0 Water
depth (m), and c) centered sediment percent total organic carbon on odds ratio for seagrass
occurrence in grid cell. Values are calculated assuming average values for all co-variables not
represented in plot and across all shorelines with the reference sediment type for model 1 25
Figure 9. Diagnostic plots for model 2 predicting seagrass presence/absence by shoreline segment.
a) Pearson residual versus predicted value, b) Pearson residual by sediment class, and
c) Pearson residual versus centered percent total organic carbon 28
Figure 10. Receiver operating characteristic (ROC) curve for final model predicting shoreline
presence/absence of seagrass based on a) original model 2 fit with full data set (ROC = 0.9546),
b) same model with ten-fold model cross-validation (ROC = 0.9547), and c) original model 3 fit for
model with full data set incorporating historic 1999 eelgrass presence/absence (ROC = 0.9753) 31
Figure 11. Diagnostic plots for model 4 predicting average shoreline occupancy, a) Pearson residuals
vs predicted values, with loess curve superimposed showing no trend, b) Q-Q plot (standardized
deviance residuals versus theoretical quantiles), and c) scale-location plot (square root of
standardized deviance residuals versus predicted values), and d) residuals vs leverage with
three labeled outliers 33
Figure 12. Interaction plots showing effect of shoreline wind-derived wave energy by sediment
class for models 4 and 5 predicting average shoreline occupancy a) without (4) and
b) with historic seagrass presence (5) as predictors. 124561213 = Gravel, Sandy gravel,
Gravel-sand-silt, Sand, Gravelly sand, Sand-silt-clay, and Gravel-silt-clay, 7 = Silty sand, and
81011 = Silty, Sandy silt, and Clay silt 35
Figure 13. ROC curve showing fit of model 4 prediction of average shoreline occupancy
after removal of three outliers 36
Figure 14. Diagnostic plots for model 5 predicting average shoreline occupancy including
historic 1999 eelgrass predictors, a) Residuals vs predicted with loess plot overlaid.
b) Normal Q-Q plot, and c) Scale-Location plot. Note the presence of three outliers 38
Figure 15. ROC curve for model 5 a) before andb) after removal of three outliers 39
VI
-------
Figure 16. Interaction of sediment class and maximum wave mixing depth on minimum depth of
seagrass occurrence (model 6). SED4 class definitions are: 12 = sand-slit-clay; 810 = silty and
sandy-silt; 5 = sand; 7 = silty sand 40
Figure 17. Diagnostic plots for model predicting minimum seagrass depth of occurrence
(model 6) a) Residuals versus predicted values, b) Q-Q plot of residuals(dashed line shows
expectation for normal distribution), and c) residuals versus leverage, highlighting presence
of three outliers (labelled) 41
Figure 18. Smoothing functions in generalized additive model 7 to predict maximum depth of
seagrass occurrence by shoreline position, a) Centered loglO shoreline maximum average Secchi
depth (m), b) Centered loglO shoreline maximum unsewered residences on high infiltration
soils/km2, and c) Centered loglO shoreline maximum sediment total organic carbon (%) 43
Figure 19. Interactive effects of centered sediment total organic carbon and shoreline maximum
of seasonal minimum Secchi depth on maximum seagrass depth at shoreline index I (model 8).
Sediment TOC has a stimulatory effect at high Secchi depths but an inhibitory effect at low Secchi
depths. Secchi depth for each plot is indicated by brown vertical line in each tan horizontal header. ..45
Figure 20 Figure 20. Historic occurrences of seagrass in Narragansett Bay compiled by Kopp (1995) 56
Figure 21. Predicted 2010 maximum tidal currents (m/sec) during ebb tides (from NOAA tidal currents
web site: http://tidesandcurrents.noaa.gov/curr_pred.html). Points represent predicted values at
tidal stations. Shaded areas were interpolated within the 0-5 meter depth zone by inverse distance
weighting (2 point interpolation with shorelines as barriers to interpolation) 57
Figure 22. Probability of fine sediment resuspension based on USGS WAVES model for growing
season. Average probability of resuspension for limited nearshore areas where the average
probability greater than zero is 4.9% 58
Figure 23. Projected change in cumulative distribution function for probability of occurrence (x)
following 40% reduction in total N without (simulation 1) and with (simulation 2) sediment
recovery, a) All shorelines combined, b) Projected short-term effect of 40% TN reduction for
8 most favorable shorelines (favorability in order of red - orange - yellow - light green - dark green
- light blue - dark blue - indigo) with solid lines indicating condition before TN reductions and
dashed lines indicating projected condition after TN reductions, c) Long term recovery following
40% reduction in total N for 8 most favorable shorelines assuming recovery of sediment percent
organic carbon to reference levels 63
VII
-------
Tables
Table 1. Potential Transplant Suitability Index (PTSI) data listing 3
Table 2. Range of predictors of seagrass habitat for sites at which seagrass is presenter absent
within Narragansett Bay buffer zone. Shoreline event measure is an index assigned to each
10-meter shoreline segment based on distance from the beginning of a shoreline route 13
Table 3. Initial full models evaluated for each seagrass endpoint. Second- and third-order terms
were added to models in later stages only in cases where plots of residual error versus
predictors showed evidence of nonlinearities 14
Table 4. Random effects associated with shorelines for GLMM model 1 predicting seagrass
presence/absence by grid cell. See Figure 3 for map of shoreline codes 22
Table 5. Fixed effects for model 1 predicting seagrass presence/absence by grid cell 23
Table 6. Random effects for GLMM model 2 predicting shoreline presence/absence for seagrass,
not accounting for historic seagrass presence 29
Table 7. Fixed effects for GLMM model 2 predicting shoreline presence/absence for seagrass,
not accounting for historic seagrass presence 29
Table 8. Fixed effects model 3 components predicting shoreline presence and including historic
1999 seagrass predictors 30
Table 9. Fixed effects in model 4 predicting average seagrass occurrence along transects
perpendicular to shoreline, not incorporating historic seagrass presence 34
Table 10. Fixed effects for model 5 predicting average seagrass presence/absence along transects
perpendicular to shoreline, including effects of historic seagrass presence. Model coefficients
are compared for model fit with full data set and with data set minus three outliers 37
Table 11. Fixed effects for model 6 predicting seagrass minimum depth of occurrence by
shoreline distance 40
Table 12. Fixed effects in generalized linear model 8 predicting maximum depth of occurrence
for seagrass at shoreline index i 44
Table 13. Summary of model performance at different scales, based on resubstitution or 10-fold
cross-validation and with or without incorporation of data on historic seagrass coverage 48
Table 14. Potential effect of independent variables across range of predictors. For sediment types,
the range of coefficients is given. Effects are expressed in terms of In (odds ratio). Maximum effect
for Secchi depth is at an intermediate value of the predictor, not the maximum 51
Table 15. Average % silt + clay in McMaster surficial sediment samples from Narragansett Bay
by Shepard's sediment class and estimated reference level of sediment total organic carbon
based on Pelletier etal. (2011) 61
VIII
-------
Acronyms
AIC Aikake's Information Criterion
ALU Aquatic life use
ANN Artificial neural network model
AUC Area Under the Curve
BART Bay Assessment and Response Team
BRT Boosted Regression Tree
C Centigrade
CBP Chesapeake Bay Program
C-CAP Coastal Change Analysis Program
CMAQ Community Multi-Scale Air Quality (Model]
CRM Coastal Relief Model
CVI Coastal vulnerability index
DEM Digital elevation model
DIN Dissolved inorganic nitrogen
EDM Estuary Data Mapper
EPA Environmental Protection Agency
ESI Environmental Sensitivity Index
ESRI Environmental Systems Research Institute, Inc
GAM Generalized additive model
GAMM Generalized additive mixed model
CIS Geographic information system
GLM Generalized linear model
GLMM Generalized linear mixed model
Kd Vertical light attenuation coefficient
LIDAR Light Detection And Ranging
MA Massachusetts
MARS Multi-adaptive regression spline
MCMC Markov Chain Monte Carlo
MGET Marine Geospatial Ecology Tools
MHHW Mean higher high water
MHW Mean high water
MLLW Mean lower low water
MLW Mean low water
MSL Mean sea level
MTL Mean tide level
N Nitrogen
NASA National Aeronautical Space Administration
NAVD North American Vertical Datum
NDBC National Data Buoy Center
NGDC National Geophysical Data Center
NGVD National Geodetic Vertical Datum
NOAA National Oceanic and Atmospheric Administration
IX
-------
ORD Office of Research and Development
P/A Presence/Absence
POC Particulate organic carbon
PSU Practical salinity unit
PTSI Preliminary transplant suitability index
R Software package
RIDEM Rhode Island Department of Environmental Management
RIGIS Rhode Island Geographic Information System
RMSE Root mean square error
ROC Receiver operating characteristic
SAC Spatial autocorrelation
SAV Submerged aquatic vegetation
SD Secchi depth
SS Suspended solids
SST Sea Surface Temperature
STD Standard deviation
TN Total nitrogen
TOC Total organic carbon
TSI Transplant suitability index
US United States
US EPA U.S. Environmental Protection Agency
USGS U.S. Geological Survey
VIF Variance inflation factor
WAVES USGS program to calculate wave statistics
WBID Water body identification code
WEMo Wave Energy Model
WGS World Geodetic System
WWTP Wastewater treatment plant
Zc Light compensation depth
-------
Executive Summary for Coastal Managers
We developed an approach for creating statistical models to predict seagrass
presence/absence at the scale of individual grid cells (10 m x 10 m] as well as a number of
endpoints relative to seagrass presence/absence along transects perpendicular to the
shoreline: presence/absence, relative frequency along transects, and minimum and
maximum depth of occurrence. Unlike many of the statistical models that have been
developed to predict seagrass occurrence, ours take into account the nonrandom
distribution of seagrass (patchiness], nonlinear effects of light availability, salinity (as an
indicator of total N gradients], sediment organic carbon, and parameter interactions.
We tested this modeling approach using data for distribution of seagrass (Zostera marina]
in Narragansett Bay, RI. Models developed were most robust for predictions of shoreline
occurrence rather than at the 10 m x 10 m grid scale. Based on model results, multiple
factors affect seagrass success. The minimum depth of occurrence is influenced by both
wave energy (wave mixing depth] as well as sediment particle size, with finer sediments
more susceptible to disturbance by wave action. The maximum depth of occurrence is
influenced by both transparency and sediment organic carbon, an indicator of past
eutrophication. Depth limits decrease as sediment organic carbon increases due to
increased energetic demands for seagrass to counteract effects of increased toxicity of
anaerobic sediments (e.g., sulfide toxicity]. Our models allowed us to distinguish multiple
modes of action for nitrogen effects on seagrass distribution: 1] shading by phytoplankton
affecting Secchi depths, 2] shading and/or competition by periphyton and macroalgae
mediated by nitrogen concentrations, and 3] effects of sediment organic carbon on
minimum light requirements. Incorporation of historic distribution of seagrass patches did
not improve model predictions, suggesting that patches may exist in a state of dynamic
equilibrium with a lag time for recovery of disturbed sites.
Our models detected significant differences in the probability of seagrass occurrence by
shoreline even after we factored out the effect of site characteristics. This could be
explained by hysteresis effects related to tidal currents. Tidal currents in Narragansett Bay
are strong enough to resuspend fine sediments, thus limiting establishment of new
seagrass patches, but not strong enough to damage established seagrass beds. This
suggests that better monitoring is needed for assessment of turbidity plumes near the
bottom of the water column, not just water column transparency. If hysteresis exists, this
has implications for recovery strategies as well. Specific restoration measures such as co-
restoration of shellfish beds (to reduce suspended sediments and effects of tidal currents
and wave action] or use of existing or constructed coastal barriers to limit effects of wave
action and tides might improve probabilities of initial colonization success and the
initiation of positive feedback effects.
Our predictive models have multiple potential applications: identification of aquatic life use
zones for setting nutrient criteria for areas of potential seagrass habitat, prioritization of
areas and strategies for seagrass restoration, and projection of potential benefits of
management actions. We applied our model to predict the potential recovery of seagrass
given a 40% decrease in total N loading from wastewater treatment plants and
atmospheric deposition (assuming an equivalent reduction in water column
XI
-------
concentrations}. Based on the current model, the colonized area for all shorelines
combined following a 40% reduction in TN loads (and concentration] would increase from
12% of area in the 0 to 5 meter depth zone to about 63% of area in the short term and
slightly more in the long term (as sediment organic carbon levels recover}. Adaptive
management will need to take into account different projections for short-term versus
long-term recovery due to the multi-decadal persistence of organic carbon in sediments
and effects on minimum light requirements for seagrass.
Abstract
Restoration of ecosystem services provided by seagrass habitats in estuaries requires a
firm understanding of the modes of action of multiple interacting stressors including
nutrients, climate change, coastal land-use change, and habitat modification. Often,
managers have used the reported historic depth limits of seagrass to project the future
distribution of seagrass in response to nitrogen load reductions. In general, these
predictions are based on empirical or modeled estimates of the influence of phytoplankton
production in the water column on the light environment, and do not account for the
interaction of multiple factors. We explored the application of generalized linear mixed
models ( GLMMs} and generalized additive mixed models (GAMMs) to describe the simple
and interactive effects of environmental factors on the distribution of a common seagrass,
Zostera marina, in Narragansett Bay, Rhode Island. We used a random shoreline effect to
account for "founder" (random colonization or extinction] effects. We provide several
strategies to overcome three challenges in developing empirical species distribution
models to describe and predict seagrass distribution in estuaries: the fine-scale patchiness
of seagrass distributions with attendant problems of spatial autocorrelation; the large
areas of interest for model development and application entailing significant memory
demands for modeling; and the potential co-variance of multiple interacting factors
affecting seagrass. We developed a spatial framework describing the coordinates of spatial
autocorrelation in estuarine systems, with the main axis parallel to the shoreline and a
secondary axis perpendicular to the shoreline. We demonstrated an approach to
incorporate a term for residual autocorrelation in GLMMs first introduced by Crase (Crase,
B., Liedloff A.C., and B.A. Wintle. 2012}. To account for anisotropy in the system, we
calculated zonal averages of residual errors within rectangular boxes oriented parallel to
the shoreline along the longer main axis. We successfully dealt with covariance of
influential factors by centering variables, by using multiple strategies to describe the
interaction of the light environment and wave energy with depth, and by excluding
correlated variables where necessary. We predicted seagrass distribution at the scale of
10-meter grid cells, as presence/absence or average presence/absence associated with
shoreline locations spaced at 10-meter intervals, and minimum or maximum depth of
distributions at those locations. Prediction of seagrass absolute or average
presence/absence at shoreline locations was very robust, with area-under-the-curve (AUC]
values associated with Receiver Operating Characteristic (ROC] curves of 0.95 - 0.98
following 10-fold cross-validation of models. Random shoreline effects varied over several
orders of magnitude, probably tied to the distribution of tidal currents. Tidal currents are
XII
-------
weak enough to allow persistence of existing seagrass beds, but strong enough to interfere
with successful recolonization through resuspension of fine sediments. For the model
predicting seagrass presence/absence at the grid cell scale, the most influential predictor is
Secchi depth, followed by (in order}: shoreline isolation, sediment percent total organic
carbon, sediment type, and salinity. The least influential variable is water depth greater
than average wave mixing depth. For the model predicting presence of seagrass at
shoreline locations, the most influential predictor is sediment type, followed by sediment
percent total organic carbon (at low Secchi depth], then salinity (as an indicator of
downstream gradients in water column total nitrogen}. As demonstrated in other recent
studies, sediment total organic carbon interacts with light availability by increasing energy
requirements and the light compensation point for seagrass. For all shorelines combined,
our model predicts that following a 40% reduction in TN loads (and concentration] the
colonized area would increase from 12% of area in the 0 to 5 meter depth zone to about
63% of area in the short term and slightly more over subsequent decades as sediment
organic carbon recovers. Finally, we provide data sources for application of this approach
to other U.S. estuaries, with much of the data available through EPA's Estuary Data Mapper
application (www.epa.gov/edm}.
Keywords: eelgrass; seagrass; Zostera marina; estuary; Narragansett Bay; species
distribution model; generalized linear mixed model; spatial autocorrelation
XIII
-------
Chapter 1. Introduction
1.1 Purpose
Seagrasses are essential in providing valuable ecosystem services but are in decline
globally (Orth etal. 2006}. Compton etal. (2011} identified loss of seagrass habitat as one
of the most costly impacts related to nitrogen loading based on the relationship between
loss of submerged aquatic vegetation and fishery declines in estuaries. Restoration of
ecosystem services provided by seagrass habitats in estuaries requires a firm
understanding of the modes of action of multiple interacting stressors including nutrients,
climate change, coastal land-use change, and habitat modification. Managers often use the
reported historic depth limits of seagrass to project the future distribution of seagrass in
response to nitrogen load reductions based on empirical or modeled estimates of the
influence of phytoplankton production in the water column on the light environment
(Dennison et al. 1993}. However, this approach does not account for the interaction of
multiple factors, so more comprehensive models are needed.
Recent discussions among Environmental Protection Agency (EPA} Office of Water, EPA
Regions, EPA Office of Research and Development (ORD}, and states have highlighted the
need for states to refine aquatic life use (ALU} definitions. The Chesapeake Bay Program
(CBP} has provided one example of tailoring ALUs to reflect different expectations and
habitat support functions for different zones within an estuary (US EPA 2003}. The
Maryland Coastal Bays National Estuary Program has used seagrass potential habitat as a
target for seagrass distribution and as a normalizing factor to describe seagrass coverage
by estuarine segment (Wazniakand Hall 2005}. Their mapping of potential habitat
considers only two factors: depth and % silt/clay in sediment. An extension and refinement
of the CBP conceptual model for ALU definitions by region and habitat zone for different
estuary types could help foster state efforts to refine ALUs. Defining appropriate ALUs
provides the foundation for setting water quality criteria (including nutrient criteria} by
establishing targets or expectations for ecosystem condition in the absence of pollution.
Support for submerged aquatic vegetation (SAV} defines one of the potential habitat uses in
the CBP model for aquatic life use (US EPA 2003}. To extend this model to estuaries in
other regions, we need to define habitat constraints (essentially a habitat suitability index}
for seagrass species dominant over different geographic ranges along the US coast using
readily available data. The regulatory objective of interest is to provide a scientific
framework to support nutrient criteria and restoration plans by defining spatially-explicit
targets for expected condition (seagrass habitat presence/absence} in estuaries. Many
water quality criteria are based on an initial definition of the use (including ALU} that
should be supported. Once areas suitable for seagrass habitat are identified, the state
agencies responsible for setting criteria can determine appropriate targets for seagrass
coverage within a system.
-------
The purpose of this report is to:
1) assess predictive modeling approaches for seagrass presence-absence,
2) discriminate the effects of nutrient enrichment from effects of other stressors
and co-factors on seagrass occurrence to foster improved management of
estuarine systems, and
3) provide users with easy access to input data for predictive models of seagrass
habitat on regional and national scales.
Our modelling approach can be applied to other estuaries, although specific data sources,
variables included, and best models selected may differ across systems. We illustrate the
approach outlined for predictive seagrass habitat models with a case study using data for
the Narragansett Bay estuary in Rhode Island. The resultant models include both variables
related to different mechanisms for nutrient effects on seagrass as well as other cofactors
and stressors which limit the distribution of seagrass. We show how the model can be
applied to assess the potential improvement in seagrass coverage following nutrient load
reductions. We designed the approach to be generic and adaptable to other estuaries. To
facilitate this process, we provide information on data sources and examples of R
programming code in the appendices to this report.
1.2 Conceptual Model
Some of the main factors affecting seagrass presence/absence are water depth, light
availability, temperature, salinity, energy regime, substrate type, sediment sulfide content,
and macro-algal coverage (Burkholder et al. 2007, de Boer 2007, Koch et al. 2007, Lee et al.
2007, Ralph et al. 2007, Touchette 2007}. Many of these factors can be readily monitored in
the field or calculated from measured parameters, mapped, and recorded as digital
geographic information system (CIS] coverages, while others (sulfide concentration,
macro-algal coverage] are not generally available from monitoring programs. Site-specific
disturbances (wasting disease, wetland dredge/fill activities, boat anchors, grazing,
hurricanes] can also limit the distribution of seagrass species but are not predictable, and
so must be accounted for on a case-by-case basis (Short et al. 1987,1988, Neckles et al.
2005, Rivers and Short 2007, Oakley et al. 2013]. We focus our habitat suitability model on
factors for which georeferenced data are readily available, and for which optima or
thresholds for seagrass growth and survival are available from the literature: wave energy,
light availability at the sediment surface, substrate type (particle size, organic matter
content], salinity, and temperature.
1.3 Existing Approaches
1.3.1 Preliminary Transplant Suitability Index:
Short and Burdick (2005] developed a software application based on their report on a
preliminary transplant suitability index (PTSI] for screening potential seagrass transplant
sites. This modeling software takes data availability into account, and requires some on-
site field measurements to assess light and bioturbation. The modeling process was
performed in three stages: the PTSI was developed using readily available data; field data
were collected to assess light conditions and bioturbation; and then the final Transplant
Suitability Index (TSI] score was developed. The PTSI development process uses a set of
-------
parameters (Table 1} and is calculated by multiplying the ratings for each parameter. The
higher the score, the better the likelihood that seagrass will survive transplantation at that
site. The TSI does not include sites that were rated zero. Sites close to existing seagrass
beds receive lower scores not because they have a low likelihood of success but because
transplantation to these sites will not increase the geographic extent of existing seagrass.
Table 1. Potential Transplant Suitability Index (PTSI) data listing.
Parameter
Historical eelgrass
distribution
Current eelgrass distribution
Proximity to natural eelgrass
bed
Sediment type
Wave exposure (calculation]
Water depth
Water quality
PTSI rating
1: previously unvegetated
2: previously vegetated
0: currently vegetated
1: currently unvegetated
0: < 100m
1: > 100m
0: rock or cobble
1: > 70% silt/clay
2: cobble free with < 70%
silt/clay
0: > mean + 2 STD
1: < mean + 2 STD
0: shallow
1: shallow edge of bed
2: average of bed
1: deep edged bed
0: poor
1: fair
2: good
Reference
Fonsecaetal. (1998}
Orthetal. (1994}
Kenworthy & Fonseca
(1977}
Short etal. (1987, 1993}
Koppetal. (1994},
Murphey & Fonseca
(1995}, Fonsecaetal.
(1998}
Short etal. (1993}
Batiuk etal. (1992},
Dennison etal. (1993},
Costa etal. (1996}
1.3.2 Predictive modeling approaches
Predictive modeling approaches for seagrass habitat range from simple empirical
approaches to more detailed mechanistic modeling approaches. The simplest approach
involves predicting presence/absence of seagrass using generalized linear models (GLM;
van der Heide et al. 2009} or generalized additive models (GAM} (Downie et al. 2013}. The
latter allow nonlinear effects to be incorporated into models. One variant of these models
-------
has been proposed to predict potential habitat based on presence only, e.g., using the
Maxent software (Downie et al. 2013}. More complex statistical approaches include
Bayesian models which can incorporate prior knowledge of parameter distributions
(March et al. 2013], and boosted regression trees (BRTs), which can handle both nonlinear
effects and parameter interactions for a large number of predictive variables (Crase et al.
2012}. Some investigators have incorporated not only the potential effects of the physical
environment on individual plants, but also emergent properties influencing the accelerated
growth of seagrass patches (Kendrick et al. 2005}. The most complex process-based
modeling approaches incorporate specific mechanisms underlying growth, loss rates, and
interactions between seagrass and their physico-chemical environment. These processes
can be incorporated into dynamic cellular automata models that can mimic the spread of
seagrass patches with physical and biological feedback loops (Wortman et al. 1997,
Fonseca et al. 2000b, Fonseca et al. 2004, Schonert and Milbradt 2005}.
1.4 Study Site (Narragansett Bay, RI)
Narragansett Bay is one of the largest estuaries in southern New England, running from
north to south along the state mid-line. The upper eastern portion of the estuary, including
part of Mount Hope Bay, is in Massachusetts. The Bay's surface area is approximately 380
km2 with 618 km of tidal shore (Figure 1}. Much of the Bay is shallow and well mixed but
depth varies, and descends to 40 m in some of the deepest channels.
Narragansett Bay's eutrophication problems began in the early to mid 1900s when
increased development and sewering led to increased nitrogen loading to the Bay (Nixon
and Pilson 1983}. More recent data suggest, however, that the bay's nutrient
concentrations are declining (Oviatt et al. 1995}. The bay receives much of its nitrogen
inputs from wastewater treatment plants near the head of the Bay, and nutrient
concentrations exhibit a strong north to south gradient (Oczkowski et al. 2008}. In general,
Narragansett Bay remains well-mixed vertically under normal climatic conditions, and has
a mean hydraulic residence time between 10 and 40 days depending on the time of year
and on location in the Bay (Howarth 1988, Pilson 1985}. Significant runoff events,
particularly during neap tides, can induce the onset of stratification in the Bay by
strengthening vertical salinity gradients.
Seagrasses were widespread in Narragansett Bay throughout the 1800s and early 1900s
(Kopp et al., 1994}. During the early to mid-1900s a precipitous decline in seagrasses was
observed globally and in the New England Region. Many research studies suggest that
decline may have been due to increases in environmental contaminants related to
eutrophication, herbicide use, and increased turbidity. Reported instances of widespread
infection (wasting disease} from Labryinthula zosterae (slime mold} in Narragansett Bay
may have also contributed to seagrass decline (Short et al. 1987,1988}.
-------
,f™' TAUNTON
fi
NARRAGANSETT
Legend
WWTF Locations
R.I. Streams/ Rivers
1
W-<^>E 0 1.5 3
12 Kilometers
Figure 1. Map of Narragansett Bay and Rhode Island with color-coded bathymetry. WWTF = Wastewater treatment
facility.
-------
Chapter 2. Methods
2.1 Data Sources
We identified potential predictors of seagrass habitat based on information in the literature
regarding previous modeling efforts, as well as the public availability of data (Appendix A}.
Detailed information on how to access similar data for other systems is provided in
Appendix B. We selected data to be contemporaneous or within + 5 years of reported
seagrass field sampling whenever possible, focusing on data from the primary growing
season, May through October. We reviewed data and metadata from each source for
completeness and applied further geoprocessing if needed. We converted each data source
to a grid with 10 x 10 meter cells, clipped the grid to the zero to five meter depth zone in
Narragansett Bay. This depth zone represents the maximum potential extent of seagrass if
no optically-active constituents are present in the water column (i.e., no light extinction
due to chlorophyll, turbidity, or dissolved organic matter}. Model points were sampled
from the centroids of these grid cells within the 0 to 5 meter depth band (Figure 2a, b].
This yielded 18,612 grid cells with seagrass present and 2,725,719 grid cells with seagrass
absent.
2.2 Data Pre-processing
2.2.1 Salinity
We derived salinity data for Narragansett Bay from a RI Department of Environmental
Management (RIDEM] dataset collected by the Bay Assessment and Response Team
(BART] (www.ri.dem.gov/bart/netdata.htm}. There are 12 BART sampling buoys that
traverse the Bay from north to south. These data included salinity values collected at 15-
minute intervals from a variety of stations between 2003 and 2012. Weekly salinity data
from 2006 through 2012, were downloaded as excel spreadsheets, compiled, and averaged.
We imported data into ArcGIS as points and created Thiessen polygons from original point
data to fill in spatial data gaps. These data were sufficient to capture the gradient of
salinity from north to south but not detailed enough to represent variation at the scale of
smaller subembayments.
2.2.2 Temperature
We downloaded daily sea surface temperature (SST] data from NASA's multi-scale ultra-
high resolution sea surface temperature remote sensing product (MUR SST; ftp://podaac-
ftp.jpl.nasa.gov/allData/ghrsst/data/L4/GLOB/JPL/MUR/] using the EPA's Estuary Data
Mapper (EDM] interface (www.epa.gov/edm], and calculated weekly mean temperature
values for the years 2003 - 2012. In areas of the estuaries near the shoreline that were not
covered by the SST data, we applied a Euclidean allocation method in ArcMap to fill in areas
without data. Euclidean allocation was used rather than focal statistics to fill in data gaps
because Euclidean allocation uses the closest cell data to the analysis point rather than
using a large rectangle to calculate a mean value (as is the case with focal statistics]. The
resultant grid is a more realistic representation of sea surface temperature.
2.2.3 Sediment Type (Grain Size and Percent Total Organic Carbon)
We obtained sediment particle size classes as ArcGIS shapefiles from the RIDEM
Narragansett Bay Estuary Program that were based on McMaster's collection of 493
-------
surficial samples (I960] (http://www.narrbay.org/biologicaLdata.htm}. Sixteen sediment
classes were available, including two missing data classes; 'not mapped' and 'not sampled'.
Size classes included; gravel, sandy gravel, gravel sand silt, sand, gravelly sand, silty sand,
silt, sandy silt, clay silt, sand silt clay, gravel silt clay and rock. No seagrass occurred in
association with four of the original 14 sediment classes mapped; these four classes were
relatively rare across the shallow water zone of Narragansett Bay.
a)
b)
'
Figure 2a. Map of zone delineating boundaries of seagrass presence/absence grid (0 to 5 meter depth) in Narragansett Bay
and b) close-up of 10-meter grid cells with seagrass presence (red) and absence (black). In Figure 2b, the shoreline is to
the south of the grids.
We dropped the Rock class from the region of interest for modeling and combined the
Sandy Gravel and Gravel-Sand-Silt classes (from which seagrass was entirely absent] with
the Gravel class. We combined the Silt class (from which seagrass was absent] with the
Sandy Silt class. Thus up to eight classes were used in the analysis.
-------
We estimated total organic carbon was estimated across the Bay based on 119 surficial
sediment grabs collected by the US EPA National Coastal Assessment
(http://oaspub.epa.gov/coastal/coast.search}. We interpolated values to create a complete
grid within the shallow-water zone using inverse distance weighting in ArcMap 10.1
(©ESRI, Redlands, CA}.
2.2.4 Seagrass (Current/Historical)
We obtained data layers for recent and historic eelgrass coverages for Narragansett Bay
from Rhode Island Geographic Information System (RIGIS}(2000a, b, 2013} as ArcGIS
shapefiles. Both the 1999 (RIGIS 2000, ab} and the 2006 RIGIS (Currentseagrass} data
were collected by the state of RI using NOAA Coastal Change Analysis Program (C-CAP}
protocols (http://coast.noaa.gov/digitalcoast/_/pdf/ccap-products.pdf), so these datasets
should be comparable. The true color imagery was recorded and analyzed at a scale of
1:12000 for both. The 1999 data were collected on July 6th and the 2006 data were
collected on August 6th.
If eelgrass patch locations are stable from year to year, and patches progressively shrink or
expand, then historic eelgrass locations should be a good predictor of current
presence/absence. We intersected model points with historic seagrass coverage (1996-
1998} to yield a presence/absence dummy variable (egPA99}. To create a composite of
patch boundaries, we converted seagrass polygons to lines and merged these with historic
seagrass line coverages (representing seagrass patches less than 40 feet in width}. Then
we calculated distance to nearest historic seagrass patch boundary for each model point
(DistToEG99}. We also calculated the area of each nearest historic seagrass patch (AREA}
as a potential explanatory variable.
We captured the effect of more distant historic events (wasting disease incidence in the
1930s and subsequent recolonization of seagrass up through the 1960s} in the models with
a random "shoreline" effect. Because seagrass beds expand predominantly through
vegetative growth rather than through reproduction by seeds, the growth of seagrass
patches is most likely to occur along shorelines, areas of contiguous potential habitat. Each
of 18 contiguous shorelines within Narragansett Bay was assigned a code, considering
different potential current delivery vectors (Figure 3}. We assigned a SHORLIN code of-99
to minor shorelines associated with small islands not connected to the mainland by
suitable habitat (0 to 5 meters depth} as well as a dummy variable for ISOLATED status.
2.2.5 Transparency
We described transparency across the Bay using three data sources. We combined Secchi
depths measured by the Narragansett Bay Commission between 2008 and 2012
(http://snapshot.narrabay.com/app/MonitoringInitiatives/WaterClarity} with
measurements collected by the US EPA's National Coastal Assessment program between
2000 and 2006 (http://oaspub.epa.gov/coastal/coast.search}. We averaged values by the
Water Body Identification Code (WBID} estuarine segments used for assessment and listing
by RI DEM. To avoid biasing segment averages, we removed Secchi depths greater than the
maximum depth at the site of measurement from the records before averaging. We filled in
gaps in transparency data along the southern shore of Conamicut Island (Sakonnett and
Newport Bays} using Kd estimates from offshore MODIS satellite imagery. (Conamicut
-------
Island is the large island near the mouth of the western arm of Narragansett Bay, Shoreline
code 14.} We downloaded the latter data using the EDM tool and averaged over the
growing season (May - October] for the years 2008-2012. We extracted average grid cell
values for locations greater than 30 meters in water depth, and averaged values over a
swath of offshore cells parallel to the south shore of Conamicutt Island1. We estimated
Secchi depth from Kd values using an empirical relationship developed by Batuik et al.
(2000}.
11 Miles
Figure 3. Shoreline code assignment for Narragansett Bay. A value of -99 was assigned to unused segments and isolated
shorelines.
1 Algorithms for Kd are not suitable for shallower water due to interference from bottom sediments.
-------
2.2.6 Wave Exposure Data - WEMo and WAVES models
Assessment of estuarine wave energy is an important component of seagrass modeling.
Wave energy relative to fetch and current can affect seagrass habitats in a number of ways,
mainly by reducing the ability to establish and maintain effective beds, changing sediment
grain size type and or loss, and physically damaging delicate sheaths and leaves. Other
ecological processes (e.g., TOG accumulation, bioturbation etc.] are also linked to wave
energy and velocities in seagrass habitats. Particle deposition and sediment resuspension
are affected by wave energy, but these processes can be attenuated to some degree by
established seagrass beds (Garcia etal. 1999}.
A variety of tools can be utilized to spatially model fetch and wave energy inputs to
estuarine shorelines (Howes etal. 1999, Fonseca and Malhotra 2010, Rohweder etal.
2012}. We used two tools to spatially model fetch and wave energy inputs, the National
Oceanographic and Atmospheric Administration (NOAA} Wave Exposure Model (WEMo;
http://products.coastalscience.noaa.gov/wemo/} for calculation of wave energy, and the
United States Geological Survey (USGS} WAVES tool for calculation of mixing depth. Each
of these tools has varying input requirements relative to data and formatting, and all
require additional pre-processing of input data through modeling applications. WEMo
requires a variety of data inputs: closest weather station data (average wind speed and
direction}, bathymetry data (ESRI ArcGIS GRID format}, shoreline edge (ESRI ArcGIS GRID
format}, and sampling points (generated by the user}. We performed WEMo modeling
using preset (default} values/conditions as specified in Fonseca and Malhotra (2010}. No
other WEMo modeling parameters were adjusted except for fetch interrogation distance,
which was set to 5000 meters.
We calculated average and maximum wave mixing depths for the 2006-2007 growing
seasons using the USGS WAVES extension to ArcMap 10.1 (Rohweder et al. 2012}.
Chambers (1987} has predicted minimum depth of occurrence of submerged aquatic
vegetation in lakes based on calculation of mixing depths from wave action, as one half of
the wave length (Chambers 1987}. We calculated fetch and wave characteristics using time
series of wind speed and direction from the nearest National Data Buoy Center (NDBC}
station (Station PTCR1 - 8452951 - Potter Cove, Prudence Island, RI} merged with
topobathymetry grids from NOAA's Coastal Relief Model (NOAA National Geophysical Data
Center, NGDC Coastal Relief Model, Volume 1-8,
http://www.ngdc.noaa.gov/mgg/coastal/coastal.html}.
2.2.7 Distance to physical disturbance source: hardened shorelines and marinas
We evaluated two different indicators of physical disturbance related to anthropogenic
activities as predictors in regression models: distance to hardened shorelines and distance
to marinas. Structural shoreline hardening in Rhode Island includes the use of rock
revetments, bulkheads and other types of walls or groins. Much of this hardening took
place before coastal regulations existed. Many of the effects of shoreline hardening are
very localized, with wave action reflecting off of hard structures causing scour (Shipman
2010}. We obtained a CIS layer with locations of hardened shorelines from RIGIS (2003}.
Other potential sources of information on shoreline hardening (not used here} available for
models in other estuaries include the Environmental Sensitivity Index (ESI}
10
-------
(http://response.restoration.noaa.gov/maps-and-spatial-data/download-esi-maps-and-
gis-data.html] and high resolution imagery with coastal LIDAR.
Although mooring beds are permitted throughout Narragansett Bay, there are no common
CIS data layers available to document their location. As an indicator of the potential for
damage from boat anchors, we calculated the distance to nearest marina based on a data
layer obtained from RIGIS (1996} using the NEAR function in ArcGIS.
2.2.8 Unsewered residential development on high infiltration soils
Groundwater inputs to estuaries (with attendant nutrient loads] can be significant, but are
generally poorly quantified. If groundwater nutrient inputs are taken up within seagrass
beds, they may not be reflected in overlying phytoplankton concentrations. We used the
density of unsewered residential development occurring on high infiltration soils in coastal
catchments as an indicator of potential groundwater nitrogen inputs to seagrass habitats.
The Narragansett Bay Sustainability Pilot Appendix A describes the derivation of data for
this indicator (http://www2.epa.gov/sites/production/files/2013-12/documents/nbsp-
phase-i-report-appendices.pdf) (data provided by Industrial Economics, Inc.}. Grid
centroid points were spatially "joined" to the polygon coverage to add this variable to point
attributes.
2.2.9 Canada goose grazer density
Canada geese have been identified as a potential impact to seagrass growth and survival in
Narragansett Bay (Rivers and Short 2007}. The density of geese and their impact will vary
by location and seagrass type. We estimated Canada goose density by zone in Narragansett
Bay based on winter waterfowl surveys conducted in 2004 - 2006
(http://www.nbnerr.org/waterfowl.htm}. Shapefiles of survey zones were provided by
Rick McKinney, US EPA Atlantic Ecology Division, Narragansett, RI.
2.3 Statistical Model Development
2.3.1 Strategies to limit memory requirements
Due to the large geographic area covered and fine spatial resolution of seagrass
presence/absence grids, we applied several strategies to limit memory requirements for
statistical analyses. First, following the best practices outlined for species distribution
models by Maggini et al. (2010}, we restricted analyses to the geographic range and range
of predictor variables associated with points at which seagrass was present in 2006 (Table
2}. Second, we limited analyses to the 0-5 meter depth zone. Third, we generated models
at multiple scales: a} predicting seagrass presence/absence at individual grid cells (grid
presence/absence P/A; n = 518,890}, b} predicting P/A at any point along transects
perpendicular to the shoreline (shoreline segment P/A; n = 19,204}, c} predicting relative
frequency of occurrence along transects perpendicular to the shoreline (shoreline segment
frequency; n =19,204}, and d} predicting minimum and maximum depth of occurrence
along transects perpendicular to shoreline where seagrass was present (seagrass minimum
or maximum depth; n = 2,749; Figure 4}.
11
-------
2.3.2 Modeling approaches at different scales
Depending on the scale of dependent variables defined, we applied different model
approaches, different initial sets of independent variables with associated interaction
terms and R packages (Table 3, Figure Cl in Appendix C]. We evaluated GLM and GAM for
minimum and maximum depth endpoints. For all other seagrass endpoints, we developed
mixed models (including both random and fixed effects, see below}. For mixed models, we
evaluated two approaches, the first a general linear mixed model which incorporated
random effects related to potential shoreline-specific colonization or disease incidence
using the glmmPQL function in the R MASS package (http://stat.ethz.ch/R-manual/R-
patched/library/MASS/html/glmmPQL.html]. Second, where memory demands were not
too large to prevent application, we also analyzed general additive mixed models using the
function gamm in the mgcv R package to account for potential nonlinearities in response
(https://stat.ethz.ch/R-manual/R-devel/library/mgcv/html/gamm.html}.
We modelled "shoreline" as a random effect because seagrass coverage varies significantly
among Narragansett Bay shorelines. The "shoreline effect" can be interpreted as the
combined result of random colonization (i.e., of wasting disease organisms] and re-
colonization (seagrass recovery] effects. We evaluated two different types of mixed effect
models with respect to their performance in predicting seagrass presence/absence in
Presence/absence
Grid cell
Shoreline
Relative abundance
Shoreline
Figure 4. Seagrass endpoints defined at different scales: presence/absence (P/A) at individual grid cell centroids, P/A within
transect perpendicular to shoreline (assessed after points were snapped to shoreline), relative frequency of points along
perpendicular transect, minimum depth (z) of occurrence along transect, and maximum depth of occurrence along transect.
12
-------
Table 2. Range of predictors of seagrass habitat for sites at which seagrass is present or absent within Narragansett Bay
buffer zone. Shoreline event measure is an index assigned to each 10-meter shoreline segment based on distance from the
beginning of a shoreline route.
Seagrass present
(n = 14527}
Variable
mm
max
Seagrass absent
(n = 1309011}
Min max
Depth
Salinity (PSU}
Avg Secchi depth (m}
Min Secchi depth (m}
Avg SD - water depth (m}
Min SD - water depth (m}
Sediment total organic carbon (%}
Isolated shoreline
Wind wave energy Qoules/m}
Avg wave mixing depth: water depth
Max wave mixing depth: water depth
Temperature (deg C}
Shoreline event measure (m}
0.0
26.6
1.3
0.7
-2.7
-3.7
0.0
0
0
0.002
0.002
14.7
125430
5.1
28.7
3.8
3.7
3.7
3.7
3.5
1
81724
11.2
145.4
15.0
3094040
0.0
11.3
0.8
0.5
-4.2
-5.0
0.0
0
0
0.002
0.002
14.5
115340
6.6
28.7
3.8
3.7
3.8
3.7
7.7
1
91021
12.4
206.6
15.0
3957430
Narragansett Bay: general linear mixed models (GLMM} and general additive mixed models
(GAMM}. We analyzed initial model residuals for evidence of spatial autocorrelation (see
below}.
2.3.3 Selection of best-fitting models
We applied different strategies to identify best-fitting models for simple GLMs and GAMs as
compared to mixed models. For GLMs and GAMs we refined original models (Table 3}
using the step function in R. We compared model fits for GLM or GAM using AIC values and
prediction errors (Burnham and Anderson 2002, Zuur et al. 2009}. For GLMMs, this
approach was not possible because AIC values are not provided by the R packages we used.
We fit models in a manner analogous to backward stepwise regressions, sequentially
eliminating variables from initial full models (Table 3} based on lack of significance (i.e., p >
0.01}. The glmmPQL package yields approximate p-values which should be used with
caution when they are marginal (i.e., near p = 0.05}, so p-values were interpreted
conservatively using a p-value of 0.01. We evaluated final models using a 10-fold cross-
validation procedure to test for robustness of results and as a check against overfitting.
Sampling with replacement was conducted separately within each shoreline class to
provide proportional representation.
13
-------
Table 3. Initial full models evaluated for each seagrass endpoint. Second- and third-order terms were added to models in
later stages only in cases where plots of residual error versus predictors showed evidence of nonlinearities.
Seagrass model dependent variable
Model type
R package(s)
Initial set of independent variables included in
Variable
Fshorlin
cSAL
cTEMPER
fSEDn
cWIND
cPTTOC
csecchiminMax
cSDavggtrZ
cSDmingtrZ
fZgtrMXZ
cZgtrMXZ
halflen0708gsa
Max
fISOLATED
Definition
Shoreline code
Centered growing
season salinity
Centered growing
season average water
temperature
Sediment type (n
represents # classes
after lumping)
Wind wave energy
Centered sediment
percent total organic
carbon
Centered Secchi Depth
(max along transect)
Centered growing season
average Secchi depth -
water depth
Centered growing season
minimum Secchi depth -
water depth
Depth greater than wave
mixing depth (0 = FALSE,
1=TRUE)
Centered water depth
greater than wave mixing
depth
Maximum growing season
wave mixing depth
Isolated shoreline (0 =
Shore-line
Shore-line segment Minimum
Grid cell segment relative depth of
P/A P/A2 frequency1 occurrence
GLMM GLMM GLMM GLM, GAM
glmmPQL glmmPQL glmmPQL glm, gam
models prior to backward stepwise selection
Units Fixed or
random
effect
-99,1-18 R x X X
PSU F x x x
DegC F x x x
ItolS F x x x x
Joules/m F x x x x
% F x x x
Meters F
Meters F x x x
Meters F x x x
F x x x
Meters F x x x
Meters F x
F x x x
Maximum
depth of
occurrence
GLM, GAM
glm, gam
x
x
FALSE, 1=TRUE)
1 Average or optimum (minimum or maximum depending on variable) value for transect substituted for point
values
14
-------
Table 3. (Continued)
Seagrass model dependent variable
Shore-line
Shore-line segment Minimum Maximum
Grid cell segment relative depth of depth of
P/A P/A3 frequency1 occurrence occurrence
Model type
R package(s)
GLMM GLMM GLMM
glmmPQL glmmPQL glmmPQL
GLM, GAM GLM, GAM
glm, gam glm, gam
Initial set of independent variables included in models prior to backward stepwise selection
Variable Definition Units Fixed or
random
effect
cDistHdShor Centered distance to Meters
hardened shoreline
cDistMarina Centered distance to Meters
nearest marina
cUSRMARIkm2 Centered unsewered #/km2
residences on high
infiltration soils/catchment
area
x
x
cCG046avkm2
fEGPA99
cAREA
cDistEG99
Interaction
terms
Centered winter 2004-2006 #/km2
Canada goose density
Historic (1999) eelgrass
presence (0 = FALSE,
1 = TRUE)
Centered area of 1999 Meters2
eelgrass patch
Centered distance to edge of Meters
nearest 1999 eelgrass patch
cWINDxfsedn
halflen0708gsaMax x
fSEDn
csecchiminMaxx
cpttocMax
csecchiminMaxx
cUSRMARIkm2Max
cpttocMax x
F x x x
F x4 x x
F x1 x x
F x1 x x
F x x x
F
F
F
F
x
x
x
x
cUSRMARIkm2Max
1 Average or optimum (minimum or maximum depending on variable) value for transect substituted for point
values
4 Models were run both with and without historic seagrass predictors
15
-------
Table 3. (Continued)
Seagrass model dependent variable
Model type
R package(s)
Initial set of independent variables included in
Variable
Definition
csecchiavMax x cpttocMax
xcUSRMARIkm2Max
cSAL x cSDavtomaxZ
cSAL x fZgtrMXZ
cSDavtomaxZ x fZgtrMXZ
cSAL x cptTOC
cSDavtomaxZ x cptTOC
fZgtrMXZ x cptTOC
cSAL x cSDavtomaxZ x
fZgtrMXZ
cSAL x cSDavtomaxZ x
cptTOC
cSAL x fZgtrMXZ x cptTOC
cSDavtomaxZ x fZgtrMXZ x
cptTOC
cSALx cSDavtomaxZ x
fZgtrMXZ x cptTOC
models prior to backward
Units Fixed or
random
effect
F
F
F
F
F
F
F
F
F
F
F
F
Shore-line
Shore-line segment Minimum
Grid cell segment relative depth of
P/A P/A5 frequency1 occurrence
GLMM GLMM GLMM GLM, GAM
glmmPQL glmmPQL glmmPQL glm, gam
stepwise selection
x
XXX
XXX
XXX
XXX
XXX
XXX
XXX
XXX
XXX
XXX
XXX
Maximum
depth of
occurrence
GLM, GAM
glm, gam
1 Average or optimum (minimum or maximum depending on variable) value for transect substituted for point
values
16
-------
2.3.4 Model diagnostic tests
For both simple and mixed models, we used a series of diagnostic tests to check model
assumptions. Methods and example codes are described in detail in Appendix C. We
checked models for independence of predictor variables using correlation coefficients and
variance inflation factors (VIFs), for homogeneity of variance based on visual examination
of residual versus predicted value plots, for linearity of response based on visual
examination of conditional probability plots and of residuals plotted against predictors,
and for spatial autocorrelation (SAC}.
2.3.5 Development of coordinate framework to assess SAC
To assign coordinates to grid cell centroids, we first assigned a shoreline "event" measure
for each observation. We edited the shoreline polyline to facilitate development of a
shoreline "route" in ArcMap. Creation of a route in ArcMap allows one to assign distances
to points along the line relative to the start of the shoreline route, termed "events". We
eliminated areas of channel "braiding" in estuarine headwaters (and thus, divergent
flowpaths}, and "flipped" the direction of line segments in the polyline to ensure that all
segments were oriented in a common direction. We visualized problem segments by
changing the line symbol to an arrow and noting locations with adjacent arrowheads
pointing in opposite directions. We "dissolved" the shoreline to minimize the number of
features, then used the dissolved shoreline to create a shoreline route. We artificially
created equidistant points along the shoreline route by generating an event table in
Microsoft Excel, with MEAS values increasing in 10-meter increments and DISTANCE (from
line] set to zero, then used the event table and shoreline route to create shoreline events in
ArcMap. We saved the temporary shoreline events as permanent point features, then used
the NEAR function in ArcMap to assign each observation to the nearest shoreline point with
an associated shoreline distance (Figure 5}. We provided variogram or correlation
structure functions in R with the MEAS value as an "X-coordinate" and a constant of zero as
the "Y-coordinate" to allow the program to calculate interpoint distances parallel to the
shoreline.
2.3.6 Assessment of SAC ranges
Due to memory constraints with R packages and anisotropy we used two different R
packages to evaluate SAC. Anisotropy exists when spatial autocorrelation varies by
direction. Based on the elongated nature of seagrass patches and the different controls on
spread along the shoreline as compared to along depth gradients, one would expect
anisotropy in spatial autocorrelation of seagrass presence/absence. Vegetative spread of
seagrasses predominates over long range seed dispersal as the mechanism for patch
growth (Marba and Duarte 1998}. Thus, seagrass is more likely to spread in a direction
parallel to shorelines rather than perpendicular to shorelines over areas of deeper water.
To evaluate the range of spatial autocorrelation parallel to shorelines, we used spline
correlograms based on shoreline distance. We evaluated SAC in a direction parallel to
shorelines using nonparametric spline correlograms with the R spline.correlog function in
the ncf package (Bjornstad 2013}. Creation of nonparametric spline correlograms requires
no assumptions about the distribution of errors or the shape of the correlogram plots.
Thus, spline correlograms are more appropriate than standard correlograms for ecological
17
-------
Event measure along route
= x-coordinate
Distance from grid centroid
to shoreline event = y-
coordinate
Shoreline turned into "route" with
associated direction and distance along
shoreline
Seagrass grids converted to points
Points "snapped" to nearest shoreline as
shoreline "event) to determine x
coordinate
Distance from original point to nearest
shoreline determined y coordinate
Shoreline route
Figure 5. Process for assigning x-, and y-coordinates to seagrass grid cell centroids.
datasets which can exhibit a complex pattern of SAC at different scales resulting from
different processes. We created coordinate inputs to the spline.correlog function by
assigning shoreline distance as the x-coordinate and a fixed off-shore distance of zero as
the y-coordinate. The spline function encounters memory limitations when more than
-9000 points are included in the analysis. Therefore, we had to develop the spline
correlograms separately by shoreline, using a random subset of 9000 for each shoreline.
To evaluate the range of autocorrelation along depth gradients perpendicular to the
shoreline we used the CommunityCorrelogram package in R to create Mantel correlograms
(Andrus et al. 2014}. The CommunityCorrelogram package adds functionality over existing
Mantel correlogram functions by allowing directional (anisotropic] restrictions in both the
xy (surface] plane and the z (depth] plane.
2.3.7 Eliminating or incorporating SAC into models
We evaluated two complementary methods to incorporate the effects of spatial
autocorrelation into models. First, we incorporated local SAC into models using the
approach of Crase et al. (2012]. We initially fit models assuming spatial independence of
errors. We then used zonal statistics in ArcMap to estimate a local average of original
model residuals for each grid cell across the local range of SAC, using a moving window of 3
x 3 grid cells. We then incorporated the residual zonal averages into a modified predictive
model to explicitly include the effect of SAC. We again checked adjusted model Pearson
residuals for SAC patterns using the spline correlograms. The output includes confidence
intervals for the spline fit and can be used to generate 95% confidence intervals for the first
zero intercept, or the range at which spatial autocorrelation is no longer observed.
18
-------
Attempts to use the Crase et al. approach in its original form failed because the range of
SAC was much greater than the scale of the 3x3 moving window. Thus we reran the
process using the empirically-derived range of SAC to adjust the size of the moving window
for calculating residual zonal averages. We used the median value of SAC values across
shorelines to calculate zonal averages, which also corresponded to the shoreline with the
greatest abundance of seagrass.
Second, we eliminated some of the spatial autocorrelation in models by reducing
dimensionality. We assigned the associated shoreline event measure (shoreline location]
to each seagrass grid centroid, then calculated an average and maximum presence/absence
value to estimate relative seagrass frequency (along a swath extending perpendicular to
the shoreline] or simply seagrass presence/absence by 10-meter shoreline unit. For each
shoreline event with seagrass present, we also calculated a minimum depth of occurrence
and a maximum depth of occurrence. Predictive models for minimum and maximum
depths were simpler because the random shoreline effect could be dropped, yielding a
fixed effect model.
2.3.8 Strategies for eliminating multicollinearity
Development of GLM, GLMM and GAM models was constrained by the need to exclude
predictor variables that exhibited collinearity (as evidenced by Spearman correlation
coefficients of > 0.7] or multicollinearity (as evidenced by variance inflation factors (VIF] <
10; Zuur et al. 2009]. Many of the variables reported in the literature as potential
influences on seagrass growth and survival are potentially correlated with one another,
making it difficult to determine precise model coefficients. We minimized the potential for
collinearity by centering continuous variables (i.e., subtracting the mean] as per
recommendations of Zuur etal. (2009]. Centering variables also facilitates interpretation
of coefficients in logistic models because the model intercept represents the predicted
value when all independent variables are at their mean levels (and dummy variables equal
zero]. To prevent multi-collinearity, we also dropped sediment particle-size categories
with no seagrass from the model, combined sediment particle-size categories representing
a very small fraction of the model points into new classes, and adjusted the chosen
"reference" level of the sediment type factor so that it represented one of the more
common sediment types.
Effects of both transparency (Secchi depth] and wave energy (mixing depth] have to be
evaluated with respect to water depth. The ratios of Secchi depth to water depth and
mixing depth to water depth are highly correlated because they share a common
denominator. Therefore we created one predictor based on the difference between Secchi
depth and water depth and a second categorical predictor indicating whether or not mixing
depth exceeded the actual water depth. The former variable had a more even distribution
than the corresponding ratio and was less likely to create problems with outliers.
Temperature and salinity were highly correlated, so we dropped temperature from our
predictive models because it was not likely to reach critical levels in the subtidal zone of
this system, and Z. marina does not occur in exposed locations in Narragansett Bay. We
initially evaluated a subset of potential interaction terms in regression models based on
probable mechanisms of action (Table 3]. However, we eventually dropped almost all
19
-------
potential interaction terms from the GLM, GLMM, and GAM models because they had
extremely high VIF values.
2.3.9 Alternative approaches
We were limited in the range of R statistical packages available for mixed effects logistic
models that could also incorporate spatial correlation structures
(http://glmm.wikidot.com/pkg-comparison}. The most commonly used functions for
mixed modeling in R are: MASS::glmmPQL, Ime4::glmer, MCMCglmm::MCMCglmm. Of the
available R packages, glmmPQL is less memory-intensive as it applies a Penalized Quasi-
Likelihood approach rather than producing true Maximum Likelihood Estimates (Venables
and Ripley 2002}. As a result, glmmPQL does not yield an estimate of log-likelihood or
Aikake's Information Criteria (AIC}, values commonly used in comparing model fits.
Although we attempted to fit models using the glmer function in Ime4, models often failed
to converge, yielded "out of memory" errors before convergence, or proved impractical
because model runs required more than 12 hours for a single iteration. Attempts to run the
glmmADMB package in R also failed.
In summary, our approach yielded predictive models for five seagrass endpoints: grid-cell
presence/absence, shoreline presence/absence, shoreline relative frequency, minimum
depth of occurrence, and maximum depth of occurrence. We evaluated glmmPQL or GAMM
models for the first three endpoints, and GLM or GAM models for the latter two. We
incorporated SAC into predictive models using the Crase et al. (2012} approach after
modifying the moving window to account for the actual range of SAC observed in model
residuals. Finally, we minimized the potential for multi-collinearity across independent
variables. See Appendix C, Figure Cl, for a summary of the CIS and R processes necessary
to develop the GLMMs.
20
-------
Chapter 3. Results
For more detailed results, see Appendix D for exploratory analyses and description of
intermediate models.
3.1 Final Models
3.1.1 Seagrass Grid Presence/Absence
The final model fit for seagrass grid cell occupancy was:
fsavcode = cSAL + cSAL2 + cSDavgrtrZ + cSDavgrtrZ2 + cSDavgrtrZ3 + cptTOC
+ cptTOC2 + cptTOC3 + fZgtMXZav + cCG046avkm + cDstoMarin + flSOLATED (1}
+ fSED4 + PResidl4fa,
where fsavcode = seagrass presence/absence code (1/0}
cSAL = centered salinity
cSDavgrtrZ = centered average Secchi depth minus water depth
cptTOC = centered sediment percent total organic carbon
fZgtMXZav = indicator of depth greater than average wave mixing
depth(1/0}
cCG046avkm2 = centered Canada goose density/square kilometer
cDstoMarin = centered distance to nearest marina
flSOLATED = indicator of isolated shoreline (1/0}
fSED4 = sediment particle-size class
SED4_5 = Sand
SED4_7 = Silty sand
SED4_81011 = Silty, Sandy silt, Clay-Silt
SED4_6 = Gravelly sand
SED4_124 = Gravel, Sandy gravel, Gravel-sand-silt
PResidl4fa = Focal average of Pearson residual over zone of 1300 m
shoreline length x 200 m offshore distance
We calculated the residual autocorrelation term as the focal average of Pearson residuals
over a zone of 1300 meters (along shoreline direction} by 200 meters (offshore direction}.
After incorporation of the residual autocorrelation term, the final model showed minimal
spatial autocorrelation (Figure 6a} and much reduced heterogeneity of variance (Figure
6b}. Area under the ROC curve based on ten-fold cross-validation was 0.7144 (Figure 7}.
The model showed a small tendency to underpredict seagrass presence, with a mean
21
-------
residual of -0.1285, and a median value of -0.001545 (interquartile range = -0.0739 to -
0.000043}. The root mean squared error was 0.34.
Colonization rates vary dramatically among shorelines (Table 4], with exponentiated
coefficients (the associated odds ratio] ranging over nine orders of magnitude from 6.1 E-5
(Shoreline 2] to 93091 (Shoreline 17]. The exponentiated intercept for the fixed effects
portion of the mixed model is the odds ratio when all continuous variables are at mean
values (centered value = 0], fZgtMXZav = 0 (depth is less than average wave mixing depth],
fISOLATED = 0 (main shorelines], and the reference sediment class (5 = sand] is 1 (true]
(Table 5]. The corresponding probability is 0.22. The odds ratio for average conditions of
seagrass on sand substrate for isolated shorelines is 0.279 + 0.004 = 0.283. The odds ratio
increases by 1.43 (43%] for water depths greater than the wave mixing depth. The odds
ratio for sediment types 8 and 10 combined (silty and sandy silt] is not significantly
different than for the sand class, but is relatively lower for classes 6 (gravelly sand] or 12
(sand-silt-clay] and higher for classes 7 (silty sand] or 1+2+4 (gravel, sandy gravel, and
gravel-sand-silt], respectively. The conditional odds ratio is 6.65 (or 565% greater] for the
combined gravel classes. Predicted effects for unsewered residential density on high
infiltration soils, Canada goose density, and distance to nearest marina are all opposite in
direction to those expected but the magnitudes of predicted influences are negligible. The
probability of seagrass presence increases at an accelerated rate as salinity increases.
Table 4. Random effects associated with shorelines for GLMM model 1 predicting seagrass presence/absence by grid cell.
See Figure 3 for map of shoreline codes.
Shoreline
-99
2
3
6
7
8
9
10
11
12
13
14
15
16
17
Random effect
2.78
-9.70
-5.67
-9.62
-6.22
-4.14
0.42
3.00
-2.64
1.46
-4.33
0.28
6.65
9.70
11.44
Exp(R.E.)
16.2
6.1E-05
0.003
6.66E-05
0.002
0.016
1.5
20
0.07
4.3
0.01
1.3
771
16398
93091
The odds ratio for seagrass presence is negligible for water depths less than the Secchi
depth, increases above 1.0 when Secchi depth exceeds water depth by one meter, and is
predicted to peak at Secchi depths about 2.5 meters greater than water depth (Figure 8;
note that these odds ratios must be adjusted for the shoreline of interest].
22
-------
Table 5. Fixed effects for model 1 predicting seagrass presence/absence by grid cell.
Std.
(Intercept)
cSDavgrtrZ
cSDavgrtrZ2
cSDavgrtrZ3
cSAL
cSAL2
cptTOC
cptTOC2
cptTOC3
cCG046avkm
cDstoMarin
fSED4124
fSED47
fSED4810
fSED412
fSED46
fZgtMXZav
fISOLATEDl
PResidl4fa
Coeff
0.42
2.61
-0.23
-0.06
1.43
-0.57
-2.44
1.32
-0.32
0.01
-0.001
3.66
0.63
-0.38
-0.83
-1.12
0.53
-10.67
3.17
Error
1.66
0.05
0.03
0.010
0.19
0.16
0.09
0.07
0.05
0.004
0.0001
0.33
0.07
0.15
0.15
0.21
0.14
1.06
0.06
DF
518856
518856
518856
518856
518856
518856
518856
518856
518856
518856
518856
518856
518856
518856
518856
518856
518856
518856
518856
Exp(Coeff)
1.52
13.63
0.79
0.95
4.20
0.57
0.09
3.73
0.73
1.01
1.00
38.77
1.88
0.69
0.43
0.33
1.71
0.00
23.88
t-value
0.25
48.69
-7.11
-5.86
7.60
-3.65
-28.45
18.60
-5.78
2.93
-27.87
10.93
8.76
-2.51
-5.50
-5.25
3.74
-10.06
56.60
p-value
0.8001
<0.0001
<0.0001
<0.0001
<0.0001
0.0003
<0.0001
<0.0001
<0.0001
0.0034
<0.0001
<0.0001
<0.0001
0.0119
<0.0001
<0.0001
0.0002
<0.0001
<0.0001
23
-------
.B o
"to
5 10
Distance (km)
15
b)
300
al
200
son
100
Pea
o
o
L
0.0 0.2 0.4 0.6 O.S 1.0
Predicted seagrass presence in grid cell
Figure 6. Diagnostic plots for model 1 predicting seagrass grid occupancy. Spatial autocorrelation of residuals is virtually
eliminated. Heterogeneity of variance is greatly reduced, a) Spline correlogram of Pearson residuals for Shoreline 14. b) Plot
of Pearson residuals versus predicted value.
a) b)
CD
O
I-J
C\|
d'
J'
.
CO
CN
CD
O
O
1.0 0.8 0.6 0.4 0.2 0.0
Specificity
1.0 O.E
0.6 0.4 0.2
Specificity
0.0
Figure 7. ROC curve for model 1 based on a) initial fit (area under curve = 0.9767) and b) 10-fold cross-validation (area under
the curve is 0.7144).
24
-------
a)
1/1
1/1
£
00
ro
12
£
o
ro
5
O
0)
u
c
§
3
u
u
O
-2
0.7
0.6
0.5
0.4
0.3 ,
^^x
rtJJr
^-^0.1
-1.5 -1 -0.5 0 0.5 1
Centered Salinity (PSU)
b)
£
00
ro „,
S5 8
»- ?.
•2 3
£ 8
3.5
3.0
2.5
2.0
1.5
1.0
0.5
-4-2024
Centered average Secchi depth - water depth (m)
c)
1/1
1/1
ro
ab
ro
0)
M
o
*^~
o
£
M
O
0)
u
c
0)
u
u
O
-2
80
, 60
\
\ 40
\
V 2°
v. 0
-10123
Centered sediment percent total organic carbon
Figure 8. Nonlinear effects of a) centered Salinity (PSU), b) centered average Secchi depth - water depth (m), and c) centered
sediment percent total organic carbon on odds ratio for seagrass occurrence in grid cell. Values are calculated assuming
average values for all co-variables not represented in plot and across all shorelines with the reference sediment type for
model 1.
25
-------
3.1.2 Predictive models for shoreline segment P/A and minimum or maximum depth for
seagrass
3.1.2.1 Shoreline segment P/A models without Serial Autocorrelation
The form of the final model for predicting seagrass presence/absence by shoreline distance
was:
fsavcodeMaxi = csalAvi + csecchimim * cpttocMin + fZgtMXZavMi +
cCG046avknii + cDistToHdSi + cDstoMarim + cwindMim * fSED8i +
PR4wFA1370i,
where i = shoreline event measure
fsavcodeMaxi = Maximum seagrass presence/absence code (0/1}
associated with shoreline event measure i
csalAvi = Average of centered average salinity at i
csecchimim = Minimum of centered Secchi depth at i
cpttocMim = Minimum of centered percent total organic carbon at i
fZgtMXZavMi = Maximum of 0/1 indicator of depths greater than
average wave mixing depth at i
cCG046avkm2i = Minimum of centered Canada goose density
(No./km2)ati
cDistToHdSi = Maximum of centered distance to hardened shoreline
ati
cDstoMarim = Maximum of centered distance to nearest marina at i
cwindMim = Minimum of centered wind relative energy at i
fSED8i = Majority sediment class at i, including
fSED8_6 = Gravelly-sand
fSED8_7 = Silty sand
fSED8_124 = Gravel, Sandy Gravel, and Gravel-Sand-Silt
fSED8_512 = Sand and Sand-Silt-Clay
fSED8_8101113 = Silty, Sandy Silt, Clay-Silt, and Gravel-Silt-
Clay
PR4wFA1370i = Focal average over 1370 meters distance of Pearson
residual for model fit
The dummy variable for "isolated" shorelines, corresponding to small islands, was no
longer significant after spatial autocorrelation was accounted for, and some of the sediment
classes had to be collapsed from the original model to retain significance. After residual
autocorrelation was incorporated into the model, heterogeneity of variance decreased
26
-------
substantially (Figure 9a-c). The original model Area-Under-The-Curve (AUC) value for the
Receiver Operating Characteristic Curve (ROC] was 0.95. Performance of the predictive
model for seagrass shoreline presence was much more robust than performance of the
model predicting seagrass presence at the grid cell scale. Ten-fold cross-validation of the
model yielded a mean residual of -0.111, a Root Mean Squared Error (RMSE] of 0.29
(interquartile range = 0.014 - 0.155], and an ROC value of 0.9547 (see Figure 10}. Again,
the odds ratio varied by several orders of magnitude across different shorelines (Table 6}.
Higher order effects were not retained in the best GLMM model predicting shoreline
presence. In this case, the odds ratio was greatest for the reference sediment class (6 =
gravelly sand] and least for sediment classes 1 + 2 + 4 (gravel-dominated}. There was an
additional interaction term involving Secchi depth and sediment percent total organic
carbon, which tends to decrease the positive effects of transparency at high sediment
organic carbon levels (Table 7}.
27
-------
b)
ra o
OJ TrH
Q- '
O
CM'
w o
OJ
ra
01 o
Q- ,-H
O
IN
0.0
0.2
0.4 0.6
Predicted
TJ
'a
oe
o
O
(N
0.8
1.0
7 12 124 810
Sediment Class
-1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5
Centered percent total organic carbon
Figure 9. Diagnostic plots for model 2 predicting seagrass presence/absence by shoreline segment, a) Pearson residual
versus predicted value, b) Pearson residual by sediment class, and c) Pearson residual versus centered percent total organic
carbon.
28
-------
Table 6. Random effects for GLMM model 2 predicting shoreline presence/absence for
seagrass, not accounting for historic seagrass presence.
Shoreline (Intercept) exp(coeff)
-99
10
12
17
9
14
13
16
18
15
8
3
2
6
7
11
-4.30
6.55
4.91
3.34
2.78
2.47
1.58
1.37
0.67
0.21
-0.47
-2.11
-3.42
-3.47
-4.32
-5.79
0.01
698.81
136.15
28.33
16.14
11.78
4.85
3.94
1.95
1.24
0.63
0.12
0.03
0.03
0.01
0.00
Table 7. Fixed effects for GLMM model 2 predicting shoreline presence/absence for seagrass, not accounting for historic
seagrass presence.
Parameter
Coefficient Std.Error exp(Coeff) DF
t-value p-value
(Intercept)
csalAv
fZgtMXZav
csecchimin
cCG046avkm
cDistToHdS
cDstoMarin
cwindMin
cpttocMin
fSED87
fSED8124
fSED8512
fSED88101113
PR9wFA1370
Csecchimin x
cpttocMin
-1.16
1.57
1.55
0.31
0.033
-6.55E-04
-2.55E-04
-1.50E-05
-1.06
-0.56
-3.55
-0.79
-0.78
3.62
0.27
0.99
0.09
0.32
0.08
0.004
7.56E-05
5.10E-05
2.70E-06
0.06
0.21
0.48
0.22
0.23
0.06
0.10
0.3
4.8
4.7
1.4
1.0
1.0
1.0
1.0
0.3
0.6
0.03
0.5
0.5
37.5
1.3
19174
19174
19174
19174
19174
19174
19174
19174
19174
19174
19174
19174
19174
19174
19174
-1.17
17.49
4.79
3.68
8.80
-8.66
-4.99
-5.57
-17.47
-2.70
-7.47
-3.65
-3.48
57.62
2.78
0.2431
0
0
0.0002
0
0
0
0
0
0.0069
0
0.0003
0.0005
0
0.0054
29
-------
3.1.2.2 Shoreline P/A models with Serial Autocorrelation
We created an alternative model to predict shoreline presence/absence using historic
seagrass presence/absence as a predictor. The final model form was (Table 8}:
fsavcodeMaxi = fEGPA99Maxi + cAREAMax + cDistToEG99Mim + csalAvi +
csecchiminMaxi * cpttocMim + fZgtMXZavMax i + cCG046avkm2Mini +
cDstoMarinaMaxi + cwindMim * fSED6i+ PR6dwFA700,
where terms are defined as above (Table 3} with the addition of:
fSED6i = Majority sediment class at i, including
fSED6_6 = Gravelly-sand
fSED6_124571213 = Gravel, Sandy Gravel, Gravel-Sand-Silt, Sand, Sand-
Silt-Clay, and Silly-Sand
fSED6_81011 = Silty, Sandy Silt, Clay-Silt, and Gravel-Silt-Clay
PR6dwFA700i= Focal average over 700 meters distance of Pearson residual for
model fit
(3)
Table 8. Fixed effects model 3 components predicting shoreline presence and including historic 1999 seagrass predictors.
Parameter
Value exp(coeff) Std. Error DF
t-value p-value
(Intercept)
csalAv
csecchiminMax
cCG046avkm2Min
cAREAMax
cDistToEG99Min
cDstoMarinaMax
cwindMin
cpttocMin
fEGPA99Max
fZgtMXZavMax
fSED6124571213
fSED681011
PR6dwFA700
csecchiminMax x
cpttocMin
cwindMin xfSED681011
cwindMin x
fSED6124571213
-2.09
0.83
0.09
-0.04
0.000006
-0.00199
-0.00024
0.00015
-0.43
1.84
1.49
-1.09
-1.93
3.57
0.75
-0.0002
-0.00016
0.12
2.28
1.09
0.96
1.00
1.00
1.00
1.00
0.65
6.32
4.44
0.34
0.15
35.39
2.12
1.00
1.00
1.09
0.13
0.14
0.007092
2.4E-06
0.00011
0.000088
3.53E-05
0.11
0.37
0.50
0.50
0.54
0.10
0.18
4.83E-05
3.55E-05
19172
19172
19172
19172
19172
19172
19172
19172
19172
19172
19172
19172
19172
19172
19172
19172
19172
-1.9
6.2
0.6
-5.9
2.6
-18.1
-2.8
4.3
-3.9
5.1
3.0
-2.2
-3.6
35.4
4.1
-4.2
-4.6
0.0556
0
0.5475
0
0.0092
0
0.0057
0
0.0001
0
0.0028
0.029
0.0003
0
0
0
0
30
-------
Model predictions based on the original fit were slightly better than the prior model which
did not include historic eelgrass presence/absence, with a ROC value of 0.9753.
b)
CO
o
C-J
o
p
d
y*
s
1.0
0.3
0.6 0.4
Specificity
c)
o
CO
o
0.2
0.0
1.0 0.
0.6 0.4 0.2
Specificity
0.0
en •*.-
o
O
1.0
0.8
0.2
0.0
0.6 0.4
Specificity
Figure 10. Receiver operating characteristic (ROC) curve for final model predicting shoreline presence/absence of seagrass
based on a) original model 2 fit with full data set (ROC = 0.9546), b) same model with ten-fold model cross-validation (ROC =
0.9547), and c) original model 3 fit for model with full data set incorporating historic 1999 eelgrass presence/absence (ROC :
0.9753).
31
-------
3.1.2.3 Shoreline Segment Relative Frequency Models without Serial Correlation
The final model predicting shoreline relative frequency [average shoreline seagrass
occupancy along offshore transect] after incorporating residual autocorrelation was:
savcodeAv = csalAv + csalAv2 + csecchiminAv + csecchiminAv2 + cpttocAv + ,.-.
cpttocAv2 + cZgtMXZavAv + cZgtMXZavAv2 + cDistToHdShAv + cDstoMarinaAv
+ cUSRMARIkm2Av + fISOLATED + cwindAv * fSEDS + PR3xFA460,
where terms are defined as above with the addition of
PR3xFA460 = focal average of Pearson residuals with zone length of 460
meters
After residual autocorrelation was incorporated into the model, the magnitude of the
random effect term was very small compared to residual error, so the random shoreline
effects were dropped and a simpler general linear model was run. Similar to the grid-cell
occupancy model, several higher order effects were retained, but this time an additional
second-order term for depth greater than average wave mixing depth was included. In this
case the direction of effects for distance to nearest marina (positive influence] and
unsewered residential density (negative influence] was as expected, but the magnitude of
these effects was relatively small (Table 9]. Diagnostics plots showed improvements over
fits without the residual SAC term (Figures lla-d]. However, 3 outliers were apparent.
The model was fit after removing the three outliers and the same terms were retained.
Figure 12 illustrates the interaction between sediment class and wind-derived wave
energy, with wave energy effects greatest for silt-dominated classes, as compared to sand
or gravel-dominated classes. Performance of the model predicting average shoreline
occupancy by seagrass was slightly lower than predictions of presence/absence, but still
much more robust than model predictions at the grid cell scale. The area under the ROC
curve was 0.889 for model 4, and slightly less (0.8377] for the model 4 fit without 3
outliers (Figure 13].
32
-------
b)
B -
-15 -10 -505
Predicted values
Theoretical Quantiles
g -
•:/:: o.ee
Predicted values
0.10 0
Leverage
Figure 11. Diagnostic plots for model 4 predicting average shoreline occupancy, a) Pearson residuals vs predicted values, with
loess curve superimposed showing no trend, b) Q-Q plot (standardized deviance residuals versus theoretical quantiles), c)
Scale-Location plot (square root of standardized deviance residuals versus predicted values), and d) Residuals vs Leverage
with three labeled outliers.
33
-------
Table 9. Fixed effects in model 4 predicting average seagrass occurrence along transects perpendicular to shoreline, not incorporating historic seagrass presence.
Model with full data set Model without 3 outliers
Coefficients:
Estimate
Std. Error value
Pr(>|z|
Std.
Estimate Error
value Pr(>|z|
(Intercept)
csalAv
csalAv2
csecchiminAv
csecchiminAv2
cDstoMarinaAv
cwindAv
cDistToHdShAv
cUSRMARIkm2Av
cpttocAv
cpttocAv2
cZgtMXZavAv
cZgtMXZavAv2
fISOLATED
fSED581011
fSED5124561213
PR3xFA460
cwindAv x
fSED5124561213
-3.
1.
0.
-1.
1.
.52
.71
.65
.20
.07
3.89E-04
-6.53E-06
-4.83E-04
-0.013
-0.
-0.
-0.
-0.
1.
0.
-0.
.14
.58
.73
.70
.46
.30
.27
3.444E-01
-2.81E-05
0.03
0.07
0.05
0.05
0.03
1.56E-05
1.11E-06
3.11E-05
0.000
0.02
0.03
0.01
0.01
0.04
0.05
0.02
2.954E-03
1.97E-06
-103.6 <
24.5 <
12.0 <
-23.9 <
39.8 <
25.0 <
-5.9
-15.5 <
-35.5 <
-6.8
-23.0 <
-51.4 <
-51.1 <
37.0 <
6.3
-11.4 <
116.6 <
-14.3 <
2.
2.
2.
2.
2.
2.
3.
2.
2.
8.
2.
2.
2.
2.
2.
2.
2.
2.
OOE-16 ***
OOE-16 ***
.OOE-16 ***
.OOE-16 ***
OOE-16 ***
.OOE-16 ***
.51E-09 ***
OOE-16 ***
OOE-16 ***
.62E-12 ***
.OOE-16 ***
OOE-16 ***
OOE-16 ***
.OOE-16 ***
.35E-10 ***
OOE-16 ***
.OOE-16 ***
.OOE-16 ***
-3.30
1.45
0.72
-1.04
1.00
3.38E-04
-2.65E-06
-4.88E-04
-0.011
-0.43
-0.17
-0.57
-0.53
1.28
-0.007
-0.47
-3.63E-05
0.03
0.06
0.05
0.05
0.02
1.45E-05
1.22E-06
2.73E-05
0.000
0.02
0.02
0.01
0.01
0.04
0.046
0.02
1.87E-06
-107.3
22.6
14.9
-22.7
40.0
23.3
-2.2
-17.9
-30.5
-24.6
-9.0
-47.3
-47.6
34.7
-0.2
-21.4
-19.4
<2e-16
<2e-16
<2e-16
<2e-16
<2e-16
<2e-16
0.0303
<2e-16
<2e-16
<2e-16
<2e-16
<2e-16
<2e-16
<2e-16
0.8738
<2e-16
<2e-16
***
***
***
***
***
***
*
***
***
***
***
***
***
***
***
***
cwindAv xfSED581011
-7.40E-05 5.02E-06 -14.7 < 2.00E-16 ***
-9.11E-05 5.27E-06 -17.3 <2e-16 ***
34
-------
-------
3rd-
i.o
o.s
0.5 0.4
Specificity
0.2
C.C
Figure 13. ROC curve showing fit of model 4 prediction of average shoreline occupancy after removal of three outliers
3.1.2.4 Shoreline Segment Relative Frequency Models with Serial Correlation
With the effect of historic 1999 occupancy incorporated, the final model predicting average
shoreline occupancy was:
savcodeAv ~ cEGPA99Av + cAREAAv + cDistToEG99Av + csalAv + csalAv2 +
csecchiminAv + csecchiminAv2 + cpttocAv + cZgtMXZavAv + cZgtMXZavAv2 +
cDistToHdShAv + cDstoMarinaAv + cUSRMARIkm2Av + cwindAv * fSEDS +
PR4gFA460,
where terms are defined as above with the addition of
cEGPA99Av = centered average shoreline occupancy in 1999
cAREAAv = centered average 1999 eelgrass patch size by shoreline
index i
cDistToEG99Av = centered shoreline average distance to edge of
nearest 1999 eelgrass patch, and
PR4gFA460 = zonal average over 460 meters of Pearson residual for
model 4h
(5)
Relative seagrass presence increased in areas of historic 1999 seagrass, more so as historic
patch size increased and less so with distance from historic patch edge (Table 10}. As
above, diagnostic plots showed the presence of three outliers (Figure 14], but re-analysis of
model 4h without these outliers yielded the same set of predictors. Performance and
36
-------
Table 10. Fixed effects for model 5 predicting average seagrass presence/absence along transects perpendicular to shoreline, including effects of historic seagrass presence. Model coefficients
are compared for model fit with full data set and with data set minus three outliers.
Coefficients
(Intercept)
csalAv
csalAv2
cEGPA99Av
cpttocAv
cAREAAv
cwindAv
cDstoMarinaAv
cDistToHdShAv
cUSRMARIkm2Av
cDistToEG99Av
cZgtMXZavAv
cZgtMXZavAv2
csecchiminAv
csecchiminAv2
fSED581011
fSED5124561213
RR4gFA460
cwindAv xfSED581011
Full data set
Estimate Std. Error z value
Pr(>|z|
Data set minus three outliers
Estimate Std. Error z value
Pr(>|z|
-5.10
2.48
1.34
1.69
0.16
0.05
0.08
0.06
0.16
0.02
7.32E-06 6.68E-07
-3.18E-06 9.98E-07
-1.79E-04 1.76E-05
-3.50E-04 3.36E-05
-8.40E-03 2.89E-04
-1.03E-03 1.42E-05
-0.54 0.01
-0.56
-1.38
0.90
0.20
-0.03
0.457
0.01
0.05
0.03
0.05
0.03
0.004
-111 < 2.00E-16 ***
31.9 < 2.00E-16 ***
23.0 < 2.00E-16 ***
10.8 < 2.00E-16 ***
7.8 6.31E-15 ***
-6.70
3.47
2.04
0.98
0.16
0.06
0.08
0.06
0.16
0.02
11.0 < 2.00E-16 *** 1.02E-05 6.68E-07
-3.2 0.00144 ** -5.74E-06 1.01E-06
-112 < 2.00E-16 ***
41.7 < 2.00E-16 ***
32.4 < 2.00E-16 ***
6.2 5.03E-10 ***
7.3 2.75E-13 ***
15.2 < 2.00E-16 ***
-5.7 1.34E-08 ***
-10.2 < 2.00E-16 *** -2.71E-04 1.82E-05 -14.9 < 2.00E-16 ***
-10.4 < 2.00E-16 *** -2.67E-04 3.62E-05 -7.4 1.55E-13 ***
-29.0 < 2.00E-16 *** -8.00E-03 2.91E-04 -27.5 < 2.00E-16 ***
-72.8 < 2.00E-16 *** -1.52E-03 1.86E-05 -81.7 < 2.00E-16 ***
-0.47 0.01
0.01
0.05
0.03
0.05
0.03
-37.1 < 2.00E-16 ***
-40.9 < 2.00E-16 ***
-28.6 < 2.00E-16 ***
34.3 < 2.00E-16 ***
4.1 4.93E-05 ***
-1.1 0.27671
0.004 114.2 < 2.00E-16 ***
-0.34
-1.60
0.95
0.29
0.23
0.61
-32.9 < 2.00E-16 ***
-27.7 < 2.00E-16 ***
-31.2 < 2.00E-16 ***
34.8 < 2.00E-16 ***
5.8 6.23E-09 ***
7.8 4.56E-15 ***
0.01 117.9 < 2.00E-16 ***
-1.29E-06 5.58E-06
cwindAv xfSED5124561213 -1.50E-05 2.20E-06
-0.2 0.81674 1.66E-05 5.28E-06 3.1 0.00167 **
-6.8 8.90E-12 *** -2.64E-05 2.53E-06 -10.4 < 2.00E-16 ***
37
-------
robustness of the model were actually degraded by including historic seagrass presence as
a predictor. The area under the ROC curve was 0.615 for model 4h but only 0.4887 for the
same model with three outliers removed (Figure 15}.
b)
-n O -.,
-15 -10 -5
Predicted values
c)
-202
Theoretical quantiles
-15 -10 -505
Predicted values
Figure 14. Diagnostic plots for model 5 predicting average shoreline occupancy including historic 1999 eelgrass predictors, a)
Residuals vs predicted with loess plot overlaid, b) Normal Q-Q plot, and c) Scale-Location plot. Note the presence of three
outliers.
38
-------
b)
1.0
0.8
0.6 0.4
Specificity
0.2
0.0
•DO
O
(N
O
O
O
1.0 0.8 0.6 0.4 0.2 0.0
Specificity
Figure 15. ROC curve for model 5 a) before and b) after removal of three outliers.
3.1.2.5 Shoreline Segment Minimum Depth of Occurrence
Based on a comparison of AIC values, the final GLM model for predicting minimum depth of
seagrass occurrence is:
(6)
bathymmin = fSED4 * halflen0708gsmMax,
where bathymmin = minimum depth of seagrass occurrence (m}
fSED4 = sediment class
halflen0708gsmMax = maximum wave mixing depth over growing
season
Maximum wave mixing depth proved to be a better predictor of minimum seagrass depth
of occurrence than the average value or average + 2SD (Table 11}. Seagrass occurring on
sediments of class 12 (sand-silt-clay} were most sensitive to wave mixing depth (Figure
16}. Residuals showed no evidence of nonlinearities, a slight tendency for increasing
variance with the mean, and were close to a normal distribution but with three outliers
identified (Figure 17}.
39
-------
Table 11. Fixed effects for model 6 predicting seagrass minimum depth of occurrence by shoreline distance.
Coefficients: Estimate Std. Error t value Pr(>|t|)
(Intercept)
halflen0708gsmMax
fSED412
fSED47
fSED4810
fSED4810 x
halflen0708gsmMax
fSED47 x half Ien0708gsm Max
0.30
0.03
5.53
0.74
0.25
0.004
-0.04
0.09
0.01
1.30
0.11
0.25
0.019
0.01
3.4
4.7
4.3
6.5
1.0
0.2
-4.3
8.07E-04
2.89E-06
2.34E-05
1.50E-10
0.32
0.82
1.65E-05
***
***
***
***
***
fSED412 x
halflen0708gsmMax
-0.37
0.10
-3.8 0.000158 ***
0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
i i i i i i i i i i i i i i i i
3-
2-
1-
0-
-1 -
-2-
-3-
i • ii
fSED:4:12
fSED4:5
f • T
fSED4:810
in • ii
fSED4 :7
ni • i • • . • . • .
4
3
2
1
0
-1
-2
-3
0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
Max of 2007-08 growing season average wave mixing depth
perpendicular to shoreline [m)
Figure 16. Interaction of sediment class and maximum wave mixing depth on minimum depth of seagrass occurrence (model
6). SED4 class definitions are: 12 = sand-slit-clay; 810 = silty and sandy-silt; 5 = sand; 7 = silty sand.
40
-------
aj
b)
0.4 0.6 0.8 1.0 1.2 1.4 1.6
Predicted values
c)
-3-2-1012
Theoretical Quantiles
0.00
0.01
0.02
Leverage
0.03
0.04
Figure 17. Diagnostic plots for model predicting minimum seagrass depth of occurrence (model 6) a) Residuals versus
predicted values, b) Q-Q plot of residuals(dashed line shows expectation for normal distribution), and c) residuals versus
leverage, highlighting presence of three outliers (labelled).
41
-------
3.1.2.6 Shoreline Segment Maximum Depth of Occurrence
Based on a comparison of AIC values, a GAM model provided a superior fit over GLM
models in predicting maximum depth of seagrass occurrence (Figure 18}:
bathymmax = s(cLlOsecchiavMax) + s(cLlOpttocMax) +
s(cL10USRMARIkm2Max)
where s = smoothing function
r2adj= 0.2
cLIOsecchiavMax = centered logic maximum of seasonal average
Secchi depth(m]
cLIOpttocMax = centered logic maximum of sediment total organic
carbon (%}, and
cL10USRMARIkm2Max = centered logic maximum of unsewered
residential units on high infiltration soils/km2
(7)
Although the GAM model provided the best fit from a statistical viewpoint (lower AIC value
of 6838.7 as compared to 7110.6], the best GLM model may provide a more practical and
realistic option, i.e., without overfitting the data. Therefore, we are presenting the best
GLM model fit as well, determined using the step option in GLM, which had an AIC of
7110.6 as compared to the next best model with an AIC of 7112.2:
bathymmax = csecchiminMax + cpttocMax + cUSRMARIkm2Max +
csecchiminMax:cpttocMax + csecchiminMax:cUSRMARIkm2Max +
cpttocMax:cUSRMARIkm2Max,
where bathymmax = maximum depth of seagrass occurrence at shoreline
index i
csecchiminMax = centered maximum of Secchi depth seasonal
minimum at shoreline index I, and
cpttocMax = centered maximum of sediment percent total organic
carbon at shoreline index i
(8)
Although main effects for minimum Secchi depth and density of unsewered residences on
high infiltration soils were not significant, they were retained in the model because these
terms were included in significant interaction terms. Maximum depth of seagrass
42
-------
a)
r- 10
p ,_
CO
x"
m
I «
ro o
o
£ -n
o o
T— —I. I
V I
-0.3
b)
-0.2 -0.1 0.0
Centered Iog10 maximum average Secchi depth (m)
0.1
00
•*r
oo
CM
E
E
Bj o
CD
Hi IT,
<••» •* '
I I I
I I I
TF
I M i
C)
(O
00
00
x"
03
-0.5 0.0 0.5 1.0 1.5 2.0 2.5
Centered log,Q maximum unsewered residences on high infiltration soils/km2
uo
o
o
o
ti
Q.
O IT)
^ 9
J I I I
; i i
-1.0 -0.5 0.0 0.5
Centered Iog10 maximum sediment total organic carbon (%)
Figure 18. Smoothing functions in generalized additive model 7 to predict maximum depth of seagrass occurrence by
shoreline position, a) Centered logio shoreline maximum average Secchi depth (m), b) Centered logio shoreline maximum
unsewered residences on high infiltration soils/km2, and c) Centered logio shoreline maximum sediment total organic
carbon (%).
43
-------
occurrence tended to increase with sediment total organic carbon at high Secchi depths but
decrease with TOG at low Secchi depths (Figure 19}. Interaction effects involving density of
unsewered residences were significant but weak (Table 12}.
Table 12. Fixed effects in generalized linear model 8 predicting maximum depth of occurrence for seagrass at shoreline
index i.
Std.
Coefficients: Estimate Error t-value Pr(>|t|)
(Intercept)
cUSRMARIkm2Max
csecchiminMax
cpttocMax
csecchiminMax x cpttocMax
csecchiminMax x
cUSRMARIkm2Max
cpttocMax x cUSRMARIkm2Max
1.43
0.0003
-0.02
-0.26
0.63
0.004
0.003
0.02
0.0008
0.03
0.03
0.05
0.001
0.001
79.3 <
0.4
-0.6
-8.7 <
13.7 <
3.1
1.9
2.00E-16
0.72197
0.53662
2.00E-16
2.00E-16
0.00166
0.05699
***
***
***
**
44
-------
-1.0 -0.5 0.0 0.5 1.0 1.5 2.0
_E
x
3.5
3.0
2.5
2,C
1.5
1.0 '
0.5 '
i illinium linn i ii
i ii miiiiiiiiiii i n mil in
V
:;;:: •••'.•=•
iiiiiiiiiiiiiiiiiiiiiiiiii
-1.0 -05 0.0 0.5 1.0 1.5 2.0 _ __ -1.0 -0.5 0.0 0.5 1.0 1.5
Centered maximum sediment total organic carbon (%)
ii
2.C
3.5
3.0
2.5
2.0
1.5
1.0
0.5
Figure 19. Interactive effects of centered sediment total organic carbon and shoreline maximum of seasonal minimum Secchi
depth on maximum seagrass depth at shoreline index I (model 8). Sediment TOC has a stimulatory effect at high Secchi
depths but an inhibitory effect at low Secchi depths. Secchi depth for each plot is indicated by brown vertical line in each tan
horizontal header.
45
-------
Chapter 4. Discussion
4.1 Data and modeling challenges and constraints for statistical predictive
seagrass models
Development of predictive statistical models for seagrass occurrence is inherently
challenging, but the problems can be addressed through application of a framework to
describe spatial autocorrelation of seagrass patches, new techniques for incorporating SAC
terms into predictive models, and careful selection of independent variables to avoid cross-
correlations. Problems to be solved include the need to predict at fine spatial resolutions
over large areas with associated memory constraints, spatial autocorrelation of model
residuals with anisotropy, serial correlations in time series, and potential cross-correlation
of predictive variables that co-vary with water depth or with salinity gradients from head-
of-tide to estuarine mouth. We successfully applied strategies to limit difficulties with
parameter cross-correlations: using centered variables, carefully selecting the reference
sediment class against which to compare other sediment factors, choosing alternative
strategies to relate transparency and mixing depth to water depth as a difference versus
binary variable, and dropping correlated variables (e.g., temperature] of lesser interest.
After quantifying the actual range of SAC both parallel and perpendicular to the shoreline,
we were successful in incorporating terms for spatial autocorrelation into models using a
modified version of the Crase etal. (2012} approach. Although memory constraints limited
the number of R packages available for analysis, we were successful in applying the
glmmPQL package.
Our modelling approach and national data sources should be applicable to other estuaries,
with substitution of localized data sources for finer resolution data or more site-specific
information on seagrass disturbances where appropriate. With some exceptions, the
significant components of the predictive models developed here for seagrass presence are
available for estuaries across the United States: bathymetry, salinity, grain size, wave
mixing depth (derived from bathymetry and wind speed and direction], sediment percent
total organic carbon, density of unsewered residences on high infiltration soils, and historic
seagrass extent. For some variables, we used more localized sources of data because of
their improved spatial resolution, e.g., sediment grain size. While satellite imagery is
available to predict optical properties of seawater, most algorithms developed to date are
not appropriate for coastal waters (< 30 meters depth] and so we used local monitoring
data for Secchi depth. As coastal algorithms for optical parameters are improved,
appropriate satellite data will become more readily available (Keith et al. 2014]. We were
able to use NOAA's WEMo model to predict wave energy (Fonseca and Malhotra 2010];
however, this model has not been updated for versions of ArcMap beyond 9.3 so will not be
readily available to other users. It is possible that coarser estimates of relative wave
energy such as metrics within the USGS Coastal Vulnerability Index (Hammar-Klose and
Thieler, 2001] could serve as a proxy for WEMo model values; that option has not been
tested here. Other researchers will have access to an updated version of the USGS WAVES
extension for ArcMap, which provides estimates of wave mixing depth. We relied on some
local maps of disturbance indicators (distance to hardened shoreline and marinas, Canada
goose density] but these effects were marginal in models.
46
-------
It is possible to incorporate spatial autocorrelation in predictive models for seagrass,
rather than avoiding it by sampling at coarser scales to avoid SAC, and ignoring the
inherent patchiness of seagrass distributions. Seagrass currently occurs in narrow linear
patches within Narragansett Bay and thus we modeled seagrass occurrence at a fine spatial
scale (10 meter grid cells}. Even after restricting model boundaries to the zero to five
meter depth zone and to the parameter space within which seagrass is currently found in
Narragansett Bay, we were still faced with the challenge of developing models with a very
large data set (518856 10 x 10 meter grid cells or 19172 10-meter length shoreline
locations}. The magnitude of our data set restricted the statistical packages and
approaches we could use for binomial GLMMs and GAMMs due to memory and processing
time constraints. However, incorporation of a residual autocorrelation term in models
allowed us to account for spatial autocorrelation in models in a less memory-intensive
manner than would have been required by simultaneous fitting of spatial autocorrelation
model distributions (with anisotropy}. We were able to create realistic representations of
spatial patterns (patchiness} in seagrass distribution by 1} incorporating shorelines as a
random "founder" effect, 2} defining spatial coordinates relative to longshore and offshore
distance, and 3} calculating moving averages of residual autocorrelation based on zones
reflecting different ranges of spatial autocorrelation in longshore versus offshore
directions.
Unexplained variation in our models could have been related to factors we did not include
in predictive models due to lack of available data, including: biotic effects
(epiphyte/grazing interactions, macroalgae competition, bioturbation and bioirrigation
effects by fauna [Nelson 2009]}, sediment sulfide/redox (Goodman et al. 1995}, tidal range
effects on minimum depth potentially interacting with light availability (Koch et al. 2001},
wave energy associated with winter storms (Kelly et al. 2001}, historic hurricanes or
tropical storms, tidal currents, anthropogenic disturbance (shellfish dredges or rakes,
anchor beds, propeller scars; see Neckles et al. 2005, Oakley et al. 2013}, and restoration
activity. We also failed to include all potential interaction terms in the models due to
problems with variance inflation factors. Any of these factors could have contributed to
unexplained variation; yet, despite this and the other challenges inherent in predicting
seagrass distributions, our prediction of seagrass absolute or average presence/absence at
shoreline locations was very robust, with area-under-the-curve (AUC} values associated
with Receiver Operating Characteristic (ROC} curves of 0.95 - 0.98
4.2 Comparative performance of alternative modeling approaches
This was the first study to explicitly address spatial autocorrelation in eelgrass models
using mixed models. Failure to account for spatial autocorrelation in species distribution
models can inflate the significance of explanatory variables (Crase etal. 2012}. Most
researchers developing predictive models for eelgrass habitat have either ignored the
potential for spatial autocorrelation and effects on variable selection (e.g., Krause-Jensen et
al. 2003, Grech and Coles 2010} or have selected sample points at a minimum distance
apart to avoid spatial autocorrelation (Downie etal. 2013}. Although the latter strategy
takes care of the issue from a statistical standpoint, it fails to acknowledge the positive
feedback effects of seagrass on sediment stability and the light environment, and cannot
successfully predict the patchy nature of seagrass distribution. However, patch dynamics or
47
-------
cellular automata models (Wortman et al. 1997, Fonseca et al. 2000b, Fonseca et al. 2004,
Schonert and Milbradt 2005} do predict this patchiness. Incorporation of a term for spatial
autocorrelation allowed us to avoid overestimation of percent variance explained and to
reduce the number of significant predictors in our final regression models to an
appropriate number. Examples of terms dropped following addition of SAC error terms
include the shoreline isolation term in some models and the number of sediment classes
among which we could detect differences in either main effects or interactions with wind
energy or wave mixing depth.
Models predicting seagrass presence were more robust at the shoreline scale than at the
individual grid cell (Table 13}. Between original tests and 10-fold cross validation results,
model performance measurements declined considerably based on 10-fold cross-validation
at the grid-cell scale but stayed constant at the shoreline scale. Model performance was
slightly better for prediction of presence/absence as compared to average
presence/absence at the shoreline segment. Model performance for prediction of average
shoreline presence/absence was significantly degraded after information on historic
eelgrass presence/absence was included in models.
Table 13. Summary of model performance at different scales, based on resubstitution or 10-fold cross-validation and with or
without incorporation of data on historic seagrass coverage.
Dependent
variable
Grid cell P/A
Grid cell P/A
Shoreline max P/A
Shoreline max P/A
Shoreline max P/A
Shoreline max P/A
Shoreline avg P/A
Shoreline avg P/A
Shoreline avg P/A
Shoreline avg P/A
Historic
P/A
included?
No
No
No
No
Yes
Yes
No
No
Yes
Yes
Resubstitution
or 10-fold cv?
R
lOxcv
R
lOxcv
R
lOxcv
R
lOxcv
R
lOxcv
AUC
0.98
0.71
0.95
0.95
0.98
0.98
0.89
0.72
0.62
0.62
We have chosen the AUC statistic to compare performance across our models and between
our models and those of others because it is insensitive to the probability of occurrence of
the species of interest, unlike overall prediction accuracy and some other measures
(Fielding and Bell 1997}. Some researchers present only an overall prediction accuracy, so
we are unable to compare their results with ours. The AUC statistic represents the area
under a ROC curve, and varies between 0 and 1, with a value of 0 representing zero
predictive power, 0.5 representing predictions no better than chance, and 1 representing
perfect predictions. The ROC curve is a plot that demonstrates the performance of a model
predicting a binary variable. It is created by plotting the true positive rate (sensitivity}
against the false positive rate (1 - specificity} as the discrimination threshold is varied.
48
-------
Our models predicting seagrass P/A at the 10 x 10 meter grid cell scale did not perform as
well as some existing models predicting P/A at 25 x 25 meter or 50 x 50 meter grid cell
scales, but our model performance for predicting shoreline occurrence exceeded earlier
reported AUC statistics. Valle et al. (2013} calculated AUC values with 5-fold cross-
validation to compare regression approaches (GLM, ANN, MARS, GAM] with machine
learning methods (BRT, RF] to predict seagrass presence/absence at the (50m x 50m] grid
cell scale. At a grid-cell scale 25 times coarser than our model predictions, their best GLM
model performance had a mean AUC of 0.84, while their best models based on boosted
regression tree methods yielded a mean AUC of 0.94. Using an independent test set
comprised of 30% of available data, Downie etal. (2013} compared the model performance
using GAM and maximum entropy modeling (Maxent] and found model AUC values of 0.84
(GAM] and 0.80 (Maxent} at a grid cell scale of 25 x 25 meters. It is possible that our model
performance at the grid cell scale is lower than that obtained by other researchers because
of the finer scale of our grid cells, or because of differences in study area characteristics and
the predominant form of Z. marina (intertidal annual versus sub tidal perennial] found in
each setting. Our models predicting shoreline presence were robust and yielded predictive
accuracy (AUC values] equal to or better than the results of Downing et al. (2013] or Valle
et al. (2013]. It is important to note, however, that our models incorporated a random
shoreline effect which explained much of the background variability within the estuary
(see discussion below].
Although GAM models can perform better than GLM models based on AIC and
AUC values, they may overfit the data and can be replaced with GLM models incorporating
higher-order terms to capture nonlinearities. In the one case where we could derive a GAM
model for comparison with a GLM to predict maximum depth of eelgrass occurrence, the
statistical fit of the GAM model was better based on AIC values. Likewise, Valle et al. (2013]
found better performance for GAM models (mean AUC = 0.93] as compared to GLM models
(mean AUC = 0.84] in predicting eelgrass presence/absence. However, the multiple peaks
we observed in response curves had no mechanistic basis, suggesting that GAM models
were probably overfitting the data. It is likely that higher order terms in GLM models are
sufficient to capture nonlinear responses in logistic models. Addition of a second order
term can yield an S-shaped probability plot, while addition of a third-order term can
describe a unimodal response to an environmental gradient.
4.3 Relative importance of different factors in predicting seagrass P/A
4.3.1 Seagrass sensitivity to different environmental variables in Narragansett Bay
We can compare the sensitivity of seagrass survival and persistence to different
environmental variables by comparing the potential change in ln(odds ratio] over the
range of each predictor in the Narragansett Bay area of interest. For predictors involving
higher order terms, we can predict minimum and maximum potential contributions. For
interaction terms, we can examine sensitivity at high and low ends of modifying factors
(Table 14]. For the model predicting seagrass presence/absence at the grid cell scale, the
most influential predictor was Secchi depth, followed by, in order: shoreline isolation,
sediment percent total organic carbon, sediment type, and salinity. The least influential
49
-------
variable was water depth greater than average wave mixing depth. For the model
predicting presence of seagrass at shoreline locations, the most influential predictor is
sediment type, followed by percent total organic carbon (at low Secchi depth], then salinity.
4.3.2 Seagrass sensitivity to different environmental variables in other regions
Managers need to consider the energy environment of different regions and estuarine
hydrogeomorphic types and different growth forms in comparing the predictions of
seagrass models and the relative influence of different predictor variables. Most of the
regression or machine-learning models developed for predicting presence/absence of
Zostera marina are based on analyses of western Europe datasets, e.g., the Baltic Sea and
Wadden Sea (Krause-Jensen et al. 2003, van der Heide et al. 2009, Downie et al. 2013, Valle
et al. 2013,} or of the Great Barrier Reef (Grech and Coles 2010}. In many of these cases,
the energy regime (relative wave exposure, current velocity} and/or substrate are the
predominant variables predicting eelgrass presence (Grech and Coles 2010, Downie etal.
2013, Valle et al. 2013}. However, investigators predicting eelgrass presence over larger
regions have found that light availability (Krause-Jensen et al. 2003, van der Heide et al.
2009}, total N in surface water (van der Heide et al. 2009}, and, to a lesser extent, salinity
(Krause-Jensen et al. 2003} are the driving factors predicting eelgrass presence.
Results from the models developed for western European regions may not be comparable
to our results for Narragansett Bay because they focus on higher energy environments (or
those with more extreme physico-chemical gradients}. In addition, Zostera marina
populations in the Wadden Sea are comprised of annual forms which occur in intertidal
zones dominated by frequent disturbance and less subject to light limitation. Both the
robust perennial form of Z. marina and the more flexible annual form occurred in the
Wadden Sea prior to the incidence of wasting disease in the 1930s, but the perennial form
never recovered (van Katwijk et al. 2009}. Narragansett Bay populations of Z. marina are
currently subtidal and are presumably composed of the perennial form. Most studies of the
factors affecting eelgrass presence/absence have ignored the distinction between annual
and perennial forms of Z marina. The annual form of Z. marina has been reported in Nova
Scotia and Maine estuaries (Keddy and Patriquin 1978} and in Ninigret Pond, RI (Thorne-
Miller et al. 1983} but may be present in other New England estuaries where eelgrass
occurs in intertidal zones (e.g., Great Bay}.
50
-------
Table 14. Potential effect of independent variables across range of predictors. For sediment types, the range of coefficients is
given. Effects are expressed in terms of ln(odds ratio). Maximum effect for Secchi depth is at an intermediate value of the
predictor, not the maximum.
Model
7
7
7
7
7
7
7
7
7
8
8
8
8
8
8
8
8
8
8
8
8
8
8
Predictor
cSDavgrtrZ, cSDavgrtrZ2,
cSDavgrtrZ3
f ISOLATE Dl
cptTOC, cptTOC2, cptTOC3
fSED4
cSAL, cSAL2
cUSRMARIkm
cDstoMarin
cCG046avkm
fZgtMXZavTl
fSEDS
cpttocMin at min
csecchimin
csalAv
cpttocMin at max
csecchimin
cDistToHdS
fZgtMXZavMTRUE
cwindMin
cCG046avkm
cDstoMarin
csecchimin at max
cpttocMin
csecchimin at min
cpttocMin
csecchimin
cpttocMin
csecchimin:cpttocMin
Min
-3.00
0
-1.26
0
-1.62
-31.02
-2242.6
-3.60
0
0
-1.72
-537
0
-8572
-3.9
-328
-0.84
-1.13
Max
4.23
1
2.24
1
0.48
285.18
252.71
28.60
1
1
0.38
1965
1
73428
28.3
2174
2.16
2.37
Coeff
1.65, -0.149, -0.053
-5.42
-0.91, 0.65, -0.33
-1.28 to 1.89
1.56,0.18
0.004
-0.0005
0.013
0.35
-3. 55 to 3.62
1.57
-6.55E-04
1.55
-1.50E-05
0.033
-2.55E-04
0.31
-1.06
0.27
Min
effect
-6.14
0.00
-3.76
-1.28
-3.33
-0.12
1.12
-0.05
0.00
-3.55
1.46
-2.70
0.54
0.35
0.00
0.13
-0.13
0.08
-0.801
-0.004
Max
effect
1.09
-5.42
1.56
1.89
-0.49
1.14
-0.13
0.37
0.35
3.62
-3.05
0.60
-1.13
-1.29
1.55
-1.10
0.93
-0.55
-0.055
0.010
Effect
range
7.23
-5.42
5.33
3.17
2.84
1.26
-1.25
0.42
0.35
7.17
-4.51
3.30
-1.67
-1.64
1.55
-1.23
1.06
-0.64
0.75
0.01
Rank
1
2
3
4
5
6
7
8
9
1
2
3
4
5
6
7
8
9
9
10
4.3.3 Limitations to existing regression models to predict seagrass in U.S. estuaries
In contrast to our approach, most regression models predicting eelgrass presence and/or
cover in U.S. estuaries include a smaller subset of variables at a time, generally focusing
either on nutrients and/or light availability (Duarte 1991, Latimer and Rego 2010, Benson
et al. 2013, Kenworthy et al. 2014} or on the energy regime (Kelly et al. 2001} but not
considering additive or interaction effects (Koch 2001}. Developers of habitat or
transplant suitability indices for eelgrass have considered a combination of eelgrass
colonization sources, substrate characteristics, wave exposure, and water quality/light
51
-------
environment based on data from empirical or experimental studies but tend to give these
factors equal weight and do not consider potential factor interactions. In more focused
studies, researchers have provided empirical and/or experimental evidence for interactive
effects of light availability with sediment organic matter content (Wicks et al. 2009,
Kenworthy et al. 2014} or of salinity with nutrients (van Katwijk et al. 1999}.
4.3.4 Comparison of estimates for light compensation depth and optimum light levels
Interpretation of predictions of light compensation points (maximum depth occurrence}
and optimum light levels depends on the vertical datum of the merged topobathymetric
grid used to estimate water depth. We used NOAA's Coastal Relief Model (CRM}, in which
source bathymetric data retained their original vertical datum of either mean lower low
water (MLLW} or mean low water (MLW}, while source topographic data remained in
either North American Vertical Datum (NAVD} 88 or National Geodetic Vertical Datum
(NGVD} 29 prior to merger. Given the semi-diurnal nature of tides in Narragansett Bay, the
MLLW value is the tidal datum of interest. Mean tidal range between MLLW and Mean
Higher High Water (MHHW} varies from 0.876 to 1.55 meters at NOAA tidal stations in
Narragansett Bay having a recorded datum
(http://tidesandcurrents.noaa.gov/stations.html?type=Datums}, with an average value of
1.27 meters.
Our model for seagrass presence/absence at the grid cell scale yields comparable
information on effects of light availability on seagrass maximum depths to values in the
literature. Our model for seagrass presence at the grid scale predicts a maximum
probability of occurrence at a centered Secchi depth (average Secchi depth - water depth}
value of 2.4 meters, which corresponds to a precentered value of 2.4 + (-0.433} = 1.97
meters. Our model also predicts an odds ratio of 1 (p = 0.5} at a centered difference value
of 0.985 meters, corresponding to a pre-centered value of 0.985 - 0.433 = 0.552 meters. If
we make the assumption that eelgrass light limitation depends on the mean tidal level, then
the model-predicted compensation depth would be 0.552 - (1.27/2} = -0.083 meters
greater (0.083 meters less} than the Secchi depth, or (2.44-0.083}/2.44 = 0.96 times
average Secchi depth. Duarte's compilation of compensation points for Z. marina from the
literature yields a prediction of Zc = 1.86/K, as compared to Dennison's (1987} value of
1.62/K for northeastern U.S. estuaries, and Nielsen and colleague's (1989} value of 1.53/K.
Using Poole and Atkin's (1929} relationship between the light attenuation coefficient, K,
and Secchi depth (K = 1.7/S.D.}, these correspond to a range of 0.9 - 1.11S.D. for Zc. Our
estimate of 0.96 times average Secchi depth at an odds ratio of 1 (probability of 0.5} falls
well within the range of values reported by Duarte (1991} for the light compensation
depth.
Although our model predicted an optimum light level for seagrass growth, the shoreward
limits to seagrass are more likely related to physical disturbance. Most researchers have
only evaluated maximum depth limits for eelgrass related to light limitation, with no
estimation of optimum light levels. Krause-Jensen et al. (2003} found that the probability
of eelgrass occurrence in Danish coastal waters increased up to 60% of surface irradiance
(at about 4 meters depth}, and then cover declined at higher values. Although Krause-
Jensen suggested that although this could have been related to photoinhibition at higher
levels, our models suggest it was more likely related to the increased probability of
52
-------
physical disturbance at shallower depths. Optimum light values for eelgrass in our model 7
occurred at depths of 1.97 meters less than average Secchi depth, or approximately 2.44 -
1.97 + (1.27/2} = 1.105 meters depth relative to Mean Tidal Level (MIL] for average
Secchi. A Secchi depth of 2.44 meters corresponds to a Kd of 0.70, which would yield a light
intensity of 1% incident light at 6.6 meters depth.
Our model for seagrass shoreline occupancy incorporated an interaction term between
water depth relative to Secchi depth and sediment percent total organic carbon which is
consistent with the literature. This interaction has also been observed by Kenworthy et al.
(2014} for southeastern MA embayments, with a decrease in eelgrass compensation depth
as sediment organic matter increases. This is likely due to the increase in sulfide content of
sediments with increased organic matter content and associated toxicity (Goodman et al.
1995, Holmer and Bondgaard 2001}.
4.3.5 Seagrass sensitivity to factors other than transparency
Our models showed mixed evidence for the effect of unsewered residences on high
infiltration soils, an indicator of potential groundwater N inputs. In two cases, model
coefficients were negative, in one case positive, and in one case the variable showed weak
interaction effects with other nutrient-related variables. Valiella's model incorporates the
potential effect of unsewered residences on nitrogen loads (Valiella et al. 2004}, but
applications of his model to predict load effects on seagrass have not attempted to separate
the effects of reduced transparency related to phytoplankton biomass with other effects
related to groundwater inputs (Latimer and Rego 2010}. It is possible that groundwater
DIN inputs could favor the growth of macroalgae and/or periphyton at the expense of Z.
marina (Harlin and Thornemiller 1981, Costa 1988, Teichberg etal. 2010}. Unsewered
residence density might have shown a stronger and more consistent effect in our models if
we had accounted for variability in residence time in various subembayments within
Narragansett Bay as was done by Latimer and Rego (2010}.
The range of salinity values encompassed by effects predicted by our model (26.6 - 28.7} is
well within the tolerance ranges reported for Zostera marina populations in the literature:
5 to 35 Practical Salinity Units (PSU} in the northern hemisphere (den Hartog 1970} but
slightly greater than ranges of 14 to 22 PSU reported for the Chesapeake Bay (Wetzel and
Penhale 1983}. Salinity effects predicted by our models could have been related to either
direct negative effects of salinity and energetic costs of osmoregulation, the interactive
effects of salinity and DIN, or the correlation of the salinity gradient with the gradient of
dissolved inorganic nitrogen (DIN} and total nitrogen (TN} in Narragansett Bay (Krumholz
2012}. Salinity and nitrate can have synergistic interactions, possibly related to the
tendency to incorporate low C (high N} amino acids as tissue nitrogen levels increase in
contrast to the need to generate the high C/low N amino acid proline involved in
osmoregulation (van Katwijk et al. 1999}. Van Katwijk et al. (1999} observed a decreased
tolerance of Z. marina to eutrophication at salinity levels of 26-30 ppt. It is likely that the
negative coefficient for salinity in our models reflects the negative impacts of water column
TN, and possibly an interaction between salinity and TN effects.
We found a strong interactive effect between wave mixing depth and sediment type on
minimum depth of seagrass occurrence, with the effect strongest on gravel and sandy
53
-------
gravel sediments, which probably also represent high energy environments. However,
wave mixing depth was the least influential predictor in our eelgrass model predicting
seagrass presence at the grid cell scale. Numerous researchers have reported a negative
effect of relative wave energy on seagrass presence, both in western Europe systems
(Downie et al. 2013, Valle et al. 2013} and in the United States (Kelly et al. 2001], but these
cases represent systems with more extreme energy gradients. Fewer researchers have
tried to test the relationship of wave mixing depth with minimum depth extent of eelgrass.
Chiscano (2000 in Koch 2001} found a poor relationship between minimum depth of
occurrence and wave mixing depth in shallow portions of the Chesapeake Bay with gentle
slopes. Our models indicated that wave mixing depth in combination with sediment
particle size was a significant factor determining minimum depth of seagrass occurrence.
The effects of sediment characteristics on habitat suitability for eelgrass have historically
been examined in isolation. Most evaluations of particle size effects have been correlative
in nature (Nelson 2009], but two studies did demonstrate greater growth rates on finer
sediments than on coarse sand or sand plus gravel, possibly related to nutrient availability
(Short 1987, Thorn et al. 2001}. Earlier evaluations of the effect of sediment organic matter
on eelgrass presence tended to be correlative as well, but more recent studies have
elucidated the same interaction between light compensation point (maximum depth of
occurrence} and percent organic matter in sediment types that we observed (Kenworthy et
al. 2014}.
Somewhat surprisingly, our shoreline model for seagrass relative abundance was not as
robust after we incorporated historic seagrass presence, distance, and patch size as
predictors in the model, suggesting that seagrass patch distribution may represent a
dynamic equilibrium with significant transition probabilities. In examining a 14-year time
series of intertidal eelgrass distributions in the Ems estuary in the Netherlands, Valle et al.
(2013} found transition probabilities of 12.7% for colonized areas and 12.9% for areas that
had disappeared, suggesting that patches are relatively mobile in this system even though
total area is relatively stable. Valle et al. (2013} characterize these intertidal populations as
r-selected in response to the high disturbance environment they are found in. Likewise,
Kelly et al. (2001} predicted that 16% of the seagrass area in their system was highly
susceptible to acute storm events. These studies suggest that predictive models will have
an inherent error rate and possibly a tendency to overpredict seagrass presence, i.e., a
significant proportion of suitable habitat may be unoccupied during any given year.
Moreover, recolonization rates maybe slower in the Northeast, leading to less patch
predictability from year to year. Based on observed trajectories, Neckles et al. (2005}
projected recovery times of over 10 years for dragged eelgrass patches in Casco Bay, Maine.
Recovery rates may also depend on the availability of seed sources and the presence of
annual versus perennial forms of Z. marina (Jarvis and Moore 2010} which will determine
the rate of vegetative spread versus sexual reproduction. The poor contribution of location
of historic seagrass patches to prediction of current seagrass patch locations suggests that
Narragansett Bay is not in a state of equilibrium and that patch locations may vary from
year to year, perhaps due to lags in recovery from disturbance.
54
-------
Random effects of shoreline code showed extreme differences in probability of seagrass
colonization, possibly due to hysteresis effects. This was surprising, given the widespread
historic occurrence of seagrass in Narragansett Bay prior to the incidence of wasting
disease in the 1930's (Figure 20}. It is possible that this discrepancy in historic and current
habitat extent is due, in part, to effects of hysteresis related to positive feedbacks of
seagrass on sediment stabilization (van derHeide etal. 2007}. Van derHeide etal. (2007}
used a simple model relating Z. marina growth rates in the Wadden Sea to light availability
which was affected by background levels of suspended solids (SS} in the water column, an
increase in SS towards the sediment surface reflecting resuspension, and the reduction in
resuspension of SS as Z. marina biomass increases and dampens tidal currents and wave
action. This simple model predicts that seagrass and bare sediment represent alternative
stable states in coastal waters. Seagrass can persist until dramatic reductions due to a
disturbance or disease, after which recovery may be unlikely or impossible in some zones
because of the resuspension of fine sediments in the shallow nearshore zone which reduce
light availability.
Tidal currents estimated for Narragansett Bay are sufficiently strong to resuspend fine
sediments and reduce the probabilities of eelgrass re-establishment following disturbance,
but not strong enough to damage existing eelgrass beds. Oviatt and Nixon (1975}
measured sediment resuspension rates in Narragansett Bay 8-20 times greater than
sediment deposition rates, with resuspension rates greatest 1 meter above the sediment
and decreasing towards the water surface. Likewise, Collins (1976} found an increase in
transparency in upper waters of the Bay from north to south, but a corresponding increase
in a bottom turbidity plume which extended 1 meter above the sediment in shallower
stations to up to 4 meters above the sediments at deeper stations. Tidal currents or
maximum orbital velocities associated with wave action greater than 15 to 30 cm/sec are
sufficient to resuspend the fine sediments (clay, silt and find sand} predominant in much of
Narragansett Bay (Oviatt and Nixon 1975}. Maximum tidal currents during ebb flow in
Narragansett Bay are typically above the 15-30 cm/s range in the mid to lower Bay, except
for selected protected areas (Figure 21}. Note that most thresholds for physical damage to
mature eelgrass by currents are somewhat higher than those required for fine sediment
resuspension, in the range of 40 - 180 cm/sec, beyond the maximum tidal currents
predicted for the bay (Nelson 2009}. Therefore, once established, eelgrass would be
expected to survive throughout much of the bay given adequate light, sediment quality and
protection from wave energy. Maximum wave orbital velocities calculated for the 2009
growing season were much lower, with areas of potential resuspension very close to the
shore and average probability of suspension for fine sand equal to 6.5 % for those areas
(Figure 22}.
55
-------
Legend
•!! Kopp Historical SAV
6 9 12 Kilometers
i i I i i i I
Figure 20. Historic occurrences of seagrass in Narragansett Bay compiled by Kopp (1995).
56
-------
Legend
PredCurrent2010alb
Spd_Ebb
* 0.01-0.15
O 0.15-0.30
• 0.30-2.70
SpEbidw2pt3d
0.01-0.15
0.15-0.3
0.3-2.70
024 8 Kilometers
I i i i I i i i i
Figure 21. Predicted 2010 maximum tidal currents (m/sec) during ebb tides (from NOAA tidal currents web site:
http://tidesandcurrents.noaa.gov/curr_pred.html). Points represent predicted values at tidal stations. Shaded areas were
nterpolated within the 0-5 meter depth zone by inverse distance weighting (2 point interpolation with shorelines as barriers
to interpolation).
57
-------
susp_prob_gs07
Value
High: 100
Figure 22. Probability of fine sediment resuspension based on USGS WAVES model for growing season. Average probability
of resuspension for limited nearshore areas where the average probability greater than zero is 4.9%.
4.4 Model application for assessing management scenarios
Here we demonstrate an application of our model to support adaptive management, with
prediction of the effects of planned reductions in total nitrogen loadings from wastewater
treatment plants (WWTPs] and from atmospheric deposition. Our predictive models
contain four terms related to direct or indirect nitrogen effects on seagrass presence:
density of unsewered residences on high infiltration soils (indicator of groundwater N
inputs], salinity (indicator of N gradient], Secchi depth, and sediment percent total organic
carbon. The density of unsewered residences on high infiltration soils is likely correlated
with groundwater N concentrations (IEC 2012], although there may be a lag in response
due to groundwater travel time. The range of salinity values encompassed by this model is
well within the tolerance ranges reported for Zostera marina populations: 5 to 35 PSU in
the northern hemisphere (den Hartog 1970] and historically, seagrass was found in
58
-------
Narragansett Bay from the lower Providence River south to the mouth of the estuary
(Figure 20}. Thus, it is probable that the modeled positive response to salinity over the
range of 26 to 28.5 PSU in Narragansett Bay represents a response to the inversely
correlated gradient in surface water nitrogen concentrations from head-of-tide to the
mouth of the estuary (Krumholz 2012}. Nitrogen can exert negative effects on seagrass
independent of light reductions from phytoplankton shading, e.g., through stimulatory
effects on epiphytes (Costa 1988} and macroalgae (Teichberg etal 2010} which reduce
light availability and increase anoxia in the sediments. Costa estimated that the depth of
eelgrass growth in Buttermilk Bay (a well-flushed sub embayment within Buzzards Bay}
decreased by 9 cm for every 1 |iM increase in dissolved inorganic nitrogen due to increased
shading by periphyton growth (not phytoplankton}. In Narragansett Bay, an increase in
salinity from 26 to 28.5 ppt corresponds to a reduction in TN from 0.508 to 0.419 mg N/L.
If Costa's relationship for a well-flushed estuary holds in Narragansett Bay, this change in
TN would yield an estimated increase in maximum depth of seagrass of 1.2 meters in
response to reduced periphyton shading.
Model coefficients for Secchi depth capture the response of seagrass to increased light
availability. Following improvements in wastewater treatment in the 1970's, the
transparency of Narragansett Bay improved as TSS loadings from wastewater treatment
plants (WWTPs} decreased (Borkman and Smayda 1988}. Turbidity values are now
uniformly low along the main axis of the bay (Nu Shuttle data;
http://www.narrbay.org/d_projects/nushuttle/shuttletree.htm}. It is likely that most of
the downstream gradient in transparency in the upper water column is related to changes
in chlorophyll a concentration.
Assuming that salinity is an indicator of nitrogen loading gradients related to tidal flushing,
we can also estimate the effect of N load reductions indirectly by converting the salinity
gradient in the study area to the corresponding gradient in total N. Based on an overlay of
nutrient sampling stations from Krumholz (2012} with our salinity grid, we calculated the
relationship between the estuarine gradients in salinity and total N:
TNanavg = 1.461 - 0.0417 Sal (r2 = 0.86}
TNsmravg = 1.431 - 0.0355 Sal (r2 = 0.79}
where TNanavg = annual average total N (mg N/L} for 2006 - 2010 surveys
TNsmravg = summer average total N (mg N/L} for 2006-2010 surveys
Sal = salinity (PSU}
We can estimate the potential effect of N concentration reductions on transparency using a
regional model relating chlorophyll a to the diffuse attenuation coefficient, derived from
data for 48 estuarine sites in southeastern Massachusetts (Bensen et al. 2013}:
Chla = 5.701n[TN] + 10.53
POC = 0.11 [Chi a]+ 0.12
59
-------
K = 1.06 [POC]+0.22,
where
Chi a = chlorophyll a (mg/L]
TN = total nitrogen (mg N/L)
POC = particulate organic carbon (mg/L]
K = attenuation coefficient (nr1)
and the relationship between the light attenuation coefficient and Secchi depth (Poole and
Atkin 1929}:
Kd = 1.7/Secchi depth
Combining these yields the following equation:
Secchi depth = 1/(0.926 + 0.388 ln[TN])
The apparent light compensation point for seagrass varies with sediment percent TOC
(Kenworthy et al. 2014], so estimates of potential seagrass recovery must factor in a lag
time for reductions in sediment TOC. We can estimate a short-term response to increased
light availability by holding sediment TOC constant, and a potential long-term response
assuming reductions in TOC content down to reference levels expected based on percent
silt + clay. We downloaded the raw data from McMaster's (I960] sediment collections in
Narragansett Bay (Poppe et al. 2003], calculated the average percent (silt + clay] for each
of Shepard's sediment classes in the data set (Table 15], and then calculated an expected
reference level of sediment total organic carbon based on Pelletier etal (2011}:
VTOC = 0.11 % silt-clay + 0.556
To simulate potential short-term effects of N load reductions we estimated a projected
change in Secchi depth and a change in salinity "equivalents", the latter to evaluate
potential effects of nutrient reductions beyond transparency effects. We kept sediment
organic carbon levels constant, assuming that these might take some time to recover. To
simulate potential recovery in the long-term, we also applied the model after reducing
sediment organic carbon levels back to reference levels based on sediment particle-size
class.
60
-------
Table 15. Average % silt + clay in McMaster surficial sediment samples from Narragansett Bay by Shepard's sediment class
and estimated reference level of sediment total organic carbon based on Pelletier et al. (2011).
Avg
Shepard_class %siltclay Ref %TOC
Gravel
Sand
Gravelly sediment
Silty sand
Sandy silt
Sand silt clay
Clayey silt
Silt
7.06
9.78
31.35
40.78
68.95
72.92
90.32
96.20
0.40
0.44
0.81
1.01
1.73
1.84
2.40
2.61
We estimated potential reductions in total N loading to Narragansett Bay based on
projected mandated changes in loadings from sewage treatment plants (Krumholz 2012],
urban runoff (as a function of reduced atmospheric loads to impervious surfaces in the
watershed], and reductions in direct atmospheric loadings to the water surface related to
implementation of the Clean Air Act. Atmospheric load reductions were based on the
projected difference in annual loads for the open waters of the estuary between 2006 and
2020, based on results from CMAQ model runs (downloaded from www.epa.gov/edm}. In
2005, RIDEM mandated load reductions from sewage treatment plants at the head of
Narragansett Bay to reduce wastewater N loading to the bay by 50% by 2014 (Krumholz
2012}. After scheduled load reductions from WWTP beyond 2006 and projected decreases
in atmospheric deposition, TN loads could be reduced further by ~40%. However, actual
reductions might be lower than planned if sediment denitrification rates continue to
decrease (Krumholz 2012}. Based on the equations above, we would expect a 40%
reduction in TN concentrations to result in a 2.2 meter increase of an initial Secchi depth of
2.4 meters to 4.6 meters. (Specific increases will vary by initial Secchi depth because of the
inverse relationship.} A reduction in TN concentrations would be equivalent to an increase
in salinity of 16.6 to 19 percent in our model; we used a conservative estimate of 16.6%.
Prediction of seagrass presence requires that we choose a threshold probability level.
Many modelers choose a default level of 50% (odds ratio of 1} as an indicator of the
minimum level at which seagrass is expected to exist. Depending on the relative cost of
false positive versus false negative errors, users may choose alternative thresholds
(Fielding and Bell 1997}. Van derHeide et al. (2007} estimated the minimum density of
seagrass required for positive feedback effects promoting sediment stability. For the
Wadden Sea this corresponded to a threshold value of 30% maximum biomass. We
present estimated changes in seagrass coverage based on a range of alternative thresholds.
4.5 Model predictions of seagrass increase with decreased nitrogen loading
There are physical constraints in the amount of potential seagrass habitat in Narragansett
Bay. Based on the current model, the colonized area for all shorelines combined following
a 40% reduction in TN loads (and concentration} would increase from 12% of area in the 0
to 5 meter depth zone to about 63% of area in the short term and slightly more in the long
term (as sediment organic carbon levels recovered} (Figure 23a}. Given a threshold for
61
-------
presence/absence of 0.5, we predict recovery potential in the short term to differ
substantially among shorelines (Figure 23b]. Long term recovery (assuming return of
sediment organic carbon to reference levels] is estimated to be virtually complete for the
majority of the most favorable shorelines, but negligible for Shoreline 17 (yellow lines,
Figure 23c]. These sensitivity analyses are approximate given that background
disturbance levels might change over time and that we haven't factored in the positive
feedbacks in sediment stabilization due to seagrass growth. It is possible that specific
restoration measures such as co-restoration of shellfish beds (to reduce suspended
sediments and effects of tidal currents and wave action] or use of existing or constructed
coastal barriers to limit effects of wave action and tides would improve probabilities of
initial colonization success and the initiation of positive feedback effects (Bos and van
Katwijk 2007]. Thus, we can use our model to project ecologically significant increases in
seagrass coverage following planned reductions in total N loading based on effects of
nitrogen on light availability mediated by phytoplankton and periphyton growth.
4.6 Future improvements
Although our prediction of seagrass absolute or average presence/absence at shoreline
locations was very robust based on our current model, our models could be improved in
the future with the incorporation of finer resolution data and/or more sophisticated
modeling approaches. Accuracy of depth limits could probably be improved with better
resolution of digital elevation models (DEMs]. NOAA's National Geophysical Data Center
(NGDC] is building high-resolution Tsunami Inundation DEMs of select U.S. coastal regions.
These DEMs are referenced to a vertical tidal datum of NAVD 88 or Mean High Water
(MHW] and horizontal datum of World Geodetic System of 1984 (WGS 84]. Cell sizes will
range from 1/3 arc-second (-10 meters] to 36 arc-seconds (~1 km]
(http://www.ngdc.noaa.gov/mgg/inundation/tsunami/inundation.html]. Although not an
issue for Narragansett Bay, the improved resolution between the shoreline (mean sea level,
MSL] and MLLW (currently set to zero in merged topobathymetry grids] will be
particularly important for predicting distributions of intertidal populations. In the near
future, we will have access to maps of predicted tidal currents throughout Narragansett
Bay as part of a recently developed hydrodynamic model (Abdelrhman 2015]. An
accompanying water quality model will also include finer scale predictions of nitrogen
concentrations throughout the estuary, so that we will not have to rely on salinity as an
indicator of total nitrogen gradients (US EPA 2015].
62
-------
a)
100%
75%
§1
1 I 5°%
So
25%
b)
0%
Simulation 1
Simulation 2
..i -
0.00
100%
75%
x
c
CD
§Ji
•ft ro
So
cu o
D- _
to
50%
25%
0.25 0.50
X
0.75
1.00
0.00
1.00
100%
75%
50% -
(D O
0 --
25%
0% -"•
0.00
0.25
0.50
X
0.75
1.00
Figure 23. Projected change in cumulative distribution function for probability of occurrence (x) following 40% reduction in
total N without (simulation 1) and with (simulation 2) sediment recovery, a) All shorelines combined, b) Projected short-
term effect of 40% TN reduction for 8 most favorable shorelines (favorability in order of red - orange - yellow - light green -
dark green - light blue - dark blue - indigo) with solid lines indicating condition before TN reductions and dashed lines
indicating projected condition after TN reductions, c) Long term recovery following 40% reduction in total N for 8 most
favorable shorelines assuming recovery of sediment percent organic carbon to reference levels.
A bio-optical model is also under development for Narragansett Bay (US EPA 2015], which
could improve model inputs describing light availability. However, the biggest limitation to
model predictions of light availability to seagrass is the lack of information on the
resuspension and transport of fine particulates in the nearshore zone. Most measurements
of light attenuation (either profiles or Secchi depth measurements] have been made in
deeper waters along the main axis of the estuary. Given that the Secchi depth is less than
63
-------
the depth of highest turbidity associated with resuspended sediments near the seabed4,
S.D. values will represent only background light attenuation and may be of limited use for
predicting light availability in seagrass beds or within potential seagrass habitat prior to
recovery. Time series of light measurements from underwater HOBOs in areas of potential
seagrass habitat both inside and outside of existing seagrass beds are needed, similar to
recordings made for Massachusetts coastal bays (Kenworthy et al. 2014}.
We would have had access to more methods for incorporating spatial autocorrelation into
predictive models if we were not constrained by memory requirements. It is possible that
use of parallel computing methods (e.g., simultaneous use of multiple CPUs in an existing
quad core processor] through application of R packages such as snow (Tierney et al. 2014}
would alleviate this issue. However, we would still have to deal with the limitation of
existing packages for methods such as regression kriging in handling anisotropy in spatial
autocorrelation.
We could not include all potential interaction terms in our models because of the high
variance inflation factors generated. Boosted regression trees (BRTs) are a robust
approach for handling models with nonlinear effects and multiple interactions. Of all those
approaches tested, Valle et al. (2013} found the best performance for BRT models.
Bayesian modeling approaches also could be useful in the future to facilitate adaptive
management of habitat (March et al. 2013}. Ultimately, managers may need to develop
coupled hydrodynamic/water quality/sediment diagenesis models with feedback effects to
capture the hysteresis and lags inherent in eelgrass decline and recovery (Eldridge and
Morse 2000, van derHeide et al. 2007, delBarrio et al. 2014, Kenworthy et al. 2014}.
Our model, like any model, could potentially be improved with additional data or enhanced
analyses. Nonetheless, we believe that our predictions are sufficiently robust to support
decision-making in an adaptive management framework. Potential applications of our
model include identification of ALU zones for setting nutrient criteria for areas of potential
seagrass habitat, prioritization of areas and strategies for seagrass restoration, and
projection of potential benefits of management actions. Even with the potential for model
improvements discussed above, our model represents a significant advancement over
existing models which focus on a limited set of factors influencing seagrass growth, fail to
incorporate and correct for spatial autocorrelation, and which fail to incorporate
interaction of environmental variables or nonlinear effects.
4 This turbidity maximum is based on depth profiles, and should not be confused with the "turbidity maximum"
that occurs in the upper estuary reaches of some estuaries.
64
-------
Chapter 5. Conclusions
Our pilot project addressed the multiple challenges associated with development of species
distribution models for seagrass: the fine-scale patchiness of seagrass distributions with
attendant problems of spatial autocorrelation, the large areas of interest for model
development and application entailing significant memory demands for modeling, and the
potential co-variance of multiple interacting factors affecting seagrass. The fine-scale of
data necessary to describe the patchy nature of seagrass distributions coupled with the
large model application area for management decisions create significant memory
demands which can limit the number and type of R packages that can be applied in
practice. However, we found that addition of a residual autocorrelation term in logistic
regression models, as suggested by Crase et al. (2013} virtually eliminated the presence of
spatial autocorrelation in model residuals. While Crase et al. limited their residual
autocorrelation term to zonal averages calculated among adjacent grid cells, we expanded
the focal average to cover the range of spatial autocorrelation evident from correlograms.
The use of a rectangular zonal average allowed us to account for anisotropy in spatial
errors, with a range of up to 1320 meters in the axis parallel to the shoreline but only 200
meters in the offshore direction. As predicted, incorporation of a spatial autocorrelation
term in regression models reduced the number of significant variables included, generally
reducing the number of distinct sediment type responses and/or interactions detected.
We were able to deal with most but not all of the issues related to cross-correlation of
potential explanatory variables. Centering variables prior to incorporating them into
models not only kept variance inflation factors values low but also aided in interpretation
of model coefficients. Both light availability and wave energy co-vary with depth, and both
were significant explanatory variables. We were able to include both types of variables in
our models by expressing light availability effects based on the difference between Secchi
depth and seagrass bed depth and representing the effect of wave mixing depth as a binary
variable (depth > wave mixing depth}. Gradients that co-vary with salinity along the main
axis of the estuary are more problematical. Although temperature is potentially important
in affecting seagrass populations, we excluded it because the populations we were dealing
with were subtidal and the range of temperatures measured in our well-flushed system
was lower than those at which effects had been observed. We did include salinity which
covaries with nutrient concentrations along the axis of the estuary, but given the low range
of values, it is probably serving as a proxy for nutrient concentration effects not captured
by changes in transparency related to phytoplankton biomass. However, the co-varying
gradients of salinity and nutrients would make it difficult to model interactions of nutrients
and salinity, although these have been demonstrated experimentally (van Katwijk et al.
1999}.
We predicted seagrass distribution at the scale of 10 x 10-meter grid cells, as
presence/absence or average presence/absence associated with shoreline locations spaced
at 10-meter intervals, and minimum or maximum depth of distributions at those locations.
Prediction of seagrass absolute or average presence/absence at shoreline locations was
65
-------
very robust, with area-under-the-curve (AUC] values associated with Receiver Operating
Characteristic (ROC] curves of 0.95 - 0.98 following 10-fold cross-validation of models.
Although the model predicting shoreline presence was the most robust of those tested
across different scales in Narragansett Bay, other scales of resolution might work better for
other types of estuaries. In Narragansett Bay, random shoreline effects varied over several
orders of magnitude, probably tied to the distribution of tidal currents which are weak
enough to allow persistence of existing seagrass beds, but strong enough to resuspend fine
sediments that interfere with successful recolonization. For the model predicting seagrass
presence/absence at the grid cell scale, the most influential predictor of fixed effects is
Secchi depth, followed by, in order: shoreline isolation, sediment percent total organic
carbon, sediment type, and salinity. The least influential variable is water depth greater
than average wave mixing depth. For the model predicting presence of seagrass at
shoreline locations, the most influential predictor is sediment type, followed by sediment
percent total organic carbon (at low Secchi depth], then salinity (as a proxy for water
column total nitrogen}.
Multiple modes of action for nutrients can be simultaneously incorporated into empirical
models for seagrass distribution. We were able to capture the effects of nutrients on light
availability, other nutrient impacts potentially related to stimulation of periphyton and
macro-algal growth (using salinity as a proxy for total N], and longer term impacts related
to the accumulation of organic matter in the sediments. The latter was reflected in a
significant interaction between the light compensation point (maximum depth] for
seagrass and sediment percent organic carbon.
Application of our model in predictive mode suggests that different shorelines in
Narragansett Bay may have very different recovery potentials, and that interventions to
reduce tidal energy and/or sediment resuspension may be needed as part of restoration
activities. Longer-term recovery potential, related to recovery of sediment organic carbon
to levels characteristic of reference condition in New England estuaries, is greater than
predicted short-term recovery potential related to reductions in total nitrogen
concentrations in the water column and increased transparency related to reductions in
phytoplankton biomass.
Many of the georeferenced data layers required to parameterized species distribution
models of this type are publically available online through EPA's Estuary Data Mapper
application (Detenbecketal. 2009; www.epa.gov/edm], although spatial resolution of
some remotely sensed indicators of water quality for estuaries are limited. In our
application, we made substitutions in some cases, based on spatially intensive monitoring
that had been carried out in Narragansett Bay. Publically available monitoring data for
nutrients are limited to several stations, which tend to be located along the main axis of the
bay, are extremely limited in the Sakonnet River arm of the system, and are nonexistent for
shallow waters of the bay. Tidal current measurement data are limited for most estuaries
and predicted values have only been made publically available online in beta form by NOAA
very recently. Finer-scale predictions yielded by future publically available models may be
critical in identifying the best protected locations for eelgrass restoration. Our model could
be improved by better characterization of optical parameters, particularly in the near-
66
-------
benthic environment and in the shallow waters in which seagrass is typically found. Long-
term records of Secchi depth do not capture trends in the resuspension of fine sediments in
the bottom waters which may be limiting seagrass recolonization.
In spite of the methodological challenges and data limitations faced in developing a
predictive statistical model for seagrass distribution, our model performed well and was
particularly robust for predictions of shoreline occurrence. Our models were also
successful in elucidating the effect of multiple interacting stressors in determining the
distribution of seagrass patches in Narragansett Bay, as well as factors which might limit
recovery. Finally, we were able to demonstrate multiple pathways for nutrient effects on
seagrass growth and survival in the bay, including long-term effects of eutrophication on
sediment organic carbon which could persist for decades and slow final recovery.
67
-------
References
Abdelrhman, M.A. Three-dimensional Modeling of Hydrodynamics and Transport in
Narragansett Bay. EPA/600/R-15/152. Washington, DC: U.S. Environmental Protection
Agency, 2015.
Batiuk, R. A., R. J. Orth, K. A. Moore, W. C. Dennison, J. C. Stevenson, L. W. Staver, V. Carter, N. B.
Rybicki, R. E. Hickman, S. Kollar, S. Bieber, and P. Heasly. Chesapeake Bay Submerged
Aquatic Vegetation Habitat Requirements and Restoration Targets: A Technical Synthesis,
CBP/TRS 83/92. Annapolis, MD: U.S. Environmental Protection Agency, 1992.
Batiuk, R., P. Bergstrom, M. Kemp, E. Koch, L. Murray, C. Stevenson, R. Bartleson, V. Carter, N.
Rybicki, J. Landwehr, C. Gallegos, L. Karrh, M. Naylor, D. Wilcox, K. Moore, S. Ailstock, and M.
Teichberg. Chesapeake Bay Submerged Aquatic Vegetation Water Quality and Habitat-Based
Requirements and Restoration Targets: A Second Technical Synthesis, CBP/TRS 245/00.
EPA/903/R-00/014. Annapolis, MD: U.S. Environmental Protection Agency, Chesapeake
Bay Program, 2000. http://archive.chesapeakebay.net/pubs/sav/index.html.
Benson, J. L., D. Schlenzinger, and B.L. Howes. Relationship between nitrogen concentration,
light, and Zostera marina habitat quality and survival in southeastern Massachusetts
estuaries. Journal of Environmental Management 131: 129-137 (2013}.
Borkman, D. G., and T. J. Smayda. Long-term trends in water clarity revealed by Secchi-disk
measurements in lower Narragansett Bay. ICES Journal of Marine Science 55: 668-679
(1988}.
Bos, A. R., and M. M. van Katwijk. Planting density, hydrodynamic exposure and mussel beds
affect survival of transplanted intertidal eelgrass. Marine Ecology Progress Series 336: 121-
129(2007}.
Bradley, M., K. Raposa, and S. Tuxbury. Report on the Analysis of True-color Aerial Photography
to Map and Inventory Zostera marina L. in Narragansett Bay and Block Island, Rhode Island.
Page 1-16 and 9 Mapsheets. Kington, RI: Rhode Island Natural History Survey, 2007.
Burkholder, J. M., D. A. Tomasco, and B. W. Touchette. Seagrasses and eutrophication./oi/rna/ of
Experimental Marine Biology and Ecology 350: 46-72 (2007}.
Burnham, K. P., and D. R. Anderson. Model Selection and Multimodel Inference: A Practical
Information-Theoretic Approach. 2nd edition. New York, NY: Springer, 2002.
Chambers, P. A. Nearshore occurrence of submersed aquatic macrophytes in relation to wave
action. Canadian Journal of Fisheries and Aquatic Sciences 44: 1666-1669 (1987}.
Collins, B.P. Suspended material transport: Narragansett Bay area, Rhode Island. Estuarine and
Coastal Marine Science 4: 33-44 (1976}.
68
-------
Compton, J. E., J. Harrison, R. L. Dennis, T. L. Greaver, B. H. Hill, S. J. Jordan, H. Walker, and H. V.
Campbell. Ecosystem services altered by human changes in the nitrogen cycle: A new
perspective for US decision making. Ecology Letters 14: 804-815 (2011}.
Costa, J. E. Distribution, production, and historical changes in abundance of eelgrass (Zostera
marina 1.} in southeastern Massachusetts. Ph.D. thesis, Boston University, Boston, MA,
1988.396pp.
Costa, J. E., B. L. Howes, and E. Gunn. Report of the Buzzards Bay Citizens' Water Quality
Monitoring Program 1992-1995. New Bedford, MA: Coalition for Buzzards Bay, 1996.
Crase, B., A. C. Liedloff, and B. A. Wintle. A new method for dealing with residual spatial
autocorrelation in species distribution models. Ecography 35: 879-888 (2012}.
De Boer, W. F. Seagrass-sediment interactions, positive feedbacks and critical thresholds for
occurrence: A review. Hydrobiologia 591: 5-24 (2007}.
Dejonge, V. N., D. J. de Jong, and M. M. van Katwjik. Policy plans and management measures to
restore eelgrass (Zostera marina L.} in the Dutch Wadden Sea. Helgoland Marine Research
54:151-158(2000}.
Del Barrio, P., N. K. Ganju, A. L. Aretxabaleta, M. Hayn, A. Garcia, and R. W. Howarth. Modeling
future scenarios of light attenuation and potential seagrass success in a eutrophic estuary.
Estuarine, Coastal and Shelf Science, 149: 13-23 (2014}.
denHartog, C. The Seagrasses of'the World. Amsterdam: North-Holland Publication, 1970. 275
pp.
Dennison, W. C. Effects of light on seagrass photosynthesis, growth and depth distribution.
Aquatic Botany 27: 15-26 (1987}.
Dennison, W. C., R. J. Orth, K. A. Moore, J. C. Stevenson, V. Carter, S. Kollar, P. W. Bergstrom, and
R. A. Batiuk. Assessing water quality with submersed aquatic vegetation. Bioscience 43: 86-
94(1993}.
Detenbeck, N., M. Pelletier, M. ten Brink, M. Abdelrhman, and S. Rego. "E-Estuary: Developing a
decision-support system for coastal management in the conterminous United States." In
Proceedings, 33rd IAHR 2009 Congress - Water Engineering for a Sustainable Environment.
Vancouver, BC, 2009. Pp. 5284-5292.
Downie, A-L, M. von Numers, and C. Bostrom. Influence of model selection on the predicted
distribution of the seagrass Zostera marina. Estuarine, Coastal and Shelf Science 121/122: 8-
19(2013}.
Duarte, C. M. Seagrass depth limits. Aquatic Botany 40: 363-377 (1991}.
69
-------
Eldridge, P. M., and J. W. Morse. A diagenetic model for sediment-seagrass interactions. Marine
Chemistry 70: 89-103 (2000}.
Fielding, A. H. and J. F. Bell. A review of methods for the assessment of prediction errors in
conservation presence/absence models. Environmental Conservation, 24: 38-49 (1997}.
Fonseca, M. S., and A. Malhotra. WEMo (Wave Exposure Model) User Manual, Version 4.0.
Beaufort, NC: NOAA National Centers for Coastal Ocean Science, Center for Coastal
Fisheries and Habitat Research. March 2010.
Fonseca M. S., W. J. Kenworthy, and G. W. Thayer. Guidelines for Conservation and Restoration of
Seagrass in the United States and Adjacent Waters. NOAA/NMFS Coastal Ocean Program
Decision Analysis Series, No. 12. Silver Spring, MD: NOAA Coastal Ocean Office, 1998.
Fonseca, M. S., P. E. Whitfield, W. J. Kenworthy, D. R. Colby, and B. E. Julius. Use of two spatially
explicit models to determine the effect of injury geometry on natural resource recovery.
Aquatic Conservation: Marine and Freshwater Ecosystems 14: 281-98 (2004}.
Garcia, E., T. C. Granata, and C. M. Duarte. An approach to measurement of particle flux and
sediment retention within seagrass (Posidonia oceanica} meadows. Aquatic Botany 65:
255-268(1999}.
Goodman, J. L, K. A. Moore, and W. C. Dennison. Photosynthetic responses of eelgrass (Zostera
marina L.} to light and sediment sulfide in a shallow barrier island lagoon. Aquatic Botany
50: 37-47 (1995}.
Grech, A., and R. G. Coles. An ecosystem-scale predictive model of coastal seagrass distribution.
Aquatic Conservation: Marine and Freshwater Ecosystems 20: 437-444 (2010}.
Hammar-Klose, and E. R. Thieler. Coastal Vulnerability to Sea-Level Rise: A Preliminary
Database for the U.S. Atlantic, Pacific, and Gulf of Mexico Coasts. U.S. Geological Survey,
Digital Data Series DDS-68,1 CD-ROM, 2001.
Harlin, M. M., and B. Thorne-Miller. Nutrient enrichment of seagrass beds in a Rhode Island
coastal lagoon. Marine Biology 65: 221-229 (1981}.
Holmer, M. and E. J. Bondgaard. Photosynthetic and growth response of eelgrass to low oxygen
and high sulfide concentrations during hypoxic events. Aquatic Botany 70: 29-38 (2001}.
Howarth, R. W. Nutrient limitation of net primary production in marine ecosystems. Annual
Review of Ecology and Systematics 19: 89-110 (1988}.
Howes, D., M. Morris, and M. Zacharias. British Columbia Estuary Mapping System. Version 1.
Province of British Columbia, Canada: Land Use Coordination Office for the Coastal Task
Force, Resource Inventory Committee Resources Inventory Committee, 1999. Downloaded
from: http://www.for.gov.bc.ca/ric.
70
-------
IEC. "Appendices", In Narragansett Bay Sustainability Pilot: Phase I Report, 2012.
http://www.epa.gov/research/docs/nbsp-phase-I-report-appendices.pdf.
Jarvis, J. C., and K. A. Moore. The role of seedlings and seed bank viability in the recovery of
Chesapeake Bay, USA, Zostera marina populations following a large-scale decline.
Hydrobiologia 649: 55-68 (2010}.
Keddy, C. J., and D. G. Patriquin. An annual form of eelgrass in Nova Scotia. Aquatic Botany 5:
163-170 (1978}.
Keith, D., B. Schaeffer, R. Lunetta, R. Gould, Jr., K. Rocha, and D. Cobb. Remote sensing of
selected water-quality indicators with the Hyperspectral Imager for the Coastal Ocean
(HICO} Sensor. International Journal of Remote Sensing 35: 2927-2962 (2014}.
Kelly, N. M., M. Fonseca, and P. Whitfield. Predictive mapping for management and
conservation of seagrass beds in North Carolina. Aquatic Conservation: Marine and
Freshwater Ecosystems 11: 437-451 (2001}.
Kendrick, G. A., C. M. Duarte, and N. Marba. Clonality in seagrasses, emergent properties and
seagrass landscapes. Marine Ecology Progress Series 290: 291-6 (2005}.
Kenworthy, W. J., and M. Fonseca. Reciprocal transplant of the seagrass Zostera marina L. Effect
of substrate on growth. Aquaculture 12: 197-213 (1977}.
Kenworthy, W. J., C. L. Gallegos, C. Costello, D. Field, and G. di Carlo. Dependence of eelgrass
(Zostera marina] light requirements on sediment organic matter in Massachusetts coastal
bays: Implications for remediation and restoration. Marine Pollution Bulletin 83: 446-457
(2014}.
Koch, E. W. Beyond light: Physical, geological, and geochemical parameters as possible
submersed aquatic vegetation habitat requirements. Estuaries 24: 1-17 (2001}.
Koch, M. S., S. A. Schopmeyer, 0.1. Nielsen, C. Kyhn-Hansen, and C. J. Madden. Conceptual model
of seagrass die-off in Florida Bay: Links to biogeochemical processes. Journal of
Experimental Marine Biology and Ecology 350: 73-88 (2007}.
Kopp, B .S. 1848-1994 Historical Eelgrass Beds of Narragansett Bay, 1995. Downloaded from
www.edc.uri.edu/eelgrass.
Kopp, B. S., A. M. Doherty, and S. W. Nixon. A Guide to Site Selection for Eelgrass Restoration
Projects In Narragansett Bay, Rhode Island. Technical Report to the Narragansett Bay
Project. Narragansett, RI: University of Rhode Island, 1994.
Krause-Jensen, D., M. F. Pedersen, and C. Jensen. Regulation of eelgrass (Zostera marina] cover
along depth gradients in Danish coastal waters. Estuaries 26: 866-977 (2003}.
Krumholz, J.S. Spatial and temporal patterns in nutrient standing stock and mass-balance in response
to load reductions in a temperate estuary. Ph.D. Thesis, University of Rhode Island, 2012.
71
-------
Latimer, J. S. and S. A. Rego. Empirical relationship between eelgrass extent and predicted
watershed-derived nitrogen loading for shallow New England estuaries. Estuarine, Coastal
and Shelf Science 90: 231-240 (2010}.
Lee, K., S. R. Park, and Y. K. Kim. Effects of irradiance, temperature, and nutrients on growth
dynamics of seagrasses: A review. Journal of Experimental Marine Biology and Ecology 350:
144-175 (2007}.
Maggini, R., A. Lehmann, N. E. Zimmermann, and A. Guisan. Improving the generalized
regression analysis for the spatial prediction of forest communities. Journal of
Biogeography 33: 1729-1749 (2010}.
Marba, N. and C. M. Duarte. Rhizome elongation and seagrass clonal growth. Marine Ecology
Progress Series 174: 269-280 (1998}.
March, D., J. Alos, M. Cabanellas-Reboredo, E. Infantes, A. Jordi, and M. Palmer. A Bayesian
spatial approach for predicting seagrass occurrence. Estuarine, Coastal and Shelf Science
131:206-12(2013}.
McMaster, R. L. Sediment Types, Narragansett Bay Sediments; A map ofNarragansettBay
sediments from McMaster, 1960, Published in 1993,1:24000 (lin=2000ft} scale, State of
Rhode Island and Providence Plantations, I960. Downloaded from
http://www.narrbay.org/d_downloads/D_Biological/D_habitat/sediment.htm#6.
Murphey P. L., and M. S. Fonseca. Role of high and low energy seagrass beds as nursery areas
for Penaeus duorarum in North Carolina. Marine Ecology Progress Series 121: 91-98 (2005}.
Neckles, H. A., F. T. Short, S. Barker, and B. S. Kopp. Disturbance of eelgrass Zostera marina by
commercial mussel Mytilus edulis harvesting in Maine: Dragging impacts and habitat
recovery. Marine Ecology Progress Series 285: 57-73 (2005}.
Nelson, W.G. (ed.} Seagrasses and Protective Criteria: A Review and Assessment of Research
Status, EPA/600/R-09/050. Washington, DC: U.S. Environmental Protection Agency, 2009.
Nielsen, S. L., J. Borum, 0. Geertz-Hansen, and J. Sand-Jensen. Marine bundplanters
dybdegraense. Vand Miljo 5: 217-220 (1989}.
Nixon, S. W., and M. E. Pilson. "Nitrogen in estuarine and coastal marine systems." In Nitrogen
in the Marine Environment, eds. E. J. Carpenter, and D. G. Capone, 565-648. New York:
Academic Press, 1983.
Oakley, B. A., J. D. Alvarez., and J. C. Boothroyd. Benthic geologic habitats of shallow estuarine
environments: Greenwich Bay and Wickford Harbor, Narragansett Bay, Rhode Island, U.S.A.
Journal of Coastal Research 28: 760-773 (2012}.
72
-------
Oczkowski, A., S. Nixon, K. Henry, P. DiMilla, M. Pilson, S. Granger, B. Buckley, C. Thornber, R.
McKinney, and J. Chaves. Distribution and trophic importance of anthropogenic nitrogen in
Narragansett Bay: An assessment using stable isotopes. Estuaries and Coasts 31: 53-69
(2008}.
Orth, R. J., M. Luckenbach, and K. A. Moore. Seed dispersal in a marine macrophyte:
Implications for colonization and restoration. Ecology 75: 1297-1939 (1994}.
Orth, R. J., T. J. B Carruthers, W. C. Dennison, C. M. Duarte, J. W. Fourqurean, K. L. Heck, Jr., A. R.
Hughes, G. A. Kendrick, W. J. Kenworthy, S. Olyarnik, F. T. Short, M. Waycott, and S. L.
Williams. A global crisis for seagrass ecosystems. Bioscience 56: 987 - 996 (2006}.
Oviatt, C. A., and S. W. Nixon. Sediment resuspension and deposition in Narragansett Bay.
Estuarine and Coastal Marine Science 3: 201-217 (1975}.
Oviatt, C., P. Doering, B. Nowicki, L. Reed, J. Cole, and J. Frithsen. An ecosystem level experiment
on nutrient limitation in temperate coastal marine environments. Marine Ecology Progress
Series 116: 171-179 (1995}.
Pelletier, M. C., D. E. Campbell, K. T. Ho, R. M. Burgess, C. T. Audette, and N. E. Detenbeck. Can
sediment total organic carbon and grain size be used to diagnose organic enrichment in
estuaries? Environmental Toxicology and Chemistry 30: 538-547 (2011}.
Pilson, M. E. Q. On the residence time of water in Narragansett Bay. Estuaries 8: 2-14 (1985}.
Poole, H.H. and W.R.G. Atkins. Photo-electric measurements of the penetration of light into sea
water. Journal of the Marine Biological Association U.K. 16: 297 (1929}.
Poppe, L. J., V. F. Paskevich, S. J. Williams, M. E. Hastings, J. T. Kelley, D. F. Belknap, L. G. Ward, D.
M. FitzGerald, and P. F. Larsen. Surficial Sediment Data from the Gulf of Maine, Georges Bank,
and Vicinity: A CIS Compilation, U.S. Geological Survey Open-File Report 03-001. Woods
Hole, MA: U.S. Geological Survey Coastal and Marine Geology Program, Woods Hole Field
Center, 2003.
R Core Team. R: A language and environment for statistical computing. Vienna, Austria: R
Foundation for Statistical Computing, 2014. http://www.R-project.org/.
Ralph, P. J., M. J. Durako, S. Enriquez, C. J. Collier, and M. A. Doblin. Impact of light limitation on
seagrasses. Journal of Experimental Marine Biology and Ecology 350: 176-193 (2007}.
RIGIS. Marinas of Rhode Island: s44fsm96. Rhode Island Geographic Information System
Database, 1996 (s44fsm96 accessed from
http://www.edc.uri.edu/rigis/data/data.aspx?ISO=structure August 2014}.
73
-------
RIGIS. Eelgrass Beds in Rhode Island (Polygon}: s44nep99. Rhode Island Geographic
Information System Database, 2000a (s44nep99 accessed from
http://www.edc.uri.edu/rigis/data/data.aspx?ISO=biota August 2014}.
RIGIS. Eelgrass Beds in Rhode Island (Line}: s44ne!99. Rhode Island Geographic Information
System Database, 2000b (s44ne!99 accessed from
http://www.edc.uri.edu/rigis/data/data.aspx?ISO=biota August 2014}.
RIGIS. Hardened Shorelines in Narragansett Bay: hrdshore. RI Dept. of Environmental
Management, Narragansett Bay Estuary Program, 2003. (hrdshore accessed from
http://www.edc.uri.edu/rigis/data/data.aspx?ISO=structure August 2014}.
RIGIS. Submerged Aquatic Vegetation (SAV} in Rhode Island Coastal Waters (2006}; SAV06.
Kingston, RI: Rhode Island Geographic Information System (RIGIS} Data Distribution
System, 2013, Kingston, Rhode Island: Environmental Data Center, University of Rhode
Island (http://www.edc.uri.edu/rigis/data/data.aspx?ISO=biota, last date accessed: 1 May
2013}.
Rivers, D. 0., and F. T. Short. Effect of grazing by Canada geese Branta canadensis on an
intertidal eelgrass Zostera marina meadow. Marine Ecology Progress Series 333: 271-279
(2007}.
Roberts J. ]., B. D. Best, D. C. Dunn, E. A. Treml, and P. N. Halpin. Marine Geospatial Ecology
Tools: An integrated framework for ecological geoprocessing with ArcGIS, Python, R,
MATLAB, and C++. Environmental Modelling & Software 25: 1197-1207 (2010}.
Rohweder, J., J. T. Rogala, B. L. Johnson, D. Anderson, S. Clark, F. Chamberlin, D. Potter, and K.
Runyon. Application of Wind Fetch and Wave Models for Habitat Rehabilitation and
Enhancement Projects-2012 Update. Contract report prepared for U.S. Army Corps
Engineers' Upper Mississippi River Restoration - Environmental Management Program,
2012. 52p.
Schonert, T., and P. Milbradt. Hydro- and morphodynamic simulation considering ecological model
components on asynchronous parallel cellular automata. Computing in Civil Engineering 2005:
1-11 (2005).
Shipman, H. "The geomorphic setting of Puget Sound: Implications for shoreline erosion and
the impacts of erosion control structures" In Puget Sound Shorelines and the Impacts of
Armoring—Proceedings of a State of the Science Workshop, U.S. Geological Survey Scientific
Investigations Report 2010-5254, by Shipman, H., M. N. Dethier, G. Gelfenbaum, K. L. Fresh,
and R. S. Dinicola (eds.}, May 2009, 2010. p. 19-34.
Short, F. T. Effects of sediment nutrients on seagrasses, literature review and mesocosm
experiment. Aquatic Botany 27: 41-57 (1987}.
Short, F. T., L. K. Muehlstein, and D. Porter. Eelgrass wasting disease: Cause and recurrence of a
marine epidemic. Biological Bulletin 173: 557-62 (1987}.
74
-------
Short, F. T., B. W. Ibelings, and C. Den Hartog. Comparison of a current eelgrass disease to the
wasting disease of the 1930's. Aquatic Botany 30: 295-304(1988}.
Short, F. T., D. Burdick, J. Wolf, and G. E. Jones. Eelgrass in Estuarine Research Reserves Along the
East Coast, USA. Part I: Declines from Pollution and Disease and Part II: Management of
Eelgrass Meadows. Durham, NH: NOAA, Coastal Ocean Program Publ., 1993.
Short, F. T. The Port of New Hampshire Interim Mitigation Success Assessment Report. Report to
the New Hampshire Department of Transportation. Durham, NH: Jackson Estuarine
Laboratory, 1993.
Short, F. T., and D. M. Burdick. Eelgrass Restoration Site Selection Model V.1.2. CDROM, 2005.
Smayda, T. J., and D. G. Borkman. "Nutrient and plankton dynamics in Narragansett Bay", Ch. 15
In Science for Ecosystem-based Management, By Desbonnet, A. and B. A. Costa-Pierce (eds.},.
New York, NY: Springer, 2008.
Teichberg, M., S. E. Fox, Y. S. Olsenz, I. Valielaw, P. Martinetto, 0. Iribarnes, E. Y. Muto, M. A. V.
Petti, T. N. Corbisier, M. Soto-Jimenez, F. Paez-Osuna, P. Castro, H. Freitas, A. Zitelli, M.
Cardinaletti, and D. Tagliapietra. Eutrophication and macroalgal blooms in temperate and
tropical coastal waters: nutrient enrichment experiments with Ulva spp. Global Change
Biology 16: 2624-2637 (2010}.
Terrados, J., C. M. Duarte, L. Kamp-Nielsen, N. S. R. Agawin, E. Gacia, D. Lacap, M. D. Fortes, J.
Borum, M. Lubanski, and T. Greve. Are seagrass growth and survival constrained by the
reducing conditions of the sediment? Aquatic Botany 65: 175-197 (1999}.
Thome-Miller, B., M. M. Harlin, G. B. Thursby, M. M. Brady-Campbell, and B. A. Dworetzky.
Variations in the distribution and biomass of submerged macrophytes in five coastal
lagoons in Rhode Island, USA. Botanica Marina XXVI:231-242 (1983}.
Touchette, B. W., and J. M. Burkholder. Seasonal variations in carbon and nitrogen constituents
in eelgrass (Zostera marina L.} as influenced by increased temperature and water column
nitrate. Botanica Marina 45: 23-34 (2002}.
Touchette, B.W. Seagrass-salinity interactions: Physiological mechanisms used by submersed
marine angiosperms for a life at sea. Journal of Experimental Marine Biology and Ecology
350:194-215(2007}.
U.S. Environmental Protection Agency. Ambient Water Quality Criteria for Dissolved Oxygen,
Water Clarity and Chlorophyll a for the Chesapeake Bay and Its Tidal Tributaries.
EPA/903/R/03/002. Annapolis, MD: U.S. Environmental Protection Agency, Region 3,
Chesapeake Bay Program Office, 2003.
U.S. Environmental Protection Agency. Quantitative Models Describing Past and Current
Nutrien t Fluxes and Associated Ecosystem Level Responses in the Narragansett Bay
Ecosystem. EPA/600/R-15/174. Washington, DC: U.S. Environmental Protection Agency,
2015.
75
-------
Valiela, I., S. Mazzilli, J. Bowen, K. Kroeger, M. Cole, G. Tomasky, and T. Isaji. ELM, an estuarine
nitrogen loading model: Formulation and verification of predicted concentrations of
dissolved inorganic nitrogen. Water Air and Soil Pollution 157: 365-391 (2004}.
Valle, M., M. M. van Katwijk, D. J. de Jong, T. J. Bouma, A. M. Schipper, G. Chust, B. M. Benito, J. M.
Garmendia, and A. Borja. Comparing the performance of species distribution models of
Zostera marina: Implications for conservation. Journal of Sea Research 83: 56-64 (2013}.
van derHeide, T., E. H. van Nes, G. W. Geerling, A. J. P. Smolders, T. J. Bouma, and M. M. van
Katwijk. Positive feedbacks in seagrass ecosystems: Implications for success in
conservation and restoration. Ecosystems 10: 1311-1322 (2007}
van der Heide.T., E. T. H. M Peeters, D. C. R Hermus., M. M van Katwijk., J. G. M Roelofs, and A. J.
P. Smolders. Predicting habitat suitability in temperate seagrass ecosystems. Limnology and
Oceanography 54: 2018-2024 (2009}.
van Katwijk, M.M., G. H. W. Schmitz, A. P. Gasseling, and P.H. van Avesaath. Effects of salinity
and nutrient load and their interaction on Zostera marina. Marine Ecology Progress Series
190:155-165(1999}.
van Katwijk, M.M., A.R Bos, V. N. de Jonge, L. S. A. M. Hanssen, D. C. R. Hermus, and D.J. de Jong.
Guidelines for seagrass restoration: Importance of habitat selection and donor population,
spreading of risks, and ecosystem engineering effects. Marine Pollution Bulletin 58: 179-
188(2009}.
Venables, W. N., and B.D. Ripley, B. D. Modern Applied Statistics with S. Fourth Edition. New
York: Springer, 2002.
Wazniak, C. E., and M. R. Hall (Eds}. Maryland's Coastal Bays: Ecosystem Health Assessment
2004, DNR-12-1202-0009. Annapolis, MD: Maryland Department of Natural Resources,
Tidewater Ecosystem Assessment, 2005.
Wetzel, R. L., and P. A. Penhale. Production ecology of seagrass communities in the lower
Chesapeake Bay. Marine Technology Society Journal 17: 22-31 (1983}.
Wicks, E. C., E. W. Koch., J. M. O'Neil, and K. Elliston. Effects of sediment organic content and
hydrodynamic conditions on the growth and distribution of Zostera marina. Marine Ecology
Progress Series 378: 71-80 (2009}.
Wortmann, J., J. W. Hearne, and J. B. Adams. A mathematical model of an estuarine seagrass.
Ecological Modeling 98: 137-149 (1997}.
Zuur, A. F., E. N. leno, N. J. Walker, A. A. Saveliev, and G. M. Smith. Mixed Effects Models and
Extensions in Ecology with R. New York, NY: Springer, 2009.
76
-------
R-package references
Andrus, J. M., T. Kuehlhorn, L. F. Rodriguez, A. D. Kent, and J. L. Zilles. CommunityCorrelogram:
Ecological Community Correlogram. R package version 1.0, 2014. http://CRAN.R-
project.org/package=CommunityCorrelogram.
Bjornstad, Ottar N. ncf: spatial nonparametric covariance functions. R package version 1.1-5,
2013. http://CRAN.R-project.org/package=ncf.
Brostrom, G. glmmML: Generalized linear models with clustering. R package version 1.0, 2013
http://CRAN.R-project.org/package=glmmML.
Fournier, D. A., H. J. Skaug, J. Ancheta, J. lanelli, A. Magnusson, M. Maunder, A. Nielsen, and J.
Sibert. AD Model Builder: Using automatic differentiation for statistical inference of highly
parameterized complex nonlinear models. Optimization Methods Software 27: 233-249
(2012}.
Fox, J. and J. Hong. Effect displays in R for multinomial and proportional-odds logit models:
Extensions to the effects package. Journal of Statistical Software 32: 1-24 (2009}.
Signorell, A. Includes R source code and/or documentation previously published by: Anderegg,
N., T. Aragon, A. Arppe, B. Bolker, S. Champely, D. Chessel, L. M. D. Chhay, H. C. Doran, S.
Dray, C. Dupont, J. Enos, C. Ekstrom, M. Elff, J. Fox, T. Galili, M. Gamer, J. Gross, G/
Grothendieck, F. E. Harrell, Jr., M. Hoehle, C. W. Hoffmann, M. Huerzeler, R. J. Hyndman, V. A.
Iglesias, D. Kahle, M. Kohl, M. Korpela, J. Lemon, M. Maechler, A. Magnusson, D. Malter, A.
Matei, D. Meyer, Y. Min, M. Naepflin, D. Ogle, S. Pavoine, R. Rapold, W. Revelle, V. E. Seshan,
G. Snow, M. Smithson, W. A. Stahel, Y. Tille, A. Trapletti, K. Ushey, J. VanDerWal, B. Venables,
J. Verzani, G. R. Warnes, S. Wellek, P. Wolf, and A. Zeileis. DescTools: Tools for descriptive
statistics. R package version 0.99.7, 2014. http://CRAN.R-project.org/package=DescTools.
Sing, T., 0. Sander, N. Beerenwinkel, and T. Lengauer. 2005. ROCR: Visualizing classifier
performance in R. Bioinformatics 21: 7881 (2005}.
Tierney, L., A. J. Rossini, N. Li, and H. Sevcikova. snow: Simple Network of Workstations. R
package version 0.3-13, 2013. http://CRAN.R-project.org/package=snow
Venables, W. N., and B. D. Ripley. Modern Applied Statistics with S. Fourth Edition. New York,
NY: Springer, 2002.
Wang, Xiao-Feng, and B. Wang. 2011. Deconvolution estimation in measurement error models:
The R Package decon. Journal of'Statistical Software 39: 1-24 (2011}.
Wickham, H. ggplot2: Elegant graphics for data analysis. New York, NY: Springer, 2009.
77
-------
Wood, S. N. Fast stable restricted maximum likelihood and marginal likelihood estimation of
semiparametric generalized linear models. Journal of the Royal Statistical Society (B) 73: 3-
36(2011}.
Wood, S. and F. Scheipl. gamm4: Generalized additive mixed models using mgcv and Ime4. R
package version 0.2-3, 2014. http://CRAN.R-project.org/package=gamm.
78
-------
Appendix A. National and Regional Data
Sources for Seagrass Habitat Model
Development
A- 1
-------
Table A-l. National and regional data sources for seagrass habitat model development
Parameter
Current
Seagrass
Data Source
Estuary Data Mapper
Sediments
Temperature
Rl DEM Narragansett
Bay Estuary Program
Estuary Data Mapper
Original Source
Rhode Island Geographic
Information System (RIGIS)/
Bradley, M., K. Raposa, and S.
Tuxbury. 2007. Report on the
Analysis of True Color Aerial
Photography to Map and
InventoryZostera marina L. in
Narragansett Bay and
Block Island, Rhode Island.
Page 1-16 and 9 Mapsheets.
Rhode Island Natural History
Survey, Kingston, Rl.
McMaster(1960)
Source
Date
2006
Salinity
Bathymetry
MURSST-Multi scale Ultra
high Resolution Sea Surface
Temperature (NASA - Jet
Propulsion Laboratory-
California Institute of
Technology)
R.I. Department of
Environmental Management.
Bay Assessment Response
Team.
Estuary Data Mapper NOAA Coastal Data Model
Estuary Data Mapper
1960
Description
Data layer of Submerged Aquatic Vegetation (SAV) in Rhode
Island Coastal Waters, 2006 from Rhode Island Geographic
Information System (RIGIS). SAV presence/absence coded
as: Present (1), Absent (0). Grid cell size 10m.
Link
http://www.edc.uri.edu/rigis/data/
Between 1988 and 1992 the Narragansett Bay Project
developed an extensive listing of CIS data layers for
applications involving Narragansett Bay. All data were
compiled at the University of Rhode Island's Environmental
Data Center in 1993 Prior to documenting, all data were
reviewed for spatial and topological correctness.
Used Estuary Data Mapper (EDM) to download daily sea
surface temperature (SST) from ftp://podaac-
ftp.jpl.nasa.gov/allData/ghrsst/data/L4/GLOB/JPL/MUR/,
and developed weekly mean SST data for the years 2003 to
2012. In areas of the estuaries which were not covered by
the SST data, a Euclidean distance method was used to fill in
areas without data.
Downloaded State of Rhode Island Department of
Environmental Management Bay Assessment & Response
Team (BART) daily salinity data and developed weekly mean
salinity data for stations in Narragansett Bay for the years
2003 to 2012. Thlessen polygons were created from the
point data to fill in all areas of Narragansett Bay.
Only one source of merged topographic and bathymetric
data is currently available in EDM, NOAA's Coastal Data
Model (Topography/Bathymetry (NOAA)
http://www.narrbay.org/biological_data.htm
http://mur.jpl.nasa.gov/index.php
http://www.narrbay.org/d_projects/
buoy/buoydata.htm
www.epa.gov/edm
A- 2
-------
Parameter
Wave Exposure
Data
Source
WeMO (Wave
Exposure Model)
Original Source
The wave exposure model is a
free software modeling tool
developed by NOAA to forecast
wave energy/exposure along
coastal areas. It provides
important habitat and
erosional information for
habitats. This software
application was developed by
NOAA but is no longer
supported.
Source
Date
2003-2006
meteorologic
al data.
Description
Data input for WeMO requires the following: a local
meteorological file (containing wind speed and direction),
bathymetry covering the area of interest and a shoreline file
which covers the boundary of the area of interest. A point
file (shape file format) must also be provided to allow the
model to perform calculations at specific points of interest.
Instructions on processing these data can be found in the
WeMO 4.0 manual available at the link provided. WeMo
model runs for this project were performed for each estuary
using default model values with a modified interrogation
distance of (5000m).
Link
http://www.csc.noaa.gov/
digitalcoast/tools/wemo
Transparency
Narragansett Bay
Commission; US
EPA's National
Coastal Assessment
program; MODIS
satellite imagery
Organic Carbon
USEPA National
Coastal Assessment
Narragansett Bay Commission;
US EPA's National Coastal
Assessment program; MODIS
satellite imagery
USEPA National Coastal
Assessment
Transparency across the bay was described using three data
sources. Secchi depths measured by the Narragansett Bay
Commission between 2008 and 2012
(http://snapshot.narrabay.eom/app/Monitoringlnitiatives/W
aterClarity) and by the US EPA's National Coastal Assessment
program between 2000 and 2006
(http://oaspub.epa.gov/coastal/coast.search) were combined
and averaged by Water Body ID (WBID) estuarine segments
used for assessment and listing by the Rhode Island DEM.
Secchi depths greater than the maximum depth at the site of
measurement were removed from the records before
averaging. Gaps in transparency data along the southern
shore of Conamicutt Island (Sakonnett and Newport Bays)
were filled in using Kd estimates from offshore MODIS
satellite imagery. The latter were downloaded using the EDM
tool, averaged over the growing season (May-October) for
nonmissing cells for the years 2008-2012. Average grid cell
values were extracted for cells greater than 30 meters in
depth, and averaged over a swath of offshore cells parallel to
the south shore of Conamicutt Island.
Total organic carbon was estimated across the Bay based on
surficial sediment grabs collected by the US EPA National
Coastal Assessment (downloaded from
http://oaspub.epa.gov/coastal/coast.search under Sediment
Grain Composition Data category). Values were interpolated
to create a complete grid within the shallow-water zone using
inverse distance weighting in ArcMap 10.1.
http://snapshot.narrabay.com/app/
Monitoringlnitiatives/WaterClarity;
http://oaspub.epa.gov/coastal/coast.search; EDM
(MODIS)
http://oaspub.epa.gov/coastal/coast.search
A- 3
-------
Data Original Source Source Description Link
Parameter Source Date
Secchi Secchi (3) data Secchi (3) data sources for varies Secchi depth was estimated from Kd values using an empirical http://www.dem.ri.gov/bart/netdata.htm;
sources for Narragansett Bay; Narragansett relationship developed by Batiuk et al. (2000). http://snapshot.narrabay.com/app/
Narragansett Bay; Bayfixed monitoring network
Narragansett Bay (RIDEM); Narragansett Bay Monitoringlnitiatives/WaterClarity; EDM
fixed monitoring Commission and MODIS data
network (RIDEM); Kd downloaded from EDM
Narragansett Bay
Commission and
MODIS data Kd
downloaded from
EDM
A- 4
-------
Appendix B. Tutorial on Finding and
Downloading Data for Seagrass Habitat
Prediction Models Using EPA's Estuary
Data Mapper
B.I. Overview of Estuary Data Mapper
The objective of EPA's Estuary Data Mapper (EDM) project is to produce an easily accessible, standalone
software product to automate the retrieval and pre-processing of GIS coverages (including remote
sensing data) and associated environmental data (e.g., tidal, hydrologic, and weather time series; water
quality and sediment quality data) to populate 1) a GIS data model for estuaries and their watersheds,
and 2) tools and models to assess, visualize, diagnose, predict, prioritize, and manage condition of
estuaries and coastal watersheds (Detenbeck et al. 2009). The EDM has been designed as a stand-alone
application requiring no other specialized software for implementation. EDM is written in C++, OpenGL,
and FLTK for extremely fast graphical visualization, user interaction and minimal memory consumption.
The interface has been designed with three tabs, the first enabling selection of an area of interest, with
background layers such as political boundaries, watersheds, estuaries, and hydrography provided for
reference and drop-down boxes to zoom in on states then estuary or watershed. The second tab allows
selection of geospatial layers and environmental data time series of interest to be selected for
download. This page also allows the user to visualize time series or other data before download. The
third tab allows the user to choose among a series of download formats (comma-delimited time series,
ASCII grid, shapefile, kmz for Google maps, png or mpg for visualization of images or time series) as well
as a location to save the files. Outputs also include a text file (edm_output.txt) containing all WCS calls
to rsigserver. This file demonstrates by example how to obtain data from the web service rather than
using the EDM GUI - so scripts and other applications can leverage it.
B- 1
-------
B.2. Installation of Estuary Data Mapper Tool
The latest version of EDM can be downloaded from www.epa.gov/edm. Right-click on the appropriate
line in the menu on the right-hand portion of the screen to download the zip file for the version of EDM
associated with your operating system and follow the instructions on this page to save as EDM.zip and
extract to a folder. Once EDM is unzipped, you can click on the EDM.bat file to start an interactive EDM
session.
apper | EMVL | US EPA
oFmpub,epa.gov
File Edit View Favorites Tools Help
Favotites
Federal Acquisition Institute... http--ordsecurity.epa.gov-r,.. © Research Management Syst, ,, ^Sponsors-Mentors Internsh... ^SHEMD Intranet - Safety -E.., >" http--rtpordh2,rtord.epa.g..,
I; ^Horne - ^r-
stfi... |©5tormw... jjj Snapsh.,, ;0 A turbid... ' g£http://s... |J|SignIn
_J Read Mai up Print - Page- Safety Tools '•*> He|p •
U.S. ENVIRONMENTAL PROTECTION AGENCY
High Performance Computing and Scientific Visualization
Search: OAII EPA ®This Area fS
I You are here: EPA Home » Programs «• Research » Envirormeftal ModeJna & VisuaBzation Lab fEMVL) » Estuary Data Mapper (EDM)
e-Estuary Project
Data Inventory
Estuary Data Mapper
I The Estuary Data Mapper (EDM) is a free, interactive graphical application under development by the US EPA that
| allows scientific researchers to quickly and easily:
:
* View and select estuary-scale geographical regions of interest - i.e., pan and zoom a map of the US coastal
states displaying watersheds, streams and estuaries, land coverage, satellite images, etc.
• Retrieve and view specified water-quality data, gauge measurements of tides and freshwater discharges,
gridded modeled nitrogen deposition, land use, topography/bathymetry, satellite data, etc., subsetted by
date-time and lon-lat rectangle. (Data streamed via web servers at EPA, NOAA, USGS, USFWS, etc. Examples
of water-quality data points include dissolved oxygen, turbidity and chlorophyll.)
* Save the resulting data subsets (as Shapefiles, grid files, etc.) for further analysis and import into the EPA
estuanne data model, ArcGIS or open source CIS, GoogleEarth, etc.
I The list of available maps and data is expanding during development and is shown in the EDM GUI. (Grayed-out
I buttons indicate a not-yet-irnplemented feature.)
EDM Mailing Lists
* EDM-Announcements. Receive periodic announcements of software updates and new data added to EDM.
Send an email to edmia'epa.gov with "EDM Announcements" in the subject line.
• EDM-Discussion. Join other EDM users in an email discussion group, Send an email to edm@epa.gov with
"EDM Discussion" in the subject line and give us permission to share your email address with other EDM
users.
nload Estuary Data
Mapper
Right-click on the
appropriate link below to
download the zip file for
your platform:
EDM Windows
EDM Mac OSX
EDM 32-bit Linux
EDM 64-bit Linux
EDM source code
Download Linked File
As... EDM.zip in your
home directory, (e.g., on
Windows this is
C:\Documents and
Settings\login\EDM.zip)
Click on EDM.zip file to
expand and it into a
folder called EDM.
Inside that folder, click
on the EDM script to
run it.
For help contact:
edm;gjepa,gov
919-541-5500®
| June 30, 2014 updated with new data:
• CMAQ CDC 2009 hourly 12km gridded deposition data CMAO CDC
• CMAQ/WDT 2001-2010 annual deposition data CMAO/WDT
* Kriged and Gridded Sediment data TNC GSM
• Coastal Vulnerability data USGS
I Screenshots:
• Multiple data layers: e.g. gridded deposition, water quality, tidal gauge stations, etc. retrieved with one click and displayed and saved in
seconds:
ft " Estuary Dala Mapper' Controller EOM'tf.ew 7SQx7SO [ 1.08929deg ^ 93.319km! ffi ( -74.9G7S3. 39.SQJ90)
'2Ge!Oata 3. Save Data Done
Figure B-l. Estuary Data Mapper home page.
B- 2
-------
B.3. Zooming to Area of Interest
The first time you start the EDM application, you will view the first tab and a default display extent:
1. Zoom Maps 2 Get Data 3. Save Data |4
• jmage MapQuest falt+i toggles to refresh)
r States Zoom To State
• Counties
• Cities
• Roads
Coastal Tributaries (main sterns) (NHD)
HUC8 (NHDPIUS) Zoom To HUC
EH Watersheds (e-Estuary Zoom To Watershed
Estuaries (NCA)
• Coastal Zone Management Areas
• Legislative Districts
: Reset Zoom
Reshape View: drag lower-right corner of window
Pan: drag left mouse button or use arrow keys
Zoom: + -,mousewheel,right-mouse,left-mouse+alt
vtr:i:n JU:3LMJJ (rtpmeta;
Figure B-2. First tab of Estuary Data Mapper with default display extent.
Although the default display extent consists of the Atlantic seaboard, EMD provides access to estuaries
and coastal watersheds across the entire conterminous United States. To facilitate finding an area of
interest, make sure the State, Watersheds (e-Estuary), and Estuanes (NCA) boundaries are turned on,
e.g., click the box next to States_until it turns white. When selections are toggled "on", the color of the
boxes next to the selections turns from black (not selected) to the color of the boundary line of interest.
To narrow down your search area, click on the Zoom to State button and select a state of interest.
B- 3
-------
1. Zoom Maps | 2. Get Data 3. Save Data 4.Done
maae MapQuest | (alt+i toggles to refresh)
P States
cSunties Alaska
Cities California
Roads Connecticut
Coastal Tributaries P.elaware
Florida
HUCB (NHDPIUS)
Georgia
watersneds (e-Esti wawai
Estuaries (NCA) Louisiana
Coastal Zone Mana Main-
Man/iand
Legislative Districts Massacnusetts
Mississippi
Reset Zoom I New Hampshire
Mew Jersey
New York
Pan: drag left rnou5eNotthCarQ||na
Zoom: +,-,rnousewnei Oregon
Pennsvlvania
South Carolina
Texas
Virginia
Washington
Version 2014070: ;rtprrn.\j!
i-or Help, edm@epa.gov yl9-541-5500
Figure B-3. Selection for zooming in to a state on first tab of Estuary Data Mapper.
Once zoomed into a state, when you click on Zoom to Estuary, you will be provided with a list of only
those estuarine boundaries appropriate for the selected state:
B- 4
-------
1_. Zoom Maps | 2. Get Data 3. Save Data 4_.Done
• Image MapQuest | (alt+i toggles to refresh)
States Zoom To State
Counties
Cities
Roads
Coastal Tributanes rrnam sterns) (TJHD)
HUCB (NHDPIUS) Zoom To HUC
• Watersheds (e-Estuary Zoom To watershed
r Estuaries (N
Coastal Zone Manaaement Area
Reset Zoom
Reshape View: drag lower-right corner of window
Pan: drag left mouse button or use arrow keys
Zoom: +,-,rnousewheel, right-mouse, left-mouse+alt
Version 2014070: ;rtpmi.\j!
Help, edm@epa.gov yl9-541-5500
Figure B-4. Process of zooming into an individual estuary in Estuary Data Mapper.
1. Zoom Maps | 2. Get Data 3. Save Data 4_.Done
image MapQuest (alt+i toggles to refresh)
P States Zoom To State
• Counties
• Cities
• Roads
Coastal Tnbiitai IPS imam sternal rNHD)
HUCS (NHDPIUS) Zoom TO HUC
• Watersheds re Estuary Zoom To Watershed
r Estuaries (NCA)
coastal Zone Manaaement Area:,
Leoislative Districts
Reset Zoom
Reshape View: drag lower-right corner of window
Pan: drag left mouse button or use arrow keys
Zoom: +,-,mousewheel,right-mouse,left-mouse+alt
For Help: edrn@epa gov 913-54 1-5500
Figure B-5. Estuary Data Mapper zoomed in to Narragansett Bay.
B- 5
-------
The display extent will shift to the boundaries of the selected estuary. The view can be reshaped using
your mouse following the instructions on the left-hand panel.
Alternatively, if you wish to identify and download data associated with the full watershed for an
estuary, you would click on Zoom to Watersheds (e-Estuary) to zoom to the associated watershed
boundary.
1. Zoom Maps | 2. Get Data 3. Save Data |4.Done
• Image MapQuest | (alt+i toggles to refresh)
P States Zoom To State.
• Counties
• Cities
• Roads
Coastal Tnoutarres (main stems} (NHD)
HUC8 (IMHDPIUS) Zoom To HUC..
(e-Estua
^= Rl Buzzards Bay
F Estuaries (IMCA) _ R| Litt|e Ngrrgggn5ett
• Coastal Zone Management R] Long Island Sound
Reset Zoom
Reshape View: drag lower-right corner of window
Pan: drag left mouse button or use arrow keys
Zoom: +,-,mousewneel, right-mouse Jeft-mouse+alt
Version 2ri407rr tipinc:;)
: edmasepa gov ':'ii9-541-5500
Figure B-6. Zooming in to an estuarine watershed in Estuary Data Mapper.
B- 6
-------
Image MapQuest (alt+i toggles to refresh)
States Zoom To State
Counties
Cities
Roads
Coastal Tribi.itaiies (main sterns) (TJHD)
HUCB (NHDPIUS) ZoomToHUC
Watersheds (e-Estuary Zoom To watershed
Estuaries (NCA)
Coastal Zone Manaaernent Area
Legislative Districts
PP'5Pt ^Ll
Reshape View: drag lower-right corner of window
Pan: drag left mouse Putton or use arrow keys
Zoom: +,-,mousewheel, right-mouse, left-mouse+alt
Help, edm@epa.gov yl 9-541-5500
Figure B-7.Get Data tab of Estuary Data Mapper.
Now you will see the full extent of the watershed of interest. To begin exploring data of interest, you
now click on the 2 Get Data tab on the top of the control panel.
B.4. Identifying Data Sources
B.4.1. System Boundaries
System boundaries that are checked on Tabl are automatically downloaded when you save other data
from that system. However, normally a smoothed version of the watershed boundaries is saved to
speed up the process of saving data. If the original higher-resolution watershed boundary is desired,
click on the Save Map Polygons/DBF box in Tab 3^ Save Data before selecting Save Data (Shapefile) and
type in an appropriate directory location to save the shapefiles to (Choose Directory Folder for Saved
Files):
B- 7
-------
Save Map Polygons/DBF (VERY SLOW;
Files1
3
. shp
.shx
.xml
36 Qvei:all_BQunds.pi:j
36 Overall_Bounds.3hp
36 Overall Bounds.shx
Figure B-8. Save Data tab on Estuary Data Mapper.
B- 8
-------
B.4.2. Seagrass
Click on the second tab to start the data discovery process. There is currently only one option for
retrieving seagrass coverages, which returns a composite coverage with data combined across multiple
sources representing the most up-to-date sources publically available as of 2013. Coding of seagrass
abundance categories was standardized as explained in the metadata and original sources are also
described in the metadata. Click on the Seagrass (State/NOAA) button and then on the Retrieve and
Show Selected Data button. In the text box at the base of the menu panel you will messages indicating
which data are being retrieved, when the data retrieval has finished, and how many records were
returned. In this case, because there is only one composite coverage, any dates you enter will be
ignored. In this example, we have manually zoomed into an area in the southern portion of
Narragansett Bay in order to see the seagrass pixels. (If you do this, you will need to remember to zoom
back out to the full extent before saving or you will only retrieve this subset.)
1. ^^^^^ 2. Get Data 3. Save Data 4.Done
Choose Scenario: Custom
Tidal Gauge Stations (NOAA)
Water Quality STORE!
Buoys
HDM Stride p TBQFS lay(l
Sediment USGS SEABED Calculated
Soil (STATSGO)
Coastal Vulnerability
Wetlands (USFWS)
Seagrass (State/NOAA)
Precipitation (PRiSM 1971-2000)
Temperature (PRISM 1971-2000 Monthly Higr
Topography/Bathymetry (NOAA)
Land Use CCAP 2006
Deposition
Date(Y/M/D)|2001 /Es/pl Days I 184
be ectec Data Timeout yui
Retrieving b'eagia^ CState/NQAA).,.
Finished retrieving data.
Retrieved 1137422 points of Seagrass (State/NC
ll |
Figure B-9. Retrieving seagrass data in Estuary Data Mapper.
B- 9
-------
B.4.3. Depth
Only one source of merged topographic and bathymetric data is currently available in EDM, NOAA's
Coastal Data Model (Topography/Bathyrnetry (NOAA)). Select this button and hit Retrieve and Show
Selected Data again. Note that if you hover your mouse over the item in the menu, you will see a popup
screen displaying the url for information on the original data source. Sources for data are also available
on the EPA EDM web page at http://ofmpub.epa.gov/rsig/rsigserver?edm/data_inventory.html.
1. Zoom Maps | 2. Get Data 3. Save Data 4.Done
Choose Scenario: Custom
Tidal Gauge Stations (NOAA) M
Water Quality STORE! | Turbidity
Buoys
IOOS
Water Temperature
Current
USGS SEABED Calculated
Org. Carbon (
CVI (-)
F Snow Labels
Soil (STATSGO)
Coastal Vulnerability
Wetlands (USFWS)
Seagrass (State/NOAA)
Precipitation (PRISM 1971-2000) Annual
Temperature (PRISM 1971-2000 Monthly High
Topography/Bathymetry (NOAA)
Land Use CCAP 2006
Deposition CMAQ | Hourly NOx
Satellite Depth tTI MODIS Attenuation
Pop. Density (ICLUS) Population A1 2010
Date(Y/M/D) |2001 /|o5/|oT Days | 184
Retrieve £ snow Selected Data J Timeout (900"
Play Delay 100
| Tirnestep
Retrieving I opoqraphy/Cdthynetry (NuAA)...
Finished retrieving data.
Retrieved 703964 points of Topography/Bathymi
jU I "±J
For Help: edm@epa.gov 919-541-5500
Figure B-10. Retrieving topobathymetry data in Estuary Data Mapper.
Users requiring finer resolution bathymetry data may wish to check local sources or the NOAA Tsunami
Inundation Digital Elevation Models (OEMs; http://www.ngdc.noaa.gov/mgg/inundation/tsunami/),
available for select regions only.
B- 10
-------
B.4.4. Transparency
Multiple data sources are available to describe transparency, including light attenuation coefficients,
Secchi depth, turbidity, chlorophyll, and total suspended solids. Data associated with grab samples or
instantaneous sensor readings for water quality parameters can be retrieved through web services
provided by the National Estuarine Research Reserve System (NERRS) program or by the joint USGS/EPA
(STORE! button) web services. Click on the STORE! button next to Water Quality to see the drop down
list which allows you to select between the NERRS and STORET web services. Clicking on the button to
the right of STORET activates the dropdown list to allow the user to choose a water quality parameter.
At this time, only one parameter can be selected at a time, but users can choose to download results
then select another parameter. Retrieval of water quality values requires that the user select a starting
date and number of days of record (up to 365) for retrieval. Depending on the time of day, available
bandwidth, demand for web services, and amount of data retrieved, these requests can take up to a few
minutes. The user may need to increase the timeout parameter to the right of the Retrieve & Show
Selected Data button to up to 900 seconds to prevent the tool from timing out and returning a message
indicating no data points are available.
1. Zoom Maps 2 Get Data 3. Save Data 4.Dane
Choose Scenario: Custom
Soil (STATSGO)
Coastal Vulnerability
Wetlands (USFWS)
Seagrass (State/NOAA)
Precipitation (PRISM 1971-2000) Annual
Temperature (PRISM 1971-2000 Monthly Hign
Topography/Bathymetry (NOAA)
Land Use CCAP 2006
Deposition
CMAQ
Hourly NOx
_.
Satellite Depth 0
MODIS Attenuation
Pop. Density (ICLUS) Population A1 2010
iate(YAVD) |2005 /|o?/|oT Days 1184
Retrieve & Show Selei
Tirnestep
Play Delay 100
Retrieving Water Quality...
Finished retrieving data.
Retrieved 0 points of Water Quality
For Help edtTn.^epa guv 9 19-54 1-5500
1. Zoom Maps 2. Get Data 3. Save Data 4.Done
Choose Scenario: Custom
Tidal Gauge Stations (NOAA)
F Water Quality _
^- Station
Discharge
HDM Stride 11 TBQFS \\3] Crilorophyll-a
Sediment USGS SE/ Chlorophyll-a Non-Pt
r Show Labels Org.
Soil (STATSGO)
Coastal Vulnerability
Wetlands (USFWS)
Chlorophyll-b
Chlorophyll-c
Color Apparent
Color True
CDQM
Conductivity
PH
Seagrass (State/NOAA)
Precipitation (PRISM 1971-2 sa|jnjty PSIJ
Temperature (PRISM 1971-; Salinity PPM
Topography/Bathymetry (NC Dissolved Fluoride
Dissolved Iron
Land Use CC
1 Dissolved Nitrogen
Deposition CMAQ | Dissolved Oxygen
Satellite Depth 0 | MQ Dissolved Oxygen %
Pop. Density (ICLUS) P NH3/NH4
— Nitrite
Date(YMTJ) |2005 / |ui / |oT [ Nitrate
N02+N03
Retrieve & Show Selected Dat Organic Carl
Organic Nitrogen
Pla^ Organic Phosphorus
Pheophytin
Phosphate
Retrieving Water Quality... Silicate
Finished retrieving data. Temperature
Retrieved 0 points of Water QL Total Nitrogen
Total Phosphorus
For Help: edrn@epa gov 919- Secchi Disk DePtn
"don
| J
Timestep
Figure B-ll. Selecting STORET variables for data retrieval in Estuary Data Mapper.
Queries through STORET web services often fail to locate data collected through EPA's National Coastal
Assessment surveys (e.g., Secchi depth, light attenuation coefficients, TSS, chlorophyll) so users may
wish to retrieve data directly from the EPA's Environmental Monitoring and Assessment web site
(http://oaspub.epa.gov/coastal/coast.search ) for data up to 2006 or the National Coastal Assessment
site (http://water.epa.gov/type/oceb/assessmonitor/ncca.cfm ) for data from 2008 or later (due to be
B- 11
-------
added by the end of 2014). Selected parameters from EPA's NCA surveys will be added to EDM in the
near future to fill the gap in WQ web services.
Remotely sensed light attenuation coefficients based on the MODIS satellite are also available through
NASA web services. Caution should be exercised in using these data for shallow systems (< 30m depth)
as algorithms were developed for the open ocean and do not include corrections for bottom reflectance.
Therefore, data should be checked against in situ measurements to determine if a relationship exists.
Algorithms for ocean color parameters are currently under development and are being tested for
shallow and more turbid coastal systems (Keith et al. 2014) but these data are not yet readily available
through web services.
• 1 Estuary Data Mapper; Controller
1 Zoom Maps
Choose Scenario'
3 Save Data 4-Done
Tidal Gauge Stations (NOAA)
Water Quality NERRS
Buoys NDE
Water Temperature
CBOFS |lay[i Water Temp
I Show Labels
Sand(%)
Soil (STATSGO) Organic (%wt)
Coastal Vulnerability CVI <-)
Wetlands (USFWS)
Seagrass (State/NOAA)
Precipitation (PRISM 1971-2000} Monthly
Temperature (PRISM 1971-2000 Monthly High
Topography/Bathymetry (NOAA)
Land Use CCAP 2006
Deposition NAOP | TotN
r Satellite Depth 0
MODIS Attenuation
I Pop Density (ICLUS) Population A1 2010
Date(Y/M/D)|2010/|o?/|o7 Days] 185
Retneve & Show Selected Data TimeoutI3GOO
Play Delay 100
Timestep
Finished retrieving data.
Retrieved 644809 points of Topography''Bath;,
Retrieved 4915265 points of Satellite
j
-------
B.4.5. Energy Environment
B.4.5.1. Wave Energy Model Inputs
There are a few options available for the user to calculate Relative Wave Energy. Previously, users could
use the NOAAWave Energy Model (WEMO; http://products.coastalscience.noaa.gov/wemo/) to
calculate fetch and relative wave energy based on input data including system boundaries, merged
topo-bathymetry, and wind data. Unfortunately WEMO is no longer being supported and is not
compatible with ArcMap versions later than 9.3. Recently, the USGS WAVE extension for ArcMap was
upgraded for use with ArcMap 10 and could be used in a similar fashion. See Sections B.4.1 and B.4.3
for information on downloading data on system boundaries and merged topo-bathymetry. Wind data
can be downloaded from the NERRS web service (wind speed and direction or from various NOAA web
services for buoy data (jiuoys button, using IOOS, NDBC or NERACOOS dropdown options).
B.4.5.2. Relative Exposure
USGS has recently calculated a Coastal Vulnerability Index, including a component related to wave
energy (http://woodshole.er.usgs.gov/project-pages/cvi/). Unfortunately values have been collapsed
onto an ordinal scale (1-5), leading to a loss of information. For some systems such as the Narragansett
Bay (below), the range of values may be restricted, making this index less useful as a predictive tool.
However, for other systems with a wider range of energy environments, the Wave Rank score might
prove useful.
1. Zoom Maps 2. Get Data 3. Save Data 4.D
Choose Scenario: Custom
Tidal Gauge Stations NOAA
Water Quality NERRS
Buoys NERACOOS
HDM Stride!"
Sediment
US'JS SEABED i acua en
Soil (STATS6O)
Coastal Vulnerability
Wetlands (USFWS)
Wave Rank [1-5]
Length (km)
CVIf-)
Seagrass (State/NOAA cv, RankM.4
Precipitation (PRISM 1! CVI Risk
Temperature (PRISM 1 Mean Wave Height j
Topography/Bathymetr
Wave Risk
Lana use Mean Tide Range (m
Deposition CMAQI Tide Rank [1-5]
ellite Depth 0 ! Tlde Risli
Pop. Density (ICLUS)
Sea Level Rise (mrn/yr;
Sea Level Rank [1-5]
Sea Level Risk
Date(Y/M/D) 12010 / US' /
: Slope (%)
Retrieve £ Show Selects 3l°Pe Rank I1'5]
Slope Risk
Erosion/Accretion (m/yr)
Timestep ^ | Erosion Rank [1-5]
Erosion Risk
Geomorphology Rank 1-
Retnevmg Coastal Vulne Geomorph0|0gy Risk
Finished retrieving data. —
Retrieved 252 points of Coastal Vulnerability
For Help: edrn:;:'e la :L> ^'i ^-541-5500
Figure B-13. Retrieving Coastal Vulnerability Index data in Estuary Data Mapper.
B- 13
-------
B.4.5.3 Current velocity
Current velocity data are available at some Buoy locations; however, these stations tend to be sparse.
For a limited number of systems, a more complete coverage of the current velocity environment can be
obtained from model output. JHydrorneteorological model outputs for NOAA's Operational Forecast
Systems are available for seven systems. Choose system using dropdown menu by HDM button for the
layer of interest.
1. Zoom Maps | 2. Get Data 3. Save Data 4.Done
Choose Scenario: Custom
Tidal Gauge Stations (NOAA)
Water Quality NERRS | Wind Speed
Buoys IOOS
HDM Stride
Sediment
Soil (STATSGO)
Coastal vulnerability
Wetlands (USFWS)
Seagrass (State/NOAA)
Precipitation (PRISM 1971-2000)
Temperature (PRISM 1971-2000 Monthly High
Topography/Bathymetry (NOAA)
Land Use CCAP 2006
Deposition
Date(Y/M/D)|2012/|09/|l2 Daysp
Retrieve £ Show Selected Data :l Timeout I'?OC
Retrieving HDM...
Finished retrieving data.
Retrieved 935664 points of HDM
Figure B-14. Retrieving current velocity data in Estuary Data Mapper.
B.4.6. Sediment Characteristics
Sediment characteristics are available both for grab samples from US EPA National Coastal Assessments
(URL) or the USGS Seabed database (URL) and as continuous grids developed using kriging methods for
the North Atlantic coast
(https://www.conservationgateway.org/ConservationByGeography/NorthAmerica/UnitedStates/edc/re
portsdata/marine/namera/Pages/default.aspx) or Gulf coast
(http://instaar.colorado.edu/~jenkinsc/dbseabed/resources/gsmseabed/). You can choose between
these options on the Sediment selection (a), which leads you to the drop down menu for sediment
characteristics for the NCA (b), or USGS Seabed datasets (c), or for the krigged datasets (d) and e):
B- 14
-------
a)
1. Zoom Maps [ 2. Get Data 3. Save Data 4.Done
Choose Scenario: Custom
I Tidal Gauge Stations (NOAA) MS
I Water Quality NERRS | Wind Speed
I Buoys
r Show Labels USGS SEABED Calculated
I Soil (STATSGO) USGS SEABED Extracted
I coastal vulnerability USGS SEABED Parsed
Wetlands (USFWS)
TNC Kriged
GSM Gfldded
• Seagrass (State/NOAA"]
• Precipitation (PRISM 1971-2000) Annual
• Temperature (PRISM 1971-2000 Monthly High
I Topography/Bathymetry (NOAA)
• Land Use CCAP 2006
• Deposition CMAQ Hourly NOx
r Satellite Depth 0~~| MODIS Attenuation
• Pop. Density (ICLUS) Population A1 2010
Date(Y/M/D)]2010 /|o5"/f)T Days| 184
Retrieve & Show Selected Data Timeout J90Q
Timestep
Play Delay in
Version 20140702 (rtpmeta)
For Help: edrn@epa.gov 919-541-5500
1. Zoom Maps | 2. Get Data 3 Save Data 4.Done
Choose Scenario: Custom
I Tidal Gauge Stations (NOAA) MS
I Water Quality NERRS | Wind Speed
I Buoys IOOS
Current Speed
Sand
I HDM Stride 11 CBOFS |lay[l Current
I Sediment GSM Gridded
P Show Label:!
I Soil (STATSGO)
I Coastal Vulnerability
I Wetlands (USFWS)
HYG
Mudu
Wave Rant Sand
Gravel
Rock
• Seagrass (State/NOAA)
• Precipitation (PRISM 1971-2000) Annual
• Temperature (PRISM 1971-2000 Monthly High
F Topography/Bathymetry (NOAA)
• Land Use CCAP 2006
• Deposition
CMAQ
Hourly NOx
r Satellite Depth 0
MODIS Attenuation
I Pop. Density (ICLUS) Population A1 2010
Date(Y/M/D)|2010/|rj5"/(bT Days| 164
Retrieve S Show Selected Data Timeout J900
Play | Delay(lO
4 I Timestep ^ |
I
ersion 20140702 (rtpmeta)
d) For Help: edrn@epa.gov 919-541-5500
b)
1. Zoom Maps 2 Get Data | 3. Save Data |4.Done
Choose Scenario: Custom
• Tidal Gauge Stations (NOAA) MSL
• Water Quality NERRS Wind Speed
• Buoys IOOS Current Speed
I Soil (STATSGO)
I Coastat Vulnerability
I Wetlands (USFWS)
Wave Longitude (deg)
Latitude (deg)
• Seagrass (State/NOAA) c!ay *%j
• Precipitation (PRISM 1971-200C SiltClay (%)
• lemperature (PRISM 1971-200 silt (%)
I Topography/Bathymetry (NOAA
• Land Use CCAP 25tn% Pni (p
• Deposition CMAQ
P Satellite Depth 0
• Pop. Density (JCLUS)
Date(Y/M/D) (201 o /|rJ5"/[oT nay Agency
Hi 50th% Phi (phi)
•J^Ji 75th% Phi (phi)
—^37 Deviation (-)
Skewness (-)
Retrieve & Show Selected Data Timeout |goo
J
Timestep
Play Delay 10
Version 20140702 (rtpmeta)
For Help: edm@epa.gov 919-541-5500
1. Zoom Maps 2. Get Data 3. Save Data 4.Done
Choose Scenario: Custom
I Tidal Gauge Stations (NOAA) MS
I Water Quality NERRS [ Wind Speed
I Buoys iOOS
Current Speed
• HDM Stride 1 CBOFS |iay|l Current
• Sediment TNC Kriged
P Show Label:!
• Soil (STATSGO)
• Coastal Vulnerability Wave Rank [1-5]
• Wetlands (USFWS)
• Seagrass (State/NOAA)
• Precipitation (PRISM 1971-2000) Annual
• Temperature (PRISM 1971-2000 Monthly High
r Topography/Bathymetry (NOAA)
• Land Use CCAP 2006
• Deposition CMAQ Hourly NOx
r Satellite Depth 0
MODIS Attenuation
I Pop. Density (ICLUS) Population A1 2010
Date(Y/M/D)|2010 /|os'/|oT Days| 184
:::. bhuvv 'ntTle.:i"i.i DJM lineout yfiu
| Play | Delayprj
4 | Timestep ^ |
Version 20140702 (rtpmeta)
For Help edmttjep.a gov 919-541-5500
c)
1. Zoom Maps | 2. Get Data 3. Save Data ]4.Done
Choose Scenario: Custom
I Tidal Gauge Stations (NOAA) MS
I Water Quality NERRS | Wind Speed
I Buoys IOOS
Current Speed
Current
I Soil (STATSGO)
I Coastal Vulnerability
I Wetlands (USFWS)
_ Sample Key
I Seagrass (State/NOAA)
Wavf Longitude (deg)
Latitude (deg)
Water Depth (m)
Sample Top (m)
• Precipitation (PRISM 1971-2001 Samp|e Base (m)
• Temperature (PRISM 1971-200 Carbonate (%)
r Topography/Bathymetry (NOAA 0rS- Carbon (%)
• Land Use
• Deposition CMAQ
r Satellite Depth 0 I
Hi
- Mud(%)
Sand (%)
MODIS Gravel (%)
• Pop. Density (ICLUS) Popi Porosity (%)
Shear Strength (Ic
Date(Y/M/D) |2010 / ID? / |o7 Day Crit. Shear Stress
Retrieve £ Show .^e ECted r
_n
_4 | Timestep
Version 20140702 (rtpmeta)
Sorting (-)
. Folk Code
P|3V I Shepard Code
Site Name
Sample Phase
Sampler
For Help: edm@epa.gov919-541-5500
Figure B-15. Drop-down selections for selection of sediment parameters from EPA National Coastal Assessment dataset.
B- 15
-------
B.4.7. Water Quality
B.4.7.1. Temperature
EDM provides access to instantaneous temperature readings through the EPA/USGS web services
(Section 4.4) time series at fixed stations through the NERRS or various buoy web services (Section
4.5.1), and remotely sensed time series with continuous gridded coverages. Although several sources
are available in the Satellite drop-down menu, the MUR option provides grids with the finest resolution
and best coverage for many estuarine systems. Note that the color scale has been optimized to
represent the range in temperature over time, so will not do a good job of illustrating gradients in
temperature for a particular day within EDM. Users investigating the effects of temperature for
systems in which seagrasses are exposed to the air can also access local air temperatures through the
web services for NERRS or NOAA's buoy systems.
HDM Stride|1 CBOFS |lay[l Current
Sediment TNC Kriped
• Soil (STATSGO)
• Coastal Vulnerability
• Wetlands (USFWS)
• Seagrass (State/NOAA)
• Precipitation (PRISM 1971-2000) Annual
• Temperature (PRISM 1971-2000 Monthly High
• Topography/Bathymetry (NOAA)
• Land Use CCAP 2006
• Deposition CMAQ
F Satellite Depth 0
• Pop. Density (ICLUS)
Retrieve £ snow Selected MODIS Da^me SST
MODIS Nighttime SST
F MODIS Attenuation
Timestep > ,~ MODIS Fluorescence
MODIS Chlorophyll
MODIS Monthly CDOM
Retrieving Satellite ..
Finished retrieving data. SEAWlFS cnigrophyll
Retrieved 863328 points of
WOA Monthly Water Tei
WOA Monthly Salinity
For Help: edm@epa.govE WOA Monthly Oxygen
Figure B-16. Retrieval of remotely sensed temperature data with Estuary Data Mapper.
B.4.7.2. Salinity
EDM provides access to instantaneous salinity readings through the EPA/USGS web services (Section
4.4), time series at fixed stations through the NERRS, or various buoy web services (Section 4.4.1), and
remotely sensed time series with continuous gridded coverages. The latter are available through NASA
web services with data aggregated from daily to annual time steps. Different satellite coverages can be
selected from the dropdown box selection under Satellites (see below).
B- 16
-------
1. Zoom Maps 2. Get Data 3 Save Data 4 Done
Choose Scenario: Custom
Tidal Gauge Stations (NOAA) MS
Water Quality NERRS | Wind Speed
Buoys IOOS
Current Speed
HDM Stride [l~ CBOFS [lay [T~ Current
Sediment
TNC Kriged
Show Labels
Sediment
Soil (STATSGO) HYGRP
Coastal Vulnerability Wave Rank [1-5]
Wetlands (USFWS)
Seagrass (State/NOAA)
Precipitation (PRISM 1971-2QOQ) Annual
Temperature (PRISM 1971-2000 Monthly High
Topography/Bathymetry (NOAA)
Land Use CCAP 2Q06
Deposition CMAQ Hourly NOx
r SateBite Depth 0
• Pop. Density (ICLUS) Po Aquarius Weekly SE
i i— i— Aquarius Monthly S:
Retrieve S Show Selected Date MODIS Da^ime SS"
— MODIS Nighttime SE
MODIS Attenuation
MODIS Fluorescenc
MODIS Chlorophyll
— MODIS Monthly CDi
Retrieving Satellite.. MUR Water Temper
Finished retrieving data SEAWIFS Chloroph-
WOA Monthly Wate
WOA Monthly Salini
1
Timestep
Play
Retrieved 863328 points of Sat*
V
For Help: edm@epa.gov 919-! WOA Montn|y Osyg,
Figure B-17. Selection of remotely sensed salinity data for download in Estuary Data Mapper.
B.4.8. Nitrogen Concentration and Loading
B.4.8.1. Nutrient Concentrations
Nutrient concentrations for estuaries and their tributaries can be retrieved via the EPA/USGS and NERRS
web services (see Section B.4.4). As mentioned previously, neither the EPA/USGS web service nor the
online STORE! database query interface reliably allow the retrieval of data collected during NCA surveys.
See Section B.4.3 for directions on current web access to NCA data sets.
B.4.8.2. Atmospheric Loads
Data regarding nitrogen sources or loads to estuaries and their watersheds are available through the
Nitrogen menu in EDM. Atmospheric loading data for nitrogen and phosphorus can be retrieved
through EDM based on either 1) interpolation of monitoring data collected by the National Atmospheric
Deposition Network (NADP; http://nadp.sws.uiuc.edu/) or 2) modeled deposition based on results of
modeling runs of EPA's Community Multiscale Air Quality model
(http://www.epa.gov/AMD/Data/wdtData.html; http://www.epa.gov/heasd/research/cdc.html). CMAQ
deposition data are available to estimate both deposition over the estuarine watershed as well as
directly to the estuary surface area for hourly, monthly or annual time steps. NADP grids only cover
watershed deposition for an annual time step. Summaries of loading by estuary are covered below in
Section B.4.8.3.
B- 17
-------
' EsTuiiiy DaU Mapptr Ccntrolkr
1 Zoom Maps | 2 Get Dala 3. Save Data ';4.00fM
PICK Scenario: Custom
-us
aU.Dooell
STORET
MLW
Seccni DISK Deplfi
ND8C
Wind
arid*|i
Sediment
EPA MCA
roc (°i)
Organ*
CVl(-)
I SMI (STATSGO)
I Coastal yuinerawiiy
I Wetlands (USFWS)
I Seagrass (Slale/NOAA)
I CHmale (PRISM '71-'OQ) Pretip. Monthly Mean |
I Topograpnyraainymetry (NOAA)
I Land Use CCAP 2006
' Nitrogen HADP Gria
I satellite Pepth o | MODIS Daytime SST
I Population C1CIUS) Ai Yearty Pop Den
rears|~T
Retrieve & Snow Selected Data ( Timeout J300
Piaylelay
Timestep _LJstride
Show Polnl Data Labels
[Version 20150710 (rtpmeta)
For Help: edm@epa gov 919-541-5500
" Nitrogen NADP Grjd
4
Satellite CMAQ Grid
Pop. Dens Source: Gridde
-------
•J Estuary Data Mapper Centre I
1 Zoom Maps 2. Get Dat;
Choose Scenario
ar I c= | E] '[ S3
3. Save Data 4.Done
Custom
• Tidal Gauge Stations (NOAA) MSL
• Water Quality NERRS
• Buoys NDBC
Turbidity
Water Temperature
• HDM Stride 1 1 CBOFs|lay1 Water Temp
• Sediment
• Show Labels
• Soil (STATSGO)
• Coastal Vulnerability
EPA NCA
Sand (%)
Organic (%wt)
CVl (-)
• Wetlands (USFWS)
• Seagrass (State/NOAA)
• Precipitation (PRISM 1971-2000) Monthly I
• Temperature (PRISM 1971-2000 Monthly High |
• Topography/Bathymetry {NOAA)
• Land Use CCAP 2006
' T Nstrogei CMAQGrd
• Satellite Depth o
• Pop Density QCLUS)
Monthlj N
MODIS Attenuation
Population A1 2010
Date(Y/M/D) |2006 3 ''|°~i~^ '|^"^j Days | 36o
Retrieve & Show Selected Data Timeout J3600
4 | Timestep fr- |
Play | Delay |lOO
Version 20140926 (rtpmeta)
For Help: edm@epa.gov 919-541-5500
750 [ 0.70000deg - 58.161km] 9
Figure B-19. Retrieval of CMAQ monthly nitrogen deposition data with Estuary Data Mapper.
| '*_"'• Estuary Data Mapper Controller
1 Zoom Maps | 2. Get Data 3
Choose Scenario
• Tidal Gauge Stations (NOAA}
• Water Quality NERRS
• Buoys
NDBC
Date(Y/M/D) boOS J / B7 J / BT3 Days [365
• HDM Stndi
• Sediment
• Show Labels
• Soil (STATSGO)
• Coastal Vulnerability
• Wetlands (USFWS)
• Seagrass (State/NOAA)
• Precipitation (PRISM 1971-2000) Monthly
• Temperature (PRISM 1971-2000 Monthly High
• Topography/Bathymetry (NOAA)
• Land Use CCAP 2006
r Nitroge
Figure B-20. Retrieval of NADP nitrogen deposition data with Estuary Data Mapper.
B- 19
-------
B.4.8.3. Watershed Sources
Nitrogen source and loading data are available through EDM in gridded form for estuarine watersheds.
In addition, summaries of annual watershed-based loads and sources and direct atmospheric loading to
estuaries by estuary are provided in dbf files associated with estuary and watershed shapefiles.
Atmospheric deposition data can be accessed from the Nitrogen menu (a), which gives access to both
gridded data at hourly (CMAQ), monthly (CMAQ), or annual time steps (CMAQ, NADP). In addition,
summaries of annual CMAQ or NADP atmospheric N deposition and monthly CMAQ atmospheric N
deposition are available for estuarine watersheds and estuaries/subestuaries (CMAQ only). NADP
deposition grids have been interpolated for terrestrial areas only.
Selection of the CMAQ grid menu (b) provides access to a wide array of time steps and nitrogen forms.
a)
r Nitrogen j NADP Grjd
• Satellite CMAQ Grid
• Pop. Dens Source: Gridded
CMAQ Estuary
j CMAQ Subestuary
Date(Y/M/D)
Src:Watershed NADP
Src:Watershed CMAQ
Retrieve &; src:Watershed Point
• Src:Watershed NonP
Load: Estuary Spa92
Ld:Est Spa02 MRB1
Ld:Est Spa02 MRB2
Time
Version 201<
Ld:Est Spa02 MRB7
Ld:Subestuary Spa92
For Help: e Ld:Sub sPa02 MRB1
Ld:Sub Spa02 MRB2
Ld:Sub Spa02 MRB5
LdrSub Spa02 MRB7
Hourly Dry NOx
Hourly Wet NOx
Hourly NOx
Monthly Dry Ox. N
Monthly Dry Red. N
Montnly Dry N
Monthly Wet Ox. N
Monthly Wet Red. N
Monthly Wet N
Monthly Ox. N
Monthly Red. N
Month lyN
Monthly Dry Sulfur
Monthly Wet Sulfur
Monthly Sulfur
Yearly Dry Ox. N
Yearly Dry Red. N
Yearly Dry N
Yearly Wet Ox. N
Yearly Wet Red. N
Yearly Wet N
Yearly Ox. N
Yearly Red. N
Yearly N
Yearly Dry Sulfur
Yearly Wet Sulfur
Yearly Sulfur
Yearly Dry Hg
Yearly Wet Hg
Yearly Hg
Yearly AWD OXN T
Yearty AWD REDN T
Yearly AWD S T
Yearly AWD SS S
YearlyAWDEPCL
Yearly AWDEP NA
Yearly DD OXN T
Yearly DD REDN T
Yearly DD S T
Yearly DD SS S
Yearly DDEP CL
Yearly DDEP NA
Yearly TD OXN T
Yea rlyTD REDN T
Yearly TD S T
Yearly TD SS S
Yearly TDEP CL
Yearly TDEP NA
Yearly TD N
c)
d)
e)
Estuary Code
Monthly Dry N
Monthly Wet N
Monthly N
Yearly Dry N
Yearly Wet N
Yearly N
Estuary Code
Total Load
Wet Load
Total Yield
Wet Yield
Figure B-21. Nitrogen submenus for CMAQ and NADP nitrogen deposition.
B- 20
-------
Selection of NADP grid sources in the main menu then provides access to a submenu (c) allowing access
to different N fractions: total ammonia, nitrate and total N. Selection of CMAQ summaries at the estuary
(d) or watershed (e) scales provides access to choices between wet, dry or total N fractions aggregated
at monthly or annual time steps.
Additional land-based nitrogen sources (a) can be assessed with b) gridded maps or c) at the estuarine
watershed scale.
a)
F Nitrogen
• Satellite
• Pop. Dens
NADP Grid
CMAQ Grid
Source: Gridded
CMAQ Estuary
Time
CMAQ Subestuary
Date(Y/M/D)| Src:Watershed NADP
Src:Watershed CMAQ
Retrieve & j src:Watershed Point
S re: Waters tied NonP
Load: Estuary Spa92
Ld:Est Spa02 MRB1
Ld:Est Spa02 MRB2
Ld:Est Spa02 MRB5
Ld:Est Spa02 MRB7
Ld:Subestuary Spa92
For Help- e Ld:Sub Spa02 MRB1
Ld:Sub SpaQ2 MRB2
Ld:Sub Spa02 MRB5
Ld:Sub Spa02 MRB7
Version 2Q1<
b)
Crops
Manure
Fertilizer
c)
Total Load
Total Yield
Total %
d)
Estuary Code
Crop Load
Fertilizer Load
Manure Load
Crop Yield
Fertilizer Yield
Manure Yield
Figure B-22. Submenus for nitrogen source data at gridded or estuarine scales.
B- 21
-------
Finally, estuarine watershed loads and yields of nitrogen can be explored based on the output of 1992
National and 2002 regional SPARROW models (http://water.usgs.gov/nawqa/sparrow/mrb/) that have
been aggregated to the estuarine watershed scale.
a
F Nitrogen j NADP Grid
• Satellite CMAQ Grid
• Pop. Dens Source: Gridded
CMAQ Estuary
CMAQ Subestuary
Date(Y/M/D) Src:Watershed NADP
Src:Watershed CMAQ
Retrieve & 5 src:Watershed Point
Src'Water'hed NonP
Load:Estuary Spa92
I Ld:Est Spa02 MRB2
Version 201-
Ld:Est Spa02 MRB7
Ld:Subestuary Spa92
For Help- e Ld:Sub Spa02 MRB1
-j Ld:Sub Spa02 MRB2
1 Ld:Sub Spa02 MRB5
Ld-Sub Spa02 MRB7
r
' [ Estuary Code
1 Predominant Source
' Total Yield
• Total Load
. ' Natural Load
Wet Deposition Load
Pert. Total Load
Pert. Corn/Soy Load
Pert. Alfalfa Load
i_fc Pert. Wheat Load
• Per. Farm-Other Ld
Forested Land Load
Barren Land Load
Shrub Land Load
b)
Estuary Code
Tntfll Yiairt
Total Loac Estuary Code
k> Fed fnrn Total Yil ' ' . 1
n Estuary Code I
Pert, Othe Total Lo Tntii v pctuarv r>dp
Aacp inori • Fertilizl ToLal Load
Developet Atm We Manure Fertilizer Load
Municipal Urban F Mariurg' Manure Load
MuniClp' Atmos Atmos. Dep. Load
Urban 1 Forest Alder Load
Forest East Load
Developed Land Load
Municipal Pt. Load
Canada Load
Boundary Load
Figure B-23. Selection of N loading data from SPARROW models.
B- 22
-------
B.5. Downloading Data
The data displayed in the view window as the result of reference boundary layers selected under Tab L
Zoom Maps, data selected under Tab 2 (Get Data) can be downloaded to the directory of the user's
choice in Tab 3. Save Data. Data can be downloaded in a variety of formats for later viewing (.png, .mpg,
.kml) or import into decision support modeling applications (shapefile, ASCII grid). Metadata are
automatically provided in associated .txt or .xml files. See the message box at the bottom for an
indication of when data downloads have been completed and Files Listing box for a list of files that have
been downloaded. Time series of remote sensing data can be saved either as a series of ASCII grids, one
for each date, or as a single shapefile. The latter output provides a more compact format for
downloading and subsequent data calculations; each date of data of the time series is provided in a
separate column of the associated dbf file. Users will need to pay attention to missing value indicators
(e.g., -9999) in using these data for later calculations.
1. Zoom Maps | J^^lM^- Save Dataja.Done
Choose Directory/Folder For Saved Files:
c:\ternp
Save Data (Shapefile) | Save Data (ASCII Grid) |
Save lrnage(png)| Save Movie(mpg) ] SaveKML
Save Map Polygons/DBF (VERY SLOW)
09:50 CllAQ_DEP.txt
09:50 CMAQJI.dbf
9:50 CHAQJI.pr:
9:50 CHAQ_M.shp
9:50 CHAQJI.shx
:SO CHAQJT.txt
pcj
9:50 estuaries shp
shx
9: 50
9:50
9:50
9:50
9:50
9:50
9:50
9:50
9:50
9:50
9:50
9:50
9:50
9:50
34 09:50
HEBRSj
(ERRS
dbf
shx
Civerall_Bounds.pr j
Overall_Bounds.shp
Overall_Eiound3. shx
states.prj
state3.shp
states.shx
watersheds.pcj
watersheds.shp
watersheds.shx
watersheds.xml
Finished retrieving data.
Retrieved 84 points of Water Quality
Retrieved 594 points of Deposition
For Help: eflrn@epa.gov 919-541-5500
Figure B-24. Save Data tab in Estuary Data Mapper, illustrating saving shapefiles.
B- 23
-------
B.6. Terminating an EDM Session and EDM Updates
Users should end an EDM session using Tab 4. Done. Exiting in this manner will save user settings so
that the next time EDM is started the user will automatically be zoomed in to the area of interest
selected in the last session run. The user can also set preferences for updating EDM versions on this tab
by selecting among the options given:
1. Zoom Maps | 2 Get Data 3. Save Data 4 Done
Quit |
Check for updated version (20140702):
• Never
• Start
r QUI
Check For Update Now
To learn more please visit our website:
http:/Aw/w epa.gov/edm/
Or contact edm@epa.gov
Finished retrieving data
Retrieved 84 points of Water Quality
Retrieved 594 points of Deposition
-------
Appendix C. R Packages and Commands
Used in Development of Predictive
Seagrass Habitat Models
C-l
-------
Convert seagrass polygons to 10m x 10m grids
Extract with 0 to 5 meter depth grid mask
Extract centroids from grid -> point shapefile
Join points with other attribute polygons (shapefiles) or extract values from attribute
grids to add predictor variables to seagrass point shapefile
Use NEAREST command to associate seagrass point with closest shoreline event (10-meter spacing)
Use NEAR command to calculate distance to nearest shoreline
Import dbf file into R for statistical analysis
Exploratory analysis
Run initial glrnmPQL models without SAC
Diagnostic tests {e.g., VIF) -> eliminate cross-correlated variables and interaction terms
Stepwise elimination of nonsignificant model predictors
Export dbf file with model residual errors
Join residuals with original shapefile using coordinates
Evaluate range of SAC along shoreline using spline correlograms
Evaluate range of SAC offshore using communitycorrelogram
Calculate zonal averages of residuals based on SAC ranges
Import dbf file to R
Rerun glmmPQL models with SAC term included
Export dbf file with model residual errors
Check for elimination of SAC using spline correlograms
Final model diagnostic checks
Figure C-1 General sequence of CIS and R analyses to create generalized linear mixed models to predict seagrass
presence/absence.
C-2
-------
In the following examples, R command lines are preceded by a ">" and are typed in
boldface. Responses from the console are shown in regular typeface.
C.I. Exploratory analysis (graphics package)
Spine plots
> # Explore linearity of binary responses
> spineplot(fsavcode~cSAL, data=nb9cc)
> spineplot(fsavcode~cUSRMARIkm2, data=nb9cc)
C.2. General linear mixed models with diagnostics plots (MASS, rms
packages)
> library (MASS)
># Define formula
> fo.glmmlS <- formula(fsavcode ~ cSAL + cSALZ + cSDavgrtrZ + cSDavgrt_l + cSDavgrt_2 + cptTOC +
cptTOCZ + cptTOCS + fZgtMXZavT + cCG046avkm + cDstoMarin + fISOLATED + fSED4 + PResidl4fa)
> # Run generalized linear mixed model
> GLMM15w.nb9cc <- glmmPQL(fo.glmml5, data = nb9ccl4fa!300x200, weights = weight, random
=~11 fSHORLIN, family = "binomial")
iteration 1
iteration 2
iteration 3
iteration 4
iteration 5
iteration 6
iteration 7
iteration 8
> # summarize results
C-3
-------
> summary(GLMM15w.nb9cc)
Linear mixed-effects model fit by maximum likelihood
Data: nb9ccl4fa!300x200
AlCBICIogLik
NA NA NA
Random effects:
Formula: ~1 | fSHORLIN
(Intercept) Residual
StdDev: 6.373279 0.9092619
Variance function:
Structure: fixed weights
Formula: ~invwt
Fixed effects: fsavcode ~ cSAL + cSAL2 + cSDavgrtrZ + cSDavgrt_l + cSDavgrt_2 + cptTOC + cptTOC2 +
cptTCO + fZgtMXZavT + cCG046avkm + cDstoMarin + fISOLATED + fSED4 + PResidl4fa
Value Std.Error DF t-value p-value
(Intercept) 0.4202571.6594588518856 0.25325 0.8001
cSAL 1.434985 0.1887963 518856 7.60071 0.0000
cSAL2 -0.5672100.1552607518856 -3.65327 0.0003
cSDavgrtrZ 2.612055 0.0536428 518856 48.69350 0.0000
cSDavgrt_l -0.2334610.0328527518856 -7.10632 0.0000
cSDavgrt_2 -0.0561920.0095881518856 -5.86057 0.0000
cptTOC -2.4438260.0858941518856-28.45162 0.0000
cptTCO 1.3156780.0707202518856 18.60401 0.0000
cptTCO -0.3165800.0547752518856 -5.77962 0.0000
C-4
-------
fZgtMXZavTl 0.5337630.1427521518856 3.739090.0002
cCG046avkm 0.0113500.0038698518856 2.93287 0.0034
cDstoMarin -0.0014540.0000522518856-27.86994 0.0000
fISOLATEDl -10.666443 1.0603103 518856 -10.05974 0.0000
fSED46 -1.119823 0.2133335 518856 -5.24917 0.0000
fSED47 0.6322200.0722058518856 8.75582 0.0000
fSED412 -0.833132 0.1516006 518856 -5.49558 0.0000
fSED4124 3.657733 0.3345729 518856 10.93254 0.0000
fSED4810 -0.3775930.1502137518856 -2.51370 0.0119
PResidl4fa 3.1729770.0560561518856 56.60364 0.0000
Correlation:
(Intr) cSAL cSAL2 cSDvgZ cSDv_l cSDv_2 cptTOC cpTOC2 cpTOC3
cSAL -0.063
cSAL2 -0.069 0.931
cSDavgrtrZ -0.007 0.093 0.048
cSDavgrt_l -0.036 0.043 0.038-0.324
cSDavgrt_2 0.032 -0.080 -0.065 -0.129 -0.844
cptTOC -0.007 0.170 0.141 -0.268 0.049 0.024
cptTOC2 -0.030 0.084 0.022 0.104 0.043 -0.074 0.069
cptTCO 0.003 -0.024 0.022 0.010 -0.048 0.045 -0.636 -0.427
fZgtMXZavTl -0.091 0.052 0.054 0.020 0.066 0.003 -0.008 0.013 -0.007
cCG046avkm -0.051 0.558 0.578 0.034 0.084-0.071 0.056 0.091 0.016
cDstoMarin -0.016 0.318 0.331 -0.095 -0.078 0.089 0.261 -0.060 0.000
fISOLATEDl -0.194 0.040 0.048 -0.092 0.063 -0.037 0.064 -0.034 0.006
fSED46 -0.005 -0.011 0.025 0.000 -0.088 0.106 -0.049 0.066 -0.003
C-5
-------
fSED47 -0.002 -0.239 -0.171 0.149 -0.030 -0.009 -0.097 0.103 -0.126
fSED412 -0.013 -0.062 -0.023 -0.039 0.024 -0.011 -0.109 -0.020 -0.074
fSED4124 0.003 -0.095 -0.085 0.063 -0.017 0.005 -0.097 0.031 -0.022
fSED4810 0.006-0.029-0.011 0.136-0.175 0.114 0.149-0.012-0.128
PResidl4fa 0.020 0.051-0.031 0.485-0.109-0.013-0.372 0.242-0.060
fZMXZT CCG046 cDstMr fISOLA fSED46 fSED47 fSED412 fSED4124 fSED48
cSAL
cSAL2
cSDavgrtrZ
cSDavgrt_l
cSDavgrt_2
cptTOC
cptTCO
cptTOC3
fZgtMXZavTl
cCG046avkm 0.114
cDstoMarin -0.006 0.296
fISOLATEDl -0.001 0.016 0.047
fSED46 -0.002 0.104 0.241 -0.026
fSED47 -0.022 -0.005 -0.089 -0.079 0.284
fSED412 -0.022-0.012-0.095 0.020 0.069 0.239
fSED4124 0.033-0.050-0.136-0.093 0.094 0.201 0.104
fSED4810 -0.031 0.004 0.245 -0.048 0.277 0.370 0.038 0.110
PResidl4fa 0.0330.037-0.243-0.1450.0190.152-0.063 0.107 0.118
C-6
-------
Standardized Within-Group Residuals:
Min Ql Med Q3 Max
-82.898446484 -0.049605288 -0.007560305 -0.001533338 361.324250323
Number of Observations: 518890
Number of Groups: 16
> fixed.effects(GLMM15w.nb9cc, family = binomial)
(Intercept) cSAL cSAL2 cSDavgrtrZ cSDavgrt_l cSDavgrt_2 cptTOC cptTOC2
cptTOC3
0.420257164 1.434985095 -0.567209800 2.612055480 -0.233461322 -0.056191589 -2.443825859
1.315678476 -0.316580172
fZgtMXZavTl cCG046avkm cDstoMarin fISOLATEDl fSED46 fSED47 fSED412
fSED4124 fSED4810
0.533762818 0.011349608 -0.001454164-10.666443279 -1.119823166 0.632220198 -0.833132486
3.657733496 -0.377592950
PResidl4fa
3.172976995
> # Output random effects
> random.effects(GLMM15w.nb9cc, family = binomial)
(Intercept)
-99 2.7822797
2 -9.7044680
3 -5.6673015
6 -9.6165676
7 -6.2225803
8 -4.1356439
9 0.4226314
C-7
-------
10 3.0034085
11 -2.6426397
12 1.4561820
13 -4.3301970
14 0.2790325
15 6.6477236
16 9.7049135
17 11.4413276
18 6.5818993
Diagnostics examples
> # Check assumption of collinearity using variance inflation factor
> library(rms)
>vif(GLMM15w.nb9cc)
> # Check for heterogeneity of variance and patterns in residuals
> nb9cc7fal300x200$PResid!4wfa <- residuals(GLMM14w.nb9cc, type = "pearson")
> nb9cc7fal300x200$Predict!4wfa <- predict(GLMM14w.nb9cc, data = nb9ccfa, type = "response")
> plot(x = nb9cc7fal300x200$Predict!4wfa, y = nb9cc7fal300x200$PResid!4wfa, main = "Pearson
Residuals vs Predicted")
> plot(x = nb9cc7fal300x200$cSDavgrtrZ, y = nb9cc7fal300x200$PResid!4wfa)
> plot(x = nb9cc7fal300x200$cSAL, y = nb9cc7fal300x200$PResid!4wfa)
> plot(x = nb9cc7fal300x200$cptTOC, y = nb9cc7fal300x200$PResid!4wfa)
> plot(x = nb9cc7fal300x200$cSAL, y = nb9cc7fal300x200$PResid!4wfa)
> plot(x = nb9cc7fal300x200$cSDavgrtrZ, y = nb9cc7fal300x200$PResid!4wfa)
> spineplot(fsavcode~Predict!4wfa, data=nb9cc7fa!300x200)
> # Export for calculating smoothed residual
C-8
-------
> # Fit spline correlogram
> # Create separate dataset for each shoreline to plot spline correlograms
>nb9cc7fal300x200.PResidl4w<-data.frame(nb9cc7fa!300x200[,c(7,15,16,106)])
>str(nb9cc7fal300x200.PResid!4w)
'data.frame': 518890 obs. of 4 variables:
$ fSHORLIN : int 111111 11 11111111 11 11...
$ ShLnDist : int 2815150 2815150 2815150 2815150 2815150 2815150 2815150 2815150 2815150
2815150...
$ Distance : int 0 0 0 0 0 0 0 0 0 0 ...
$ PResidl4wfa: atomic -0.0186 -0.0117 -0.0114 -0.0114 -0.0114 ...
..- attr(*, "label")= chr "Standardized residuals"
C.3. General linear mixed model with interaction plots (effects package)
> fo.glmSy <- formula(savcodeAv ~ csalAv + csalAvZ + csecchiminAv + csecchiminAvZ + cpttocAv +
cpttocAvZ + cZgtMXZavAv + cZgtMXZavAvZ + cDistToHdShAv + cDstoMarinaAv + cUSRMARIkmZAv +
fISOLATED + cwindAv * fSEDS + PR3xFA460)
> GLMSy.avgPAbyShLnDist <- glm(fo.glm3y, data = avgPAbyShLnDistwPResid3xFA460, family =
"binomial", weights = Nweight)
> # effects interaction plots for fSEDS * cwindAv
> library(effects)
> plot(effect("cwindAv:fSED5",GLM3y.avgPAbyShLnDist))
C.4. Generalized additive mixed models (mcgv package)
> fo.gamm4w4w <- formula(fsavcodeMa ~ s(csalAv) + s(csecchimin) + s(cpttocMin) +
csecchimin:cpttocMin + fZgtMXZavM + cCG046avkm + cDistToHdS + cDstoMarin + cwindMin * fSED4 +
PR4wFA1370)
> GAMM4wfa4w.maxPAbyShLnDist <- gamm(fo.gamm4w4w, data = maxPAbyShLnDistwPR4wFA1370m,
random=list(fSHORLIN=~l),weights = weight, family = "binomial")
C-9
-------
C.5. Spline correlogram evaluations (ncf package)
> # Spline correlogram
>load("D:\\savhabitat\\workingmodels\\NB9\\nb9ccl4fal300x200_ShLnl4.rda")
> library(ncf)
> # Randomly select 9000 to evaluate SA at coarser scales
> nb9ccl4fal300x200.PResidl5w.ShLn!4.sub9000 <-
nb9ccl4fal300x200.PResidl5w.ShLnl4.subset[sample(l:nrow(nb9ccl4fal300x200.PResidl5w.ShLn!4
.subset),9000);]
># Generate spline correlogram using shoreline distance as x coordinate and setting y (Distance) to 0
> fitl5wShLn!4 <- spline.correlog(nb9ccl4fal300x200.PResidl5w.ShLn!4.sub9000$ShLnDist,
nb9ccl4fal300x200.PResid!5w.ShLnl4.sub9000$Distance,
+ nb9ccl4fal300x200.PResid!5w.ShLnl4.sub9000$PResidl5wfa,
+ w = NULL, df = NULL, type = "boot",resamp = 100, npoints = 300, save = FALSE, filter = FALSE, fw = 0,
max.it = 25, xmax = 15000, latlon = FALSE, na.rm = FALSE,
+ quiet = FALSE)
1 of 100
2 of 100
99 of 100
100 of 100
> plot.spline.correlog(fitl5wShLn!4)
> # export elements of spline correlogram fit for plotting because ncf doesn't allow adjustment of y-
axis
> fitl5wShLn!4.predictedy <- fitl5wShLn!4$real$predicted$y
> write.csv(fitl5wShLn!4.predictedy, file =
"D:\\savhabitat\\workingmodels\\NB9\\fitl5wShLnl4predictedy.csv")
> fitl5wShLn!4.predictedx <- fitl5wShLn!4$real$predicted$x
C-10
-------
> write.csv(fitl5wShLn!4.predictedx, file =
"D:\\savhabitat\\workingmodels\\NB9\\fitl5wShLnl4predictedx.csv")
> fitl5wShLn!4.xint <- fitl5wShLn!4$real$x.intercept
> write.csv(fitl5wShLnl4.xint, file = "D:\\savhabitat\\workingmodels\\NB9\\fitl5wShLnl4xint.csv")
> fitl5wShLn!4.bootxint <- fitl5wShLnl4$boot$boot.summary$x.intercept
> write.csv(fitl5wShLn!4.bootxint, file =
"D:\\savhabitat\\workingmodels\\NB9\\fitl5wShLnl4bootxint.csv")
> # bootstrap summaries - plot rows 3, 6, and 9 for 5%ile, median and 95%ile
> fitl5wShLn!4.bootpredy <- fitl5wShLnl4$boot$boot.summary$predicted$y
> write.csv(fitl5wShLn!4.bootpredy, file =
"D:\\savhabitat\\workingmodels\\NB9\\fitl5wShLnl4bootpredy.csv")
> plot(fitl5wShLnl4.predictedx,fit!5wShLnl4.predictedy, type = "I", main = "Spline Correlogram", sub
= "Model GLMlSw shoreline 14 random subset of 9000", xlab = "Distance (m)",
+ ylab = "Correlation")
>lines(fitl5wShLnl4.predictedx,fitl5wShLnl4.bootpredy[3,],col="blue")
>lines(fitl5wShLnl4.predictedx,fitl5wShLnl4.bootpredy[9,],col="red")
> abline(h = 0)
> plot(fitl5wShLnl4.predictedx,fit!5wShLnl4.predictedy, type = "I", main = "Spline Correlogram", sub
= "Model GLM15.nb9cc Shoreline 14 random subset of 9000", xlab = "Distance (m)",
+ ylab = "Correlation")
>lines(fitl5wShLnl4.predictedx,fitl5wShLnl4.bootpredy[3,],col="blue")
>lines(fitl5wShLnl4.predictedx,fitl5wShLnl4.bootpredy[9,],col="red")
> abline(h = 0)
C.6. Community correlograms with anisotropy (CommunityCorrelogram
package)
# Extract dataframe with x-y coordinates;
PResid4w.xy <- nb9cc.PResid4w.ShLnl4.sub!200[,c(2,4)]
C-ll
-------
head(PResid4w.xy)
# Find maximum range of distances;
max(PResid4w.xy$DistToShor)
# Figure out optimum lag size and number;
# lagmin should be greater than smallest distance (10m)
# lagmax should be less than 2/3 maximum interpoint distance
lagSelect(sampleData=PResid4w,sampleLocation=cbind(PResid4w.xy,z=0),sampleTime=NULL,Location
Names=NULL,lagmin=ll,lagmax=250,by=30,option=l,plot=T,anisotropic=T,azimuth=90,azimuthTol=0,
bandwidth=0,dipAngle=0,dipTol=0,dipBandwidth=0,distmeth='euclidean')
commcorrelogram(sampleData=PResid4w,sampleTime=NULL,sampleLocation=cbind(PResid4w.xy,z=0)
,LocationNames=NULL,option=l,metric='manter,lagNumber=25,lagSize=ll,lagTol =
5.5,numTests=99,anisotropic=TRUE,azimuth=90,azimuthTol=0,bandwidth=0,dipAngle=0,dipTol=0,dipB
andwidth=0,distmeth='euclidean',mantmeth='spearman',adj='holm',prog=TRUE,alternative='one.side
d1)
C.7. Cross-validation and ROC construction (pROC package)
Subsetting data and generating training and test data sets
> library(MASS)
> # XVAL_100514.R Cross-validation for final models
> load(file = "D:\\savhabitat\\workingmodels\\NB9\\nb9ccl4fal300x200.RDA")
> load(file = "D:\\savhabitat\\workingmodels\\NB9\\GLMM15w_nb9cc.RDA")
> # First sample each random effects group (SHORLIN) separately and combine rows, then repeat 10
times
> df.O <- nb9cc!4fa 1300x200 # original data frame
> df.2 <- df.O[df.O$fSHORLIN == 2,]
> df.7 <- df.O[df.O$fSHORLIN == 7,]
> df.6 <- df.O[df.O$fSHORLIN == 6,]
> df.9 <- df.O[df.O$fSHORLIN == 9,]
> df.14 <- df.O[df.O$fSHORLIN == 14,]
C-12
-------
> df.8 <- df.O[df.O$fSHORLIN == 8,]
> df.10 <- df.O[df.O$fSHORLIN == 10,]
> df.12 <- df.O[df.O$fSHORLIN == 12,]
> df.13 <- df.O[df.O$fSHORLIN == 13,]
> df.3 <- df.O[df.O$fSHORLIN == 3,]
> df.17 <- df.O[df.O$fSHORLIN == 17,]
> df.ll <- df.O[df.O$fSHORLIN == 11,]
> df.15 <- df.O[df.O$fSHORLIN == 15,]
> df.18 <- df.O[df.O$fSHORLIN == 18,]
> df.16 <- df.O[df.O$fSHORLIN == 16,]
> df.99 <- df.O[df.O$fSHORLIN == -99,]
>
> # create subsample 1 of 1/10 each SHORLIN n set with replacement
> df.2.1 <- df.2[sample(l:163486,size = 16347, replace = TRUE),]
> df.7.1 <- df.7[sample(l:84233,size = 8423, replace = TRUE),]
> df.6.1 <- df.6[sample(l:63230,size = 6323, replace = TRUE),]
> df.9.1 <- df.9[sample(l:58307,size = 5831, replace = TRUE),]
> df.14.1 <- df.!4[sample(l:48586,size = 4859, replace = TRUE),]
> df.8.1 <- df.8[sample(l:44745,size = 4475, replace = TRUE),]
> df.10.1 <- df.lO[sample(l:22029,size = 2203, replace = TRUE),]
> df.12.1 <- df.l2[sample(l:11176,size = 1118, replace = TRUE),]
> df.13.1 <- df.!3[sample(l:5389,size = 539, replace = TRUE),]
> df.3.1 <- df.3[sample(l:4768,size = 477, replace = TRUE),]
> df.17.1 <- df.!7[sample(l:3189,size = 319, replace = TRUE),]
> df.11.1 <- df.ll[sample(l:3025,size = 303, replace = TRUE),]
C-13
-------
> df.15.1 <- df.l5[sample(l:2354,size = 235, replace = TRUE),]
> df.18.1 <- df.!8[sample(l:2135,size = 214, replace = TRUE),]
> df.16.1 <- df.l6[sample(l:1172,size = 117, replace = TRUE),]
> df.99.1 <- df.99[sample(l:1084,size = 108, replace = TRUE),]
>dftest.l<-
rbind(df.2.1,df.7.1,df.6.1,df.9.1,df.l4.1,df.8.1,df.l0.1,df.l2.1,df.l3.1,df.3.1,df.!7.1,df.ll.l,df.l5.1,df.l
8.1, df.16.1, df.99.1)
>dftest.l$TESTK-l
>TEST1.JOINID <- dftest.l[,c("JOINID","TESTl")]
> dfall.l <- merge(df.O,TESTl.JOINID,all.x = TRUE)
> dftrain.l <- dfall.l[-which(dfall.l$TESTl == 1),]
># subsample 10
> df.2.10 <- df.2[sample(l:163486,size = 16347, replace = TRUE),] # create subsample 1 of 1/10
SHORLIN 2 set with replacement
> df.7.10 <- df.7[sample(l:84233,size = 8423, replace = TRUE),]
> df.6.10 <- df.6[sample(l:63230,size = 6323, replace = TRUE),]
> df.9.10 <- df.9[sample(l:58307,size = 5831, replace = TRUE),]
> df.14.10 <- df.!4[sample(l:48586,size = 4859, replace = TRUE),]
> df.8.10 <- df.8[sample(l:44745,size = 4475, replace = TRUE),]
> df.10.10 <- df.lO[sample(l:22029,size = 2203, replace = TRUE),]
> df.12.10 <- df.l2[sample(l:11176,size = 1118, replace = TRUE),]
> df.13.10 <- df.!3[sample(l:5389,size = 539, replace = TRUE),]
> df.3.10 <- df.3[sample(l:4768,size = 477, replace = TRUE),]
> df.17.10 <- df.!7[sample(l:3189,size = 319, replace = TRUE),]
> df.11.10 <- df.ll[sample(l:3025,size = 303, replace = TRUE),]
> df.15.10 <- df.!5[sample(l:2354,size = 235, replace = TRUE),]
C-14
-------
> df.18.10 <- df.!8[sample(l:2135,size = 214, replace = TRUE),]
> df.16.10 <- df.l6[sample(l:1172,size = 117, replace = TRUE),]
> df.99.10 <- df.99[sample(l:1084,size = 108, replace = TRUE),]
>dftest.lO<-
rbind(df.2.10,df.7.10,df.6.10,df.9.10,df.l4.10,df.8.10,df.l0.10,df.l2.10,df.!3.10,df.3.10,df.l7.10,df.ll.
10,df.l5.10,df.l8.10,df.l6.10,df.99.10)
>dftest.lO$TEST10<-l
>TEST10.JOINID <- dftest.lO[,c("JOINID","TEST10")]
> dfall.10 <- merge(df.O,TEST10.JOINID,all.x = TRUE)
> dftrain.10 <- dfall.lO[-which(dfall.lO$TEST10 == 1),]
> # Save training and test data sets so that they can be selectively reloaded to conserve memory for
model runs
> save(dftrain.l,file = "D:\\savhabitat\\workingmodels\\NB9\\dftrainl.rda")
> save(dftest.l,file = "D:\\savhabitat\\workingmodels\\NB9\\dftestl.rda")
> save(dftrain.2,file = "D:\\savhabitat\\workingmodels\\NB9\\dftrain2.rda")
> save(dftest.2,file = "D:\\savhabitat\\workingmodels\\NB9\\dftest2.rda")
> save(dftrain.3,file = "D:\\savhabitat\\workingmodels\\NB9\\dftrain3.rda")
> save(dftest.3,file = "D:\\savhabitat\\workingmodels\\NB9\\dftest3.rda")
> save(dftrain.4,file = "D:\\savhabitat\\workingmodels\\NB9\\dftrain4.rda")
> save(dftest.4,file = "D:\\savhabitat\\workingmodels\\NB9\\dftest4.rda")
> save(dftrain.5,file = "D:\\savhabitat\\workingmodels\\NB9\\dftrain5.rda")
> save(dftest.5,file = "D:\\savhabitat\\workingmodels\\NB9\\dftest5.rda")
> save(dftrain.6,file = "D:\\savhabitat\\workingmodels\\NB9\\dftrain6.rda")
> save(dftest.6,file = "D:\\savhabitat\\workingmodels\\NB9\\dftest6.rda")
> save(dftrain.7,file = "D:\\savhabitat\\workingmodels\\NB9\\dftrain7.rda")
> save(dftest.7,file = "D:\\savhabitat\\workingmodels\\NB9\\dftest7.rda")
> save(dftrain.8,file = "D:\\savhabitat\\workingmodels\\NB9\\dftrain8.rda")
C-15
-------
> save(dftest.8,file = "D:\\savhabitat\\workingmodels\\NB9\\dftest8.rda")
> save(dftrain.9,file = "D:\\savhabitat\\workingmodels\\NB9\\dftrain9.rda")
> save(dftest.9,file = "D:\\savhabitat\\workingmodels\\NB9\\dftest9.rda")
> save(dftrain.lO,file = "D:\\savhabitat\\workingmodels\\NB9\\dftrainl0.rda")
> save(dftest.lO,file = "D:\\savhabitat\\workingmodels\\NB9\\dftestl0.rda")
> save.image("D:\\savhabitat\\workingmodels\\NB9\\101214a_xval.RData")
> # For each set run model with training set and predict results for test set
> # Repeat for sample 10;
> load("D:\\savhabitat\\workingmodels\\NB9\\dftrainl0.rda")
>load("D:\\savhabitat\\workingmodels\\NB9\\dftestl0.rda")
> # Run models first using original PResid zonal averages;
> # Run model foglmmlS with training sample 10 dftrain.10
> GLMM15w.dftrain.10 <- glmmPQL(fo.glmml5, data = dftrain.10, weights = weight, random
=~11 fSHORLIN, family = "binomial")
iteration 1
iteration 2
iteration 3
iteration 4
iteration 5
iteration 6
iteration 7
iteration 8
> # Generate predictions with test sample 10, type = response dftest.10
> dftest.lO$Predict!5 <- predict(GLMM15w.dftrain.lO, newdata = dftest.10, type = "response")
> # Generate raw residuals with test sample 10 after converting factor back to original numeric value
C-16
-------
> dftest.lO$Resid!5 <- as.numeric(levels(dftest.lO$fsavcode))[dftest.lO$fsavcode] -
dftest.lO$Predict!5
> save(dftest.lO,file = "D:\\savhabitat\\workingmodels\\NB9\\dftestlOxval.rda")
> # Combine predicted values, residuals and calculate mean residual, MSE sum
>load("D:\\savhabitat\\workingmodels\\NB9\\dftestlxval.rda")
>load("D:\\savhabitat\\workingmodels\\NB9\\dftest2xval.rda")
> load("D:\\savhabitat\\workingmodels\\NB9\\dftest3xval.rda")
> load("D:\\savhabitat\\workingmodels\\NB9\\dftest4xval.rda")
> load("D:\\savhabitat\\workingmodels\\NB9\\dftest5xval.rda")
> load("D:\\savhabitat\\workingmodels\\NB9\\dftest6xval.rda")
> load("D:\\savhabitat\\workingmodels\\NB9\\dftest7xval.rda")
> load("D:\\savhabitat\\workingmodels\\NB9\\dftest8xval.rda")
> load("D:\\savhabitat\\workingmodels\\NB9\\dftest9xval.rda")
> load("D:\\savhabitat\\workingmodels\\NB9\\dftestlOxval.rda")
> #Append test files
>dftest.ltolO<-rbind(dftest.l[,-113],dftest.2[,-113],dftest.3[,-113],dftest.4[,-113],dftest.5[,-
113],dftest.6[,-113],dftest.7[,-113],dftest.8[,-113],dftest.9[,-113],
+ dftest.lO[,-113])
>#Calculate summary statistics for residuals and squared residuals
> dftest.ltolO$Resid!5_2 <- dftest.ltolO$Resid!5 * dftest.ltolO$Resid!5
> summary(dftest.ltolO$Residl5)
Min. IstQu. Median Mean 3rd Qu. Max. NA's
-1.000000 -0.073900 -0.001545 -0.128500 -0.000043 1.000000 18
> summary(dftest.ltolO$Residl5_2)
Min. IstQu. Median Mean 3rd Qu. Max. NA's
0.000000 0.000000 0.000004 0.117600 0.009783 1.000000 18
C-17
-------
ROC statistics
> library(pROC)
Type 'citation("pROC")' for a citation.
Attaching package: 'pROC
The following objects are masked from 'package:stats':
cov, smooth, var
> roc(dftest.ltolO$fsavcode,dftest.ltolO$Predictl5, plot = TRUE)
Call:
roc.default(response = dftest.ltolO$fsavcode, predictor = dftest.ltolO$Predictl5, plot = TRUE)
Data: dftest.ltolO$Predict!5 in 506133 controls (dftest.ltolO$fsavcode 0) < 12759 cases
(dftest.ltolO$fsavcode 1).
Area under the curve: 0.7144
C-18
-------
Appendix D. Exploratory Analyses and
Diagnostic Tests
D-l
-------
D.I. Models for Seagrass Grid Cell Presence/Absence
D.I.I. Exploratory Analyses
Exploratory analyses were conducted to evaluate the distribution of predictor variables,
e.g., to determine if there are potential problems with outliers. Figures Dla-r illustrate the
distribution of continuous variables used in model development. With the exception of
wave energy, variables exhibit few influential outliers. If necessary, a square root
transformation could be applied to the wave energy variable to even out the distribution;
however that would change the form of the response to wave energy.
Conditional plots for the binary response variable as a function of predictors (spine plots]
provide information on potential nonlinearities in response (before the effect of other
variables are factored out}. Spine plots show conditional probabilities of seagrass presence
as a function of predictor variables, segmented into groups (Figures D2a-l}. The width of
each band is inversely proportional to the number of observations in that band. The
inverted light color bands at the top of plots represent the relative frequency of occurrence
of current (2006} eelgrass and the darker gray band segments at the bottom represent the
relative frequency of current eelgrass absence1. Most of the continuous potential
predictors appear to be related linearly to the relative frequency of seagrass occurrence,
with two exceptions. The response to wind-generated wave mixing depth appears to be
unimodal, with maximum response at intermediate values. The effect of distance from
historic seagrass patch edge appears to drop off exponentially rather than linearly.
Seagrass probability of occurrence appears to vary among shoreline segments. Counter to
our initial hypothesis, probability of seagrass appears to be greater along isolated shoreline
segments. Seagrass shows evidence of persistence, with greater probability of occurrence
in areas of historic occurrence, and lessening likelihood of occurrence as distance from
historic patch edges increases.
1R also includes the cdplot function which will produce a smoothed version of conditional
probability plots, but these are less useful because they can include artifacts where there is
sparse representation of points along a gradient.
D-2
-------
b)
8-
c)
12345
Sediment total organic carbon (%)
26.0 26.5 27.0 27.5 28.0 28.5
Salinity (PSD)
8-
14.7
14.8 14.9
Temperature (deg C)
15.0
e)
-20 -15 -10 -505
Depth - max wave mixing depth (m)
S §-
0 20000 40000 60000 80000
Representative wave energy (joules/m)
50 100 150 200 250 300
•v'Windwave energy (Joules/m)
Figure D-l. Distribution of independent continuous variables entered into general linear models and general additive
models: a) sediment percent total organic carbon; b) salinity (PSU); c) temperature (deg C); d) Depth minus max wave mixing
depth (m); e) representative wave energy (joules/m); f) square root wind wave energy (Joules/m); g) Avg Secchi depth minus
water depth; h) Minimum Secchi depth minus water depth; i) Canada goose (numbers/km2); j) logio (Canada goose
density/km2 +1); k) distance to hardened shoreline (m); I) ^distance to hardened shoreline (m); m) density of unsewered
residences on high infiltration soils (no/km2); n) distance to nearest marina (m); o) A/distance to nearest marina; p) area of
1999 seagrass patch (m2); q) distance to edge of nearest 1999 seagrass patch (m); r) Iogl0(max wave mixing depth/water
depth).
D-3
-------
g)
h)
CD
U) O
=! CO
O
CT3
n
o.
<•£>
U- O
CN
i)
4-20 2 4
Average Secchi depth - Depth (m)
j)
-4-2024
Minimum Secchi depth - Depth (m)
o
1/T m
-a
c
03
ID
u
09
IT 8.1
,
in
T3
O)
=!
cr
o
o-
o
Q
CN
o
LO -
o
o-
O
Li'j
n
Q
0 5 10 15 20 25 30 35
Avg 2004-2006 Canada geese/km2
o
' r--
8-
V)
=! O
O Tf
O &~> ~
c
0)
I)
0.0 0.5 1.0 1.5
Log.0(Avg 2004-2006 Canada geese/km2)
o.
m
CM
o O
-S= CN
1=-. t
o
o- o.
0 500 1000 1500 2000 2500
Nearest hardened shoreline (m)
0 10 20 30 40
\Nearest hardened shoreline (m)
50
Figure D-l (cont'd)
D-4
-------
m)
n)
o
o
CO
18-
£ CN
8-
ri
0>
^
cr
0) o
-m
n_
0 50 100 150 200 250 300
Unsewered residences on high infiltration soils/km2
500 1000 1500 2000
Distance to nearest marina (m)
2500
o
o-
-
-So
£= O
to co
o
o
10 20 30 40
\Distance to nearest marina (m)
50
p)
s
q)
o
O-i
CD O
13 T—
o
O
LL if)
l~l-i-
0 2000 4000 6000 8000 10000
Distance to nearest 1999 eelgrass patch (m)
Figure D-l (cont'd)
-3-2-10 1 2
Log-,j(max wave mixing depth/water depth)
D-5
-------
b)
6 7 11 1213
Sediment class 2
810
-15 -10 -5 05 15 25
Centered wind wave energy (kJoules/m)
-1.6 -1-0.8 -0.6-0.4-0.20 0.2 0.611.63.4
Centered sediment total organic carbon (%)
e.
-4 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.54
Centered avg Secchi depth-water depth (m)
Q.
0 1
Depth > average mixing depth (T = 1)
Figure D-2. Spine plots showing conditional probabilities of seagrass presence as a function of potential independent
predictor variables. Width of band is inversely proportional to number of observations in band. Inverted light color bands at
top of plots represent relative frequency of occurrence of current (2006) eelgrass and darker gray band segments at bottom
represent relative frequency of current eelgrass absence : a) Sediment particle-size class (5 = sand, 6 = gravelly-sand, 7= silty-
sand, 11 = clay-silt, 12 = sand-silt-clay, 13 = gravel-silt-clay, 1+2+4 = gravel + sandy gravel + gravel-sand-silt, 8+10=
silty+sandy-silt); b) Sediment percent total carbon; c) centered salinity (PSU); d) centered average Secchi depth -water
depth (m); e) centered wave energy(kJoules/m); f) dummy variable indicating water depth in excess of wave mixing depth; g)
dummy variable indicating isolated shoreline segment in middle of channel; h) shoreline segment in Narragansett Bay (see
Figure 3); i) dummy variable indicating historic (1996-97) eelgrass presence/absence; j) centered distance to historic (1996-97)
eelgrass edge; k) distance to hardened shoreline (km); I) area of 1999 seagrass patch.
D-6
-------
g)
CD
tn
-S3
ID
a>
0
•>
in
Q.
in
CD
m
t>
CO
GO
o"
i)
CD
o
CD
O
(J
c
Q)
ID
CO
CO
ID
k)
s.
•>
p
CD
Isolated shoreline (T =1}
op
o
ID
O
CN
o
o
o
Q 1
1999 Eelgrass presence (1) or absence (0}
OJ
u
CO
CD
D)
ro
CO
CO
o
CD
O
j)
CN
o
o
o
-99
6 7 8 9 10 13 15
Shoreline code
o.
ro
o
c
CD
CO
"CD
o
tu
if!
CD
tn
to
CD
ta
cq
o
CD
o
o
c
CD
in
0)
tn
CD
0>
CO
4 -3 -2 -.5 .5 1.5 2.5 3.5
'} Centered distance to nearest 1999 eelgrass
o
o
4.5
CO
o
C-J
t=i
CD
o
CD
O
CD
in
CD
cn
CD
tfl
-.7 -.5 -.4 -.3 -.2 0.1 .3 .5 .8
Centered distance to nearest hardened shoreline (km)
-20000 0
Centered area of 1999 eelgrass patch (m2)
Figure D-2. (Cont'd)
D-7
-------
D.1.2. Evaluation of Model Assumptions in Preliminary Seagrass Grid P/A Models
We performed diagnostic tests to check on validity of model assumptions for the
preliminary seagrass grid presence/absence models, e.g., multicollinearity, heterogeneity
of variance, random distribution of model residuals, and spatial independence of model
residuals. We evaluated an initial generalized linear mixed model (GLMM) (Table 3} to
predict seagrass presence at the scale of 10-meter grid cells, with shoreline included as a
random effect and all potential predictors (Table D-l] and interaction terms included. To
avoid problems with multi-collinearity of independent variables we modified the initial
model to remove all interaction terms except for Wave Energy x Sediment Type. A review
of the working model residuals showed multiple problems with heterogeneity of variance
and potential nonlinearities (Figure D-3a-c). In general, values of predictor variables
associated with high probability of seagrass occurrence also showed a greater range of
Pearson residuals. (Pearson residuals are adjusted to account for higher expected variance
with the mean, so should show no pattern when graphed against predicted values or
against each predictor.}
Diagnostic tests showed evidence of spatial autocorrelation for working model residuals. A
run of the model using the quasibinomial family in place of the binomial family
demonstrated that errors were underdispersed, with a dispersion factor of 0.71. Normally
a binomial model is assumed to have a dispersion value of one because, by definition, the
variance is equal to the mean. The spline correlogram constructed based on Pearson
residuals for a random subset of 9000 observations from the working model showed
evidence of spatial autocorrelation at fine scales, with only small levels of autocorrelation
at coarser scales (Figure D-4a}. We first used horizontal swaths of 9000 points from each
of the main branches of Narragansett Bay to fine-tune the estimate of the range of spatial
autocorrelation at short spatial scales. The lower limit of the 95% confidence intervals of
the first x-intercept from bootstrapping spline fits ranged from 345 to 3428 meters (Figure
D-4b). Note that reconstruction of spline correlograms from horizontal swaths suggested
that spatial autocorrelation occurred not only at fine scales but also at intermediate and
extreme distance ranges. Because these did not appear in the correlogram constructed
from points sampled from the entire range, we assumed that these zones of apparent
spatial autocorrelation were an artifact of subsampling.
As explained earlier, we used spline correlograms based on shoreline distance to evaluate
the range of spatial autocorrelation parallel to shorelines, but used Mantel correlograms
(using the R commcorrelogram package] to evaluate the range of autocorrelation along
depth gradients perpendicular to the shoreline. Due to memory constraints, we ran the
communitycorrelogram function for a random subset of points for one shoreline at a time.
Figure D-5 shows the Mantel correlogram for Shoreline 14. Spatial autocorrelation was not
significant at distances greater than 200 meters in the offshore direction.
D-8
-------
Table D-l. Potential independent variables included in models predicting seagrass presence.
Variable
Definition
Units
fSHORLIN
cSAL
cTEMPER
fSEDn
cWIND
cPTTOC
cSDavggtrZ
cSDmngtrZ
fZgtrMXZ
cZgtrMXZ
fISOLATED
cDistHdShor
cDistMarina
cUSRMARIkm2
cCG046avkm2
fEGPA99
cAREA
cDistEG99
Shoreline code -99,1-19
Centered growing season salinity PSU
Centered growing season average water Deg C
temperature
Sediment type (n represents level of 1-13
lumping]
Wind wave energy Joules/m
Centered sediment percent total organic %
carbon
Centered growing season average Secchi Meters
depth - water depth
Centered growing season minimum Secchi Meters
depth - water depth
Depth greater than wave mixing depth (0 =
FALSE, 1 = TRUE]
Centered water depth greater than wave Meters
mixing depth
Isolated shoreline (0 = FALSE, 1=TRUE)
Centered distance to hardened shoreline Meters
Centered distance to nearest marina Meters
Centered unsewered residences on high #/km2
infiltration soils/catchment area
Centered winter 2004-2006 Canada goose #/km2
density
Historic (1999} eelgrass presence (0 =
FALSE, 1 = TRUE]
Centered area of 1999 eelgrass patch Meters2
Centered distance to edge of nearest 1999 Meters
eelgrass patch
D-9
-------
a)
s •
re
•S o
I s
19
O .
II
b)
o .
o .
n
3 o
•o °
I/I
Ss-<
go
^
£
o
6
o .
• MM*
-1.5
-1.0
Centered
-0.5
salmitv(PSU)
0.0
0,5
—i— —i— —r—
-1.0 -0.5 0.0 0.5 1.0 1.5 2.0
Centered sediment total organic carbon (%)
-2024
Centered average Seccht Depth -water depth (m)
Figure D-3. Heterogeneity of variance for residuals of initial GLMM model to predict seagrass presence at the scale of 10
meter grid cells. Residuals show evidence of heterogeneity of variance.
D-10
-------
a) b)
o
in
o"
j_
o
to O _
cu o
t
o
O
9"
p
I f^
\ ^ /A
\
P
ID
o"
f—
O
To o
t °
o
c^
ID
^3 ~
J S
i i i i i i i i ^
0 5 10 15 20 25 30 v"'
^
T=ii_j.
^^%=. _^l^ii=*^ie-J5ai ,
^~=J"g*""i"F ^^^f.-^^^^'L^ jaJ^^gtf
1 1 1 1 1 1
Distance [km) ° 2 4 G 8 10
Distance (km)
Figure D-4. Spline correlogram with 95% confidence intervals from bootstrapping based on Pearson residuals from working
version of general linear mixed model, a) Random subset of 9000 points from all shorelines, b) Random subset of 9000
points from shoreline 14.
D-ll
-------
o
"co
§
co
"c
03
E
00
o
p
o
o
o
o
o
p
o
50
100
150
200
250
lag distance
I
I
'c
CT
'co
CO
O
CD
O
CN
O
O
O
empirical
alpha
0
50
100
150
200
250
lag distance
Figure D-5. Correlogram showing spatial autocorrelation for offshore distance based on residuals from model predicting
seagrass presence/absence for Shoreline 14. Spatial autocorrelation is not significant at distances greater than 200 meters.
D-12
-------
D.2. Models for Shoreline Segment Presence/Absence
D.2.1. Exploratory Analyses of Shoreline P/A Relationships
Figures D-6a-p display exploratory spine plots showing general relationships between
single predictors and seagrass presence/absence by shoreline distance. Percent shoreline
occupied by seagrass ranges from less than 1 percent (shoreline 3} to almost 60%
(shoreline 17], with occupancy rates almost twice as high for isolated as compared to main
shorelines (Figures D-6a, b}. Shoreline occupancy ranges from almost 60 to over 80
percent for zones near larger 1999 eelgrass patches, and falls off exponentially with
distance to nearest 1999 eelgrass patch (Figures D-6d}. As a single predictor Secchi depth
performs poorly along the entire gradient, although occupancy rates are greater than 60%
for the highest decile of minimum transparency values (Figures D-6e, f). Shoreline
occupancy shows an apparent unimodal response to salinity and percent sediment total
organic carbon, but does not appear to vary with average temperature (Figures D-6g, h, i}.
Occupancy increases steadily as the proportion of depths greater than average wave mixing
depths increases. Shoreline occupancy is low for gravel- dominated sediment classes (1-2-
4 and 6} (Figure D-61}. Occupancy is highest for the top four deciles of distance to
hardened shoreline (Figure D-6o}.
D-13
-------
b!
CD
O
-99 2 367 S 9 10 13 14 15 18
Shoreline code
Isolated shoreline (T = 1)
-20000 0
Max centered seagrass patch area (m2)
e)
-2.5 -2.0 -1.5 -1.0 -0.50 1.5 2.5 4.0
Min centered distance to 1999 eelgrass (km)
.4-1-0.6 0 0.20.6
Max centered avg Secchi depth (m)
-1 -0.4 -0.2 0 0.20.61
Max centered seasonal min Secchi depth (m)
Figure D-6. Spine plots showing conditional probabilities of seagrass shoreline presence (1) versus absence (0) as a function
of potential independent predictor variables: a) Shoreline (see Figure 3); b) shoreline isolation; c) centered maximum 1999
eelgrass patch size (m2); d) centered min distance to 1999 eelgrass patch (km); e) maximum centered time-averaged Secchi
depth at i (m); f) maximum centered seasonal-minimum Secchi depth (m); g) centered average salinity (PSU); h) minimum
centered average temperature (deg C); i) minimum centered sediment percent total organic carbon; j) minimum centered
unsewered residences on high infiltration soils/km2 watershed area; k) minimum centered wind energy (kJoules/m); I)
maximum centered depth - average wave mixing depth; m) dominant sediment particle-size class; n) minimum centered
average Canada goose density (geese/km2); o) maximum centered distance to hardened shoreline; and p) maximum
centered distance to marina (km). Width of band is proportional to number of observations in band. Inverted light color
bands at top of plots represent relative frequency of occurrence of current (2006) eelgrass and darker gray band segments at
bottom represent relative frequency of current eelgrass absence at shoreline unit i. For continuous variables, values are
grouped into subsets by decile to facilitate viewing patterns.
D-14
-------
g)
1
-1.8-1.7-0.5 0.1 0.2
Average centered salinity (PSU)
co
CD
o
'.£>
O
CN
O
Win centered temperature (deg C)
-1.2 -0.8 -0.6 -0.4-0.20 0.2 0.4).6 1 1.4
Min centered sediment total organic carbon (%)
•40 -20 02040180
Min centered unsewered residences on hi inf n soils/km2
CO
CD
O
(N
O
CO
-10 -5 0 5 15 40
Min centered wind wave energy (kJoules/m)
Figure D-6 (Cont'd).
p
o
-4 -2 -1.5 -0.5 0.5 1 1.5 2 2.5
Max centered depth - avg wave mixing depth (m)
D-15
-------
ml
6 7 13 8/10
Sediment class
p) -4 -2 0 2 6 20
Min centered 2004-06 Canada geese/km2
-0.6 -0.4 -0.2 0 0.2 0.61.00
Max centered distance to hardened shoreline (km)
-0.4 -0.20.41.22.2
Maximum centered distance to marina (km)
Figure D-6 (Cont'd)
D.2.2. Test of Model Assumptions During Initial Model Development for Shoreline P/A
Prior to fitting a model for shoreline occupancy, we screened predictors for potential cross-
correlations. Minimum shoreline average temperature was correlated with maximum
shoreline average Secchi depth and shoreline average salinity (r > 0.70}. Shoreline average
salinity was correlated with maximum shoreline average Secchi but not with maximum
shoreline Secchi depth minima. Thus we dropped temperature and maximum shoreline
average Secchi depth, but retained maximum shoreline Secchi depth minimum as
predictors.
D-16
-------
The original model fit after dropping insignificant terms and prior to accounting for spatial
autocorrelation included the following main effects: salinity, availability of depths greater
than average mixing depth, shoreline isolation, distance to hardened shorelines, distance to
nearest marina, density of Canada geese, and two interaction terms: wind-generated wave
energy x sediment particle- size class and Secchi depth x sediment percent total organic
carbon. The density of unsewered residences on high infiltration soils was not a
significant predictor. Model diagnostic tests showed evidence of heterogeneity of variance,
and spline correlograms with 95% confidence interval by shoreline showed evidence of
strong spatial autocorrelation (Figures D-7a, b}. The minimum x-intercept, corresponding
to the minimum range at which spatial autocorrelation of residuals was undetectable, was
1365 meters. Subsequent peaks in the correlogram have similar breadth, and probably
represent autocorrelation between different seagrass patches, as compared to correlation
within a given patch. In many cases, residuals appeared to be greater for classes or ranges
associated with shoreline occupancy (Figures D-7c, d}.
D.2.3. Test of Model Assumptions During Initial Model Development for Shoreline Relative
Abundance
The final model fit to predict average shoreline occupancy was more complex,
incorporating both higher-order terms to account for nonlinear responses and two- and
three-way interaction terms:
savcodeAv = csalAv + csalAv2 + csalAv3 + csecchiminAv + csecchiminAv3 +
cpttocAv2 + cpttocAv3 + cZgtMXZavAv + cZgtMXZavAv2 + cDistToHdShAv + (Dl)
cDstoMarinaAv + cUSRMARIkm2Av + flSOLATED + csalAv x cpttocAv +
csecchiminAv x cpttocAv + csalAv x csecchiminAv x cZgtMXZavAv + cwindAv x
fSEDS)
with terms defined as above (Table D2) and added variables defined as:
fSEDS = sediment classes
SED5_124561213 = Gravel, Sandy gravel, Gravel-sand-silt, Sand, Gravelly
sand, Sand-silt-clay, Gravel-silt-clay
SED5_7 = Siltysand
SED5_81011 = Silty, Sandy Silt, and Clay-Silt
D-17
-------
Based on diagnostic plots of Pearson residuals versus individual predictor variables,
nonlinearities appear to have been accounted for in the final model. However, the model
still demonstrated heterogeneity of variance and spatial autocorrelation of residuals
(Figure D-7a-d).
bi
75 co
CL
o
J5 p
OJ CD
o
O
0.2 0.4 0.6
Predicted shoreline presence
5 10
Distance (km)
15
M-
0.
"™ CO-
TS
C/5
£ <£>-
c
0
(n ^._
CD
•1
CL
-------
o
O
c)
0.0 0.2 0.4 O.S
Predicted shoreline occurrence
C.S
d)
5 10
Distance (km)
15
. SI
I
i«:.|i
11*
-0.5 0.0 0.5 1.0 1.5 2.0
Centered average Secchi depth minimum (m)
-1.0 -0.5 0.0 0.5 1.0 1.5 2.0
Centered percent total organic carbon
Figure D-8. Diagnostics for model 3j. a) Pearson residuals versus predicted value, b) spline correlogram of Pearson residuals
for Shoreline 14. The minimum x-intercept (within 95% confidence interval) was 420 meters. Higher order terms appear to
have accounted for most of the nonlinearities in relationships, as illustrated by plots of Pearson residuals versus c) centered
average Secchi depth minimum (m) and d) centered shoreline average percent sediment total organic carbon.
D-19
-------
D.2.4. Test of Model Assumptions During Initial Model Development for Maximum Seagrass
Depth
Models predicting minimum and maximum depth of occurrence by shoreline index were
developed for the subset of records corresponding to shoreline occupancy. Models
predicting minimum depth of occurrence were created with a subset of data after excluding
records with average or maximum mixing depths of zero. The restriction of model
predictions to shoreline occurrence transects allowed us to drop the random shoreline
effect from models and compare the fit of general linear models with general additive
models because fitting GAMs is much less memory intensive than fitting GAMMs.
The initial GLM model evaluated to predict maximum depth of seagrass occurrence
included main effects and both two- and three-way interaction terms for shoreline
maximum Secchi depth average, shoreline maximum percent sediment total organic
carbon, and density of unsewered residences on high infiltration soils:
bathymmax = csecchiavMax + cpttocMax + cUSRMARIkm2Max + csecchiavMax x
cpttocMax + csecchiavMax x cUSRMARIkm2Max + cpttocMax x
cUSRMARIkm2Max + csecchiavMax x cpttocMax x cUSRMARIkm2Max
(D2)
After we dropped insignificant interaction terms, the resultant model was
bathymmax = csecchiavMax + cpttocMax + cUSRMARIkm2Max + csecchiavMax x
cpttocMax + csecchiavMax x cUSRMARIkm2Max
(D3)
Residual plots showed an increasing variance with the mean, possibly peaking at
intermediate values, so we fit a generalized additive models and general linear models with
second- and third-order terms to try to capture nonlinearities in response. The best GLM
model with higher order terms and interactions was:
bathymmax = csecchiavMax2 + csecchiavMax3 + cpttocMax + cpttocMax2 +
cUSRMARIkm2Max + cUSR
csecchiavMax x cpttocMax
(D4)
cUSRMARIkm2Max + cUSRMARIkm2Max2 + cUSRMARIkm2Max3
However, this model still showed evidence of heterogeneity of variance, so we refit GLM
models after loglO-transformation of predictors. No terms were dropped from the log-
transformed model with 3-way interactions and VIF terms were less than 10:
D-20
-------
bathymmax = cL10secchiavMax*cL10pttocMax*cL10USRMARIkm2Max) (D5)
However, based on comparison of AIC values, a GAM model provided a superior fit
compared to the GLM model with log-transformed predictors:
bathymmax = s(cL10secchiavMax] + s(cLlOpttocMax) + s(cL10USRMARIkm2Max)) (D6)
where s = smoothing function
D-21
-------
Appendix E. Quality of the Data and
Limitations on Use of the Data
E-l
-------
We used our Narragansett Bay pilot application of the statistical modelling approach
to illustrate how a predictive model could be developed to assess factors affecting
seagrass growth. Our goal was to distinguish between nutrient and nonnutrient
factors affecting seagrass growth and survival and to elucidate different
mechanisms of action for effects of nutrients on seagrass. We used the best
publically available data sets to describe environmental variables that affect
seagrass growth and survival in Narragansett Bay to support our predictive
statistical model based on: temporal matches to 2006 seagrass maps, spatial extent,
completeness, and spatial resolution. We filled in gaps in Secchi depth at the
southern end of the western arm of Narragansett Bay using offshore remotely
sensed estimates of light attenuation coefficients (from > 30 meters depth],
assuming these were similar to nearshore values. When only point data were
available, e.g., for salinity, Secchi depth, and wave energy, we created a continuous
grid by interpolation through use of Theissen polygons, the Euclidean function in
ArcMap (filling in gaps along shoreline], or inverse distance weighting. Some
variables were not available so we substituted indicators. This included use of
density of coastal residences on high infiltration soils as an indicator of potential
groundwater N inputs, use of distance to nearest marina as an indicator of potential
physical disturbance from mooring beds, and use of salinity gradients as an
indicator of total N gradients in final scenarios. Data on some variables potentially
affecting seagrass distribution were simply not available, e.g., sulfide concentrations
in sediment porewater, actual measurements of tidal current velocity as compared
to estimated values, and measurements of turbidity or light attenuation near the
sediment interface (as opposed to upper water column]. Any of the limitations to
data availability or completeness of model inputs could have influenced the
accuracy of our model predictions, but are unlikely to have produced biased model
results. We based model projections of future condition following nutrient load
reductions on the assumption that nutrient concentrations will decline in
proportion to load reductions and that space-for-time substitutions are appropriate
for model development. The model can be improved in the future as more complete
data or modelled estimates become available, e.g., for tidal currents and for
dissolved inorganic N and total N concentrations across the bay.
E-2
-------
United States
Environmental Protection
Agency
PRESORTED STANDARD
POSTAGE & FEES PAID
EPA
PERMIT NO. G-35
Office of Research and Development (8101R)
Washington, DC 20460
Official Business
Penalty for Private Use
$300
------- |