Sea-level rise and other influences on decadal-scale salinity variability in a coastal plain estuary lable at ScienceDirect Estuarine, Coastal and Shelf Science 157 (2015) 79e92 Contents lists avai Estuarine, Coastal and Shelf Science journal homepage: www.elsevier.com/locate/ecss Sea-level rise and other influences on decadal-scale salinity variability in a coastal plain estuary Andrew C. Ross a, *, Raymond G. Najjar a, Ming Li b, Michael E. Mann c, Susan E. Ford d, Brandon Katz a, 1 a Department of Meteorology, 503 Walker Building, The Pennsylvania State University, University Park, PA 16802, USA b Horn Point Lab, University of Maryland Center for Environmental Science, P.O. Box 775, Cambridge, MD 21613, USA c Department of Meteorology and Earth and Environmental Systems Institute, The Pennsylvania State University, University Park, PA 16802, USA d Haskin Shellfish Research Laboratory, Rutgers University, 6959 Miller Avenue, Port Norris, NJ 08349, USA a r t i c l e i n f o Article history: Received 22 July 2014 Accepted 31 January 2015 Available online 5 March 2015 Keywords: estuaries salinity climate climatic changes sea level changes USA Delaware Delaware Estuary * Corresponding author. E-mail addresses: andrewross@psu.edu (A.C. R (R.G. Najjar), mingli@umces.edu (M. Li), mann@psu.e rutgers.edu (S.E. Ford). 1 Present address: JLT Towers Re, New York, NY 10 http://dx.doi.org/10.1016/j.ecss.2015.01.022 0272-7714/© 2015 Elsevier Ltd. All rights reserved. a b s t r a c t The response of salinity in the Delaware Estuary to climatic variations is determined using statistical models and long-term (1950-present) records of salinity from the U.S. Geological Survey and the Haskin Shellfish Research Laboratory. The statistical models include non-parametric terms and are robust against autocorrelated and heteroscedastic errors. After using the models to adjust for the influence of streamflow and seasonal effects on salinity, several locations in the estuary show significant upward trends in salinity. Insignificant trends are found at locations that are normally upstream of the salt front. The models indicate a positive correlation between rising sea levels and increasing residual salinity, with salinity rising from 2.5 to 4.4 per meter of sea-level rise. These results are consistent with results from 1D and dynamical models. Wind stress also appears to play some role in driving salinity variations, consistent with its effect on vertical mixing and Ekman transport between the estuary and the ocean. The results suggest that continued sea-level rise in the future will cause salinity to increase regardless of any change in streamflow. © 2015 Elsevier Ltd. All rights reserved. 1. Introduction Salinity influences both the physical properties of an estuary and the characteristics of its ecosystem. For example, salinity is the dominant factor regulating stratification. Even small changes in the salinity of an estuary can have a significant impact on the estuary's ecosystem. For example, salinity influences the spread of oyster disease (Powell et al., 1992), the distribution and diversity of ammonia-oxidizing bacteria (Bernhard et al., 2005), and the development of phytoplankton blooms (Gallegos and Jordan, 2002). Understanding and mitigating the impacts of changing salinity are particularly important because climate change and other human activities have already stressed many estuarine eco- systems (Kennish, 2002). oss), najjar@meteo.psu.edu du (M.E. Mann), susan@hsrl. 017, USA. Many climatic and oceanic factors, including streamflow, sea level, oceanic salinity, and wind stress, have an influence on the salinity and water quality of an estuary. Streamflow determines the amount of fresh water entering the estuary. Elevated streamflows are typically associated with fresher water in the estuary; lower streamflows are associated with increased salinity in the estuary. Higher sea levels increase salinity by bringing more salt water into the estuary. Variations in oceanic salinity alter the salinity of water circulating into the estuary. Finally, wind stress may influence salinity through vertical mixing, Ekman transport and upwelling (Banas et al., 2004), and other mechanisms. Climate change has the potential to cause changes in all of these variables. Precipitation amounts, frequencies, and intensities are expected to change in many areas as a result of anthropogenic climate change, and the associated effects on streamflow may be complicated by land use and evaporation changes (Krakauer and Fung, 2008). Global mean sea level has risen significantly during the twentieth century and is expected to rise at an increasing rate through the twenty-first century (Rahmstorf, 2007; Vermeer and Rahmstorf, 2009; Church et al., 2013). Meanwhile, changes in Delta:1_given name Delta:1_surname Delta:1_given name Delta:1_surname Delta:1_given name Delta:1_surname mailto:andrewross@psu.edu mailto:najjar@meteo.psu.edu mailto:mingli@umces.edu mailto:mann@psu.edu mailto:susan@hsrl.rutgers.edu mailto:susan@hsrl.rutgers.edu http://crossmark.crossref.org/dialog/?doi=10.1016/j.ecss.2015.01.022&domain=pdf www.sciencedirect.com/science/journal/02727714 http://www.elsevier.com/locate/ecss http://dx.doi.org/10.1016/j.ecss.2015.01.022 http://dx.doi.org/10.1016/j.ecss.2015.01.022 http://dx.doi.org/10.1016/j.ecss.2015.01.022 A.C. Ross et al. / Estuarine, Coastal and Shelf Science 157 (2015) 79e9280 land use and large-scale atmospheric circulation have slowed winds over much of the Northern Hemisphere land area (Jiang et al., 2009; Vautard et al., 2010), although wind speeds have also increased in some areas (Hartmann et al., 2013). Regardless of the causes, salinity change could be detrimental to many estuaries. This study focuses on the salinity of the Delaware Estuary on the United States East Coast. Over 8 million people live within the Delaware River basin (Sanchez et al., 2012), and the estuary contains the largest freshwater port in the United States (Philadelphia) (Kauffman et al., 2011). The Delaware River and Es- tuary provide a significant amount of freshwater to New York City and Philadelphia. Salt intrusion into the Philadelphia area water supply can occur during periods of high salinity (Hull and Titus, 1986). These factors prompt the Delaware River Basin Commis- sion to carefully regulate the position of the salt front in the estuary. In addition, species in this estuary are typically sensitive to salinity; and interactions between parasites and hosts, as well as predators and prey, are often influenced by salinity. For instance, eastern oysters, a keystone species in the estuary, are host to two parasites that cause important diseases: Perkinsus marinus (dermo disease) and Haplosporidium nelsoni (MSX disease). The distribution of both parasites, as well as that of the oyster, is restricted by low salinity, but oysters can tolerate much lower salinity than either parasite, thereby providing low-salinity refuges from disease in the upper estuary (Haskin and Ford, 1982; Bushek et al., 2012; Ford et al., 2012). Because of the importance of the estuary and river for shipping, drinking water, and fishing, a number of studies have examined the physical properties of the estuary. Salinity is higher in the center of the estuary and lower near the shores. The estuary is weakly to partially stratified, and the lateral salinity difference is typically greater than the vertical difference (Wong and Münchow, 1995; Wong, 1995). However, significant vertical stratification can occur in the main channel (Aristiz�abal and Chant, 2013). The estuary experiences two tides per day as a result of a large principal lunar semidiurnal (M2) constituent (Wong, 1995). Salinity variability in the estuary produced by tidal advection is larger than the vari- ability caused by streamflow (Garvine et al., 1992). Sea level and circulation also vary on the subtidal time scale primarily as a result of wind forcing (Wong and Garvine, 1984). Several numerical modeling studies have examined the response of salinity to sea-level rise (Hull and Tortoriello, 1979; U.S. Army Corps of Engineers Philadelphia District (1997); Kim and Johnson, 2007). These studies found that salinity should increase in response to sea-level rise in most of the estuary. Numerical models have produced similar results in other estuaries, including the Chesapeake Bay (Hilton et al., 2008) and the San Francisco Bay (Cloern et al., 2011). Although numerical model simulations can be informative, they are also subject to potentially restrictive assumptions and are no substitute for long-term observations of salinity trends. For example, all modeling studies to date assume that sea-level rise has no influence on bottom topography, even though it is likely that sea-level rise causes increased shoreline erosion, which increases sediment deposition (Cronin et al., 2003). Thus, empirical methods are an essential complement to numerical models for determining the effects of climate change and sea-level rise on salinity. Ordinary linear regression is commonly applied to empirically model salinity. For example, Garvine et al. (1992) and Wong (1995) used linear regression to model the response of the salt intrusion length to streamflow in the Delaware Estuary. Marshall et al. (2011) used multiple linear regression to build predictive models of salinity in the Florida Everglades. However, when applying linear regression, care must be taken to account for issues such as cor- relation among data (autocorrelation), non-constant variance (heteroscedasticity), and non-linearity that are commonly found in water quality data (including salinity data). Autoregressive models have been applied to empirically model salinity by taking advantage of the highly autocorrelated nature of most water quality data. Using autoregressive models, Gibson and Najjar (2000) predicted the response of salinity in the Ches- apeake Bay to future changes in streamflow, and Hilton et al. (2008) tested whether sea-level rise has caused significant changes in Chesapeake Bay salinity. Saenger et al. (2006) used autoregressive models to relate river discharge to salinity and to reconstruct Ho- locene discharge and precipitation in the Chesapeake Bay watershed. Other studies have applied additive models to empirically model salinity. The additive model and the related generalized additive model expand the traditional linear regression model by modeling the response variable with one or more smooth functions with forms that are nonparametric (i.e., are not defined a priori). Several authors have recently applied these models in studies of salinity and other water quality metrics. Jolly et al. (2001) and Morton and Henderson (2008) used additive and generalized ad- ditive models to determine changes in river salinity. Autin and Edwards (2010) applied additive models to extract tidal variations from salinity, dissolved oxygen, and temperature data and found that the additive methods performed better than multiple regression. Neither additive nor autoregressive models offer a complete solution to the problems of autocorrelation, non-linear relation- ships, and heteroscedasticity commonly found in water quality data. Additive models are not typically robust against correlated or heteroscedastic errors, and autoregressive models do not handle heteroscedasticity or non-linear relationships between variables. Additive mixed models (AMMs) offer a solution to the compli- cations commonly encountered when modeling water quality time series. AMMs combine the nonparametric smooth functions of additive models with the ability to handle correlated errors and observations. AMMs are popular in many environmental fields that deal with autocorrelated and non-linear data, such as air pollution (Coull et al., 2001) and paleoclimatology (Simpson and Anderson, 2009). In this work, AMMs are applied to perform a data-driven analysis of the effects of climatic variations on salinity in the Delaware Estuary. 2. Methods 2.1. Study area and data The Delaware Estuary is located in the Eastern United States to the east of the Chesapeake Bay (Fig. 1a). The Delaware River is the primary source of river discharge to the estuary. The head of tide extends to Trenton. Salt intrusion normally extends through the lower half of the estuary (Garvine et al., 1992). Salinity in the Delaware Estuary has been measured through many different monitoring programs, including surveys, automated sensors, and boat sensors. However, records with long-term data coverage are rare. The goal of this analysis was to determine which variables have an influence on salinity over long time periods. In addition, the statistical models that were applied perform better with larger amounts of data. As a result, of the many salinity datasets that are available, the automated sensor data from the United States Geological Survey (USGS) and bottom salinity mea- surements from the Haskin Shellfish Research Laboratory (HSRL) were selected for statistical modeling. Both datasets are suitable for studying long-term trends, as they include a large number of measurements and together cover the period from the 1950s to the present. Fig. 1. a: Map of the United States East Coast. Black box indicates the region shown in b. b: Map of the Delaware Estuary, including the sea-level measurement site at Atlantic City and the locations where the USGS measured salinity and streamflow. Black box indicates the region shown in c. c: Locations of the oyster beds where the Haskin Shellfish Research Laboratory measured bottom salinity. Names correspond to the IDs in Table 1b. A.C. Ross et al. / Estuarine, Coastal and Shelf Science 157 (2015) 79e92 81 The USGS has measured near-surface salinity at five locations in the Delaware Estuary since the 1960s (Table 1a, Fig.1b). The salinity data are reported as daily averages that are computed from instantaneous measurements every 15 min. The measurements from Reedy Island Jetty contain the least amount of missing data. The USGS discontinued measurements at Ship John Shoal in 1986. Chester, Fort Mifflin, and Ben Franklin Bridge contain a large amount of missing data, particularly during winter months. When this analysis was performed, the USGS had approved the accuracy of salinity measurements through 2012, so values from Table 1 a: Summary of USGS salinity data. Columns provide station name, axial distance (distance from the mouth of the estuary), percent of non-missing monthly means during 1964e2012, and mean salinity for the five USGS salinity stations in the es- tuary. Mean salinity was calculated from all available monthly means. Axial distance was obtained from the Delaware River Basin Commission. b: Summary of oyster bed salinity data. Columns provide oyster bed ID, approximate axial and lateral distance, and mean of all available bottom salinity measurements for each oyster bed. Axial and lateral distances are approximated based on distances from the lines shown in Fig. 1c. (a) Name Axial dist. (km) Data coverage (%) Mean salinity Ben Franklin Bridge 161 77 0.12 Fort Mifflin 147 31 0.15 Chester 133 74 0.22 Reedy Island 87 87 4.4 Ship John Shoal 60 28 13 (b) ID Axial dist. (km) Lateral dist. (km) Number of obs. Mean salinity ARN 34.87 �1.82 409 11.5 MID 28.62 �0.93 303 13.6 COH 26.64 �2.98 643 13.9 SHJ 26.00 �0.44 178 14.8 SHR 21.68 �1.07 596 15.0 BEN 16.91 �2.82 491 17.0 NPT 14.69 �5.11 120 16.3 HGS 13.83 �3.45 193 17.0 NWB 13.27 �2.91 572 17.5 LDG 9.83 0.41 226 19.6 BDN 9.29 �5.91 329 17.7 EIS 7.69 �3.89 316 19.1 the earliest possible date through 2012 were used. Salinity mea- surements from all locations were converted from electrical con- ductivity to practical salinity units using the algorithm introduced by Lewis and Perkin (1981) and simplified by Schemel (2001). To focus on long-term variations and to reduce the computational time needed to fit the statistical models, the salinity data were converted from daily means to monthly means by averaging any month with at least 15 days of data. Plots of the annual cycles and anomalies in these data are shown in Fig. 2. Salinity is typically highest in late summer and early autumn and lowest in mid-spring, and streamflow follows the opposite pattern. Similarly, streamflow anomalies were lowest and salinity anomalies highest during the extensive drought in the 1960s and reversed during the wetter 1970s and 2000s. The accuracy of the USGS salinity measurements should be sufficient for statistical analysis. The earlier USGS measurements were made with a flow-through monitor. The accuracy stated by the USGS documentation for electrical conductivity (salinity) measurements from the flow-through monitor is ±3% of the full scale (Gordon and Katzenbach, 1983). The USGS also issues an annual water data report for each location, which classifies the measurements into four accuracy categories. During recent years, the USGS has typically categorized the accuracy as ±3�10% or ±10�15%. Even ±3% may be too pessimistic, however, as Katzenbach (1990) determined that the flow-through monitor was more accurate than other systems and that the monitor performed better than the stated accuracy. Furthermore, some USGS locations, including Reedy Island, have recently switched to a YSI Incorpo- rated sonde that has a stated salinity accuracy of 0.1 or 1%, whichever is greater (Mark R. Beaver, personal communication, March 5, 2007). Assuming these errors are random, they will be absorbed by the model residual and will not bias the analysis. Relocations are more difficult to account for and could affect the analysis. The USGS station at Chester moved 0.8 km upstream in April 1981, and the Ben Franklin Bridge station moved 0.09 km upstream in July 1988. Despite these relocations, the USGS has approved these data, so no attempt was made to correct for the effects of the relocations. Salinity was also measured at various oyster beds in the Dela- ware Bay by the HSRL. Haskin (1972) and Haskin and Ford (1982) Fig. 2. a: Mean annual cycles of salinity and streamflow. The annual cycles were calculated using all data after 1970 to exclude the unusual drought in the 1960s. b: Time series of anomalies. A.C. Ross et al. / Estuarine, Coastal and Shelf Science 157 (2015) 79e9282 report some results from this measurement program. From the 1950s to about 1980, HSRL salinity measurements were made by titration; thereafter, conductivity was measured using an “Autosal” Laboratory Salinometer and converted to practical salinity units. The sampling frequency was irregular; many beds were sampled once or twice per month, and sampling occurred more often in the warm season. The locations of the 12 oyster beds used in this analysis are shown in Fig.1c, and river distances and mean salinities are provided in Table 1b. Whereas the USGS measures surface salinity, the HSRL measured bottom salinity. Because the mea- surements were made intermittently from about 1950 to 1990, they were not reduced to monthly averages. Tide stage, streamflow, and day of year are included in the statistical models used to analyze these data to minimize the impact of the intermittent and instan- taneous sampling frequency. Daily averages of streamflow in the Delaware River at Trenton, NJ were obtained from the USGS. At Trenton, the flow is approxi- mately 58% of the total discharge into the estuary from land (Sharp, 1983; Sharp et al., 1986). Measurements were also obtained from the Schuylkill River near Philadelphia, PA, which accounts for an additional 15% of the total discharge. The Schuylkill gauge, which is upstream of the entrance of the river into the Delaware River, has a drainage area of 4903 km2 compared to the total of 4952 km2 in the watershed. To approximate the actual flow into the Delaware River, the Schuylkill gauge measurements were multiplied by the ratio of the total drainage area to the gauged drainage area and added to A.C. Ross et al. / Estuarine, Coastal and Shelf Science 157 (2015) 79e92 83 the corresponding day's average streamflow at Trenton (except when analyzing the salinity at the Ben Franklin Bridge, which is upstream of the Schuylkill River). No other river or stream con- tributes more than 1% to the total discharge (Sharp et al., 1986). For use in modeling monthly mean salinity at the five USGS stations, the daily total streamflows were converted to monthly means. Before modeling the instantaneous oyster bed salinity measurements, the daily streamflow measurements were smoothed with an exponential moving average with a weight co- efficient of 1/15. This accounts for the slow response of the estuary to streamflow. The weight coefficient was determining by maxi- mizing the model's log-likelihood. Goodness of fit measures also indicated that applying an exponential moving average produced better fits than other methods, such as applying a simple lag. Monthly averages of sea level at Atlantic City, NJ were obtained from the Permanent Service for Mean Sea Level (Holgate et al., 2013). Atlantic City was selected because it is just outside of the estuary. To match the monthly mean sea level with the instanta- neous oyster bed salinity measurements, the mean sea levels were interpolated to the days of the oyster bed salinity measurements using cubic splines. Because the oyster bed salinities were measured throughout the tidal cycle, it is also necessary to also account for tidal fluctuations in sea level. However, measurements of sea level with sufficient temporal resolution are not available for most of the time period of the oyster bed salinity measurements. Instead, the water level at Ship John Shoal, which is approximately in the middle of the axial length of the oyster beds, was approximated using harmonic con- stituents obtained from the National Ocean Service. Although the time difference between a high or low tide and the corresponding salinity measurement is provided in the HSRL data, the actual time of measurement is not. Therefore, to match the salinity measure- ments with a time and water level, the water level predictions were offset by the provided time difference. Then it was assumed that salinity measurements would have only taken place during the day (8AMe8PM), so the appropriate offset high or low water level in this range was selected. In the event that measurements could have occurred at either 8AM or 8PM, the two water levels were averaged. The subtidal mean sea levels discussed in the previous paragraph are not correlated with the tidal water levels, so it is permissible to include both in the statistical models. The accuracy of the oyster bed data is unknown, and errors may be present in the data, particularly since the data were read from punch cards. To eliminate outliers, the local outlier factor method (Breunig et al., 2000) was applied. Here, the outlier factor was based on how isolated an observation's salinity, log-streamflow, and sea-level values are compared to a minimum of the 15 near- est points in the three dimensional salinity-streamflow-sea-level space. The 1% of the data (43 observations) with the worst outlier factors were discarded. Although removing these outliers improved the error distributions obtained after applying the sta- tistical models later in the analysis, the results of the analysis were not significantly different. The time series of salinity at each oyster bed after removing outliers are shown in Fig. 3. As would be ex- pected, the patterns in the oyster bed bottom salinity are similar to those in the USGS surface salinity data. Salinity is typically lowest in spring and highest in autumn, and the impact of the extensive drought in the 1960s is immediately visible in the anomaly time series. Long-term measurements of offshore salinity outside of the Delaware Bay are not available. Instead, the Gulf Stream Index (Taylor, 1995), which represents the first principal component of the latitude of the north wall of the Gulf Stream, was used as a proxy for oceanic salinity. Lee and Lwiza (2008) determined that the index is a suitable proxy for salinity in the Mid-Atlantic Bight. The monthly data from 1966 to 2013 were obtained from http:// www.pml-gulfstream.org.uk/. Wind speed and direction were obtained from the North American Regional Reanalysis (Mesinger et al., 2006), which has a horizontal resolution of 32 km. Reanalysis data are advantageous because they contain no missing data or instrument biases and because they provide complete coverage over water (although they do so by blending observations with imperfect models). Wind data were taken from three reanalysis grid points over the Delaware Bay between 38.8 and 39.63�N, 74.9e75.6�W. 3-hourly wind speed and direction from the reanalysis were used to compute the meridional and zonal components of the wind stress with the equation t ! ¼ C10r ����U10��! ����U10��!, where C10 is a drag coefficient, r is the density of air, ����U10��! ���� is the wind speed at 10 m, and U10��! is the wind vector at 10 m. C10, which varies with wind speed, was calculated using the equation from Wu (1982) assuming constant air density. The wind stress components and magnitude were then averaged over the bay and by month to form monthly averages. Finally, alongshore and cross-shore wind stresses were calculated using an alongshore di- rection of south-southwest to north-northeast (the orientation of the estuary's mouth) and a cross-shore direction of east-southeast to west-northwest. Wind stresses over the shelf were also computed but not used because they were nearly identical to stresses over the Bay at the monthly time scale. The reanalysis data only cover 1979 through the present. When relating salinity to wind stress, any salinity measurements before 1979 were dropped. 2.2. Statistical models The influences of observed streamflow, sea level, wind stress, and oceanic salinity on estuarine salinity were extracted using AMMs. A separate AMM was used to model each USGS location. The oyster beds are relatively similar to each other, so one AMM was used to model all of the oyster beds together with a term included to account for random variability between beds. The basic model for surface salinity at each USGS station was Si ¼ b0 þ b1Xi þ fQðQiÞ þ fMðMonthiÞ þ εi (1) where Si is the ith monthly-mean salinity value, b0 is a constant intercept, fQ(Qi) is a spline that relates salinity to streamflow and is evaluated at the streamflow value Qi, and fM(Monthi) is a cyclic spline for capturing seasonality in salinity that is not explained by the other independent variables and is evaluated at the ith month. b1Xi is an optional term used to test the importance of another variable Xdfor example, sea level. This term is analogous to that in an ordinary or multiple linear regression model. The residual εi is assumed to follow a Gaussian distribution with zero mean and variance s2L, i.e. εi � Nð0; s2LÞ, where L functions as a weight that accounts for autocorrelation and heteroscedasticity. Note that an ordinary regression model assumes εi � Nð0; s2Þ, and therefore does not account for autocorrelation or heteroscedasticity. For the oyster bed bottom salinity data, the basic model was Sij ¼ bi þ b0 þ � b1Xij � þ b2Hij þ b3xi þ b4yi þ fQ � Qij � þ fDD � DDij � þ εij (2) where the subscript j denotes an individual observation at the ith oyster bed.In thismodel,bi isa uniqueinterceptforeachoysterbed; it accounts for random variability in background mean salinity mea- surementsateachoysterbed. Thisvariabilitycouldoccurasa resultof other physical factors that are not in the model such as depth, http://www.pml-gulfstream.org.uk/ http://www.pml-gulfstream.org.uk/ Fig. 3. a: Mean annual cycles of bottom salinity at the oyster beds. b: Time series of salinity anomalies. A.C. Ross et al. / Estuarine, Coastal and Shelf Science 157 (2015) 79e9284 proximitytounmodeledfreshwatersources,orlocalcirculation.Allof the b parameters are analogous to the parameters in a multiple linear regression model. b0 is a constant intercept that applies to every bed. b1Xij is an optional term for some predictor X. b2 gives the slope of the response of salinity to tidal water level and Hij is the water level predicted from harmonics. b3 is the axial salinitygradient and xi is the relative axial distance for the ith oyster bed (column 1 in Table 1b). Lateral salinitygradients can be large in the Delaware Estuary (Wong, 1994), so b4yi accounts for the lateral salinity gradient and distance (column 2 in Table 1b). fQ(Qij) is a smooth function of exponential moving averaged streamflow. fDD(DDij) is a cyclic spline that relates salinity to decimal day. Finally, the residuals εij � Nð0; s2LiÞ, where Li includes only heteroscedasticity. To test the influence of sea level, oceanic salinity, and wind stress on the salinity measurements, additional terms representing these variables were inserted as the optional terms in the two basic models (Equations (1) and (2)). For example, to test the influence of sea level on the USGS salinity data, an additional term b1Hi was added to Equation (1). This term works like a traditional linear regression model, and is known as a parametric term because it has a specified form. In this example, b1 represents a slope that models the linear response of salinity to sea level. When adding additional terms to the model, all available data were used (in other words, pairwise deletion of missing values was applied). The time trend and sea level slope were tested for both the USGS salinity moni- toring locations and for the oyster beds. The influences of oceanic A.C. Ross et al. / Estuarine, Coastal and Shelf Science 157 (2015) 79e92 85 salinity and wind stress were only tested at the USGS salinity monitoring locations. 2.2.1. Smoothing functions The terms fQ, fM, and fDD in Equations (1) and (2) are smoothing functions that enable the response of salinity to the predictor var- iables to be modeled nonparametrically. For example, fQ models the influence of streamflow on salinity. The shape and amount of smoothness of the functions are determined by the model fitting algorithm. The USGS locations are notably different, so fQ and fM are determined separately for each location. The oyster beds have sufficiently similar responses to streamflow and residual season- ality, so all of the oyster beds were grouped together to determine fQ and fDD. To make the model algebraically fit, the expected value of every smoothing function is set to zero (Wood, 2006). This is a purely algebraic concern and has no practical influence on the re- sults other than centering the smoothing functions at zero. Thin plate regression splines (Wood, 2003) were applied to model the response of salinity to streamflow. These splines are considered optimal for use in additive models (Wood, 2006). To account for residual seasonal effects, a cyclic cubic spline was used to model the response of salinity to the current month (for the USGS monthly averages) or to the decimal day of the year (for the HSRL instantaneous measurements). A cubic spline is a piecewise curve composed from a number of cubic polynomials. The locations where the polynomials are joined are known as knots. At each knot, the values and first and second derivatives of the joining poly- nomials are equal. For cyclic cubic splines, the first and last knots also have equal values and first and second derivatives. This makes a cyclic spline useful for modeling data where the response should be similar at the boundaries of the predictor variable. The first and last knots in the smoothing function were placed at months 1 and 13 or decimal days 0 and 1, which causes the spline to smooth continuously from December to January. The remaining knots were placed with even spacing by the model fitting algorithm. 2.2.2. Distributions The distribution for the residuals is εiðjÞ � Nð0; s2LðiÞÞ, where the subscripts in parenthesis are used in Equation (2) but not Equation (1). The L s are defined to accommodate any autocorrelation or heteroscedasticity, and the subscript i in the second distribution indicates that the residuals are independent from values at other oyster beds. For both the USGS and HSRL salinity data, the variance was assumed to depend on a power function of streamflow such that varðεiðjÞÞ ¼ s2 ���QiðjÞ ���2d (Pinheiro and Bates, 2000). Other vari- ance structures including homoscedastic errors, variances related to the fitted values, and exponential relationships were also considered. Experiments indicated that the power function of streamflow produced the highest likelihoods for most models. Furthermore, not including the fitted values in the variance func- tion allows the use of exact procedures to find the parameter d (Pinheiro and Bates, 2000). In addition to heteroscedasticity, temporal autocorrelation is also present in the USGS data. This autocorrelation is a common problem in hydrological and climatological time series (Hirsch and Slack, 1984). A first-order autoregressive (AR1) error process was used to model the temporal autocorrelation in the USGS data. The lag-1 correlation was estimated by the model fitting algorithm. The USGS locations are modeled separately, so spatial autocorrelation is not an issue. There was typically enough time between successive observations in the HSRL data that both spatial and temporal autocorrelation were assumed negligible. Additive mixed models assume a Gaussian distribution for the errors; however, the errors may be modeled as following another distribution in the exponential family, such as the gamma or Poisson distribution. In this case, the model is known as a gener- alized additive mixed model. The gamma distribution may seem like the best option for modeling salinity since both are positive definite. However, tests comparing the model performances and the distributions of the residuals indicated that the standard ad- ditive mixed model with a Gaussian distribution was by far the best choice for modeling salinity. One consequence of this choice is that the models can predict negative salinities. In practice, however, this rarely occurred. 2.2.3. Fitting and testing Additive mixed models can be fit using maximum likelihood estimation (MLE) methods (Wood, 2006). However, MLE generally produces biased estimates of variance. A modification to MLE known as restricted maximum likelihood estimation (REML) solves this problem. However, REML makes model selection difficult because two models that have been fit with REML can only be compared under certain restrictive conditions. Mixed models contain terms for both fixed effects and random effects. Fixed effects are model parame- ters that apply to the entire population being sampled. Random effects apply to an individual unit or group that was randomly taken from the population. For example, the unique intercept term bi in Equation (2) is a random effect because the value for bi is different for each oyster bed and is assumed to be drawn from a random normal distribution. The term b1 is a fixed effect because it applies to all observations from all oyster beds. The smooth terms such as fQ are split into both fixed and random effects (Lin and Zhang, 1999; Wood, 2004). REML works by removing the fixed ef- fects from the likelihood maximization procedure (Corbeil and Searle, 1976). As a result, it is only possible to compare models fit with REML when the fixed effects are identical (Pinheiro and Bates, 2000; Wood, 2006). To resolve this issue, the significance of the model terms was tested using models fit with MLE. The resulting best model was then re-fit with REML to produce the final results. Likelihood ratio tests were applied to determine the significance of the fixed effects. The likelihood ratio test works under the assumption that twice the difference of the log-likelihoods of two nested models has a known c2 distribution (Wilks, 1938; Pinheiro and Bates, 2000). Model selection for the unique intercept (bi in Equation (2)) and other random effects, such as a separate trend for each oyster bed, was performed by including all possible fixed effects, fitting with REML, and using likelihood ratio tests to determine which random effects to include. Finally, the fixed effects were tested using MLE as above and the final best model was re-fit using REML. These methods indicated that a random intercept for each oyster bed significantly improved the model. However, other random effects did not improve the model. The models were fit and tested using the mgcv version 1.7e28 package of the open source statistical software R (Wood, 2014). Details of the methods used in this software are described in Wood (2003, 2004, 2006). 3. Results 3.1. Basic models Streamflow is often one of the primary influences on estuarine salinity. In the Delaware Estuary, the annual cycles of salinity and streamflow are clearly opposed (Fig. 2); maximum streamflow and minimum salinity both occur in April, and minimum streamflow in August precedes peak salinity in September and October. Similarly, salinity was unusually high during the drought of the mid 1960s and subsequently declined during the period of increased Fig. 5. Seasonal variations in salinity. For the oyster beds, the smooth plots the rela- tionship between salinity and decimal day of year (term fDD in Equation (2)). For the USGS data, the smooths show the relationships between salinity and month of year (term fM in Equation (1)). The gray shaded regions indicate ±2 standard errors. A.C. Ross et al. / Estuarine, Coastal and Shelf Science 157 (2015) 79e9286 precipitation and streamflow in the 1970s. There is also a modest downward trend in salinity after 1980 corresponding to a moderate increase in streamflow during this time. Equations (1) and (2) model this relationship between salinity and streamflow. The fitted smooth functions of streamflow fQ, without the optional term b1Xi(j), are shown in Fig. 4. Actual salinity predictions from the model are derived by adding the values from this smooth function with the intercept and other terms in the model; therefore, these smooth functions of streamflow represent only the modeled influence of streamflow on salinity. Smooth terms such as fQ are specified to have an expected value of 0, so negative values for the smooths are entirely valid. The smooths show that salinity and streamflow are negatively correlated, as expected. The magnitude of the marginal response of salinity to streamflow is larger under low-flow conditions. The residual seasonal variation included in the model appears similar at all locations (Fig. 5). In general, after accounting for the influence of streamflow, salinity is lowest in May and June and highest in October and November. These residual seasonal terms are not significantly changed after including the additional model terms that are discussed in the following paragraphs. The basic terms included in the oyster bed model also give the response to tidal water level and the axial and lateral salinity gra- dients. Likelihood ratio tests indicated that all three terms improved the model. The response of salinity to tidal water level, b2, is 1.1 m �1, which is smaller than the response to subtidal (monthly-mean) sea level (presented later in Table 2). This is reasonable, since the estuarine salinity field may not fully adjust to sea level over the tidal time scale but should over the monthly time scale. The axial salinity gradient (b3 ¼ �0.28 km�1) and lateral salinity gradient (b4 ¼ 0.35 km�1) are similar in the oyster bed re- gion. The axial gradient of bottom salinity is slightly weaker than the �0.337 km�1 calculated by Garvine et al. (1992) from five bottom and near-surface monitoring stations spread throughout Fig. 4. Relationship between streamflow and salinity. The black lines denote term fQ in Equations (1) and (2). The shaded gray regions indicate ±2 standard errors. the estuary. The lateral distance was defined as positive away from the eastern shore of the bay, so the positive lateral gradient in- dicates that salinity increases away from the eastern beds and to- wards the center of the bay, which is consistent with Wong (1994). The adjusted R2 value for the model fit to Reedy Island salinity is 0.73. Farther down the estuary, the adjusted R2 value for the fit to salinity is 0.57 at Ship John Shoal and 0.83 at the oyster beds. Up- stream at Chester, Fort Mifflin, and Ben Franklin bridge, the adjusted R2 values are 0.18, 0.60, and 0.55 respectively. The reduced performance at these upstream locations is a result of higher variability of salinity and the weaker relationship between salinity and streamflow under low-flow conditions. As a result of the higher variability and weaker relationship, streamflow has less predictive power under low-flow conditions, and any statistical model based on streamflow would have reduced performance. The high salinity variability may simply be the result of some sort of random or in- ternal variability, or a process that happens at longer or shorter time scales than are captured by the model. It is also possible that extremely low flows amplify the effect of sea level, wind, and other factors. For example, although the axial salinity gradient in the estuary is mostly constant, salinity levels far upstream exhibit a Table 2 Trends in streamflow-adjusted salinity and response of salinity to sea level. Bold indicates results that are significant at the 95% confidence level. Location Trend (decade�1) Response to sea level (m�1) Ben Franklin Bridge 8.9 � 10�4 (0.30) 3.8 � 10�3 (0.68) Fort Mifflin 2.5 � 10�3 (0.31) 3.3 � 10�2 (0.17) Chester 2.1 � 10�4 (0.99) 2.6 � 10�2 (0.30) Reedy Island 0.17 (1.4 � 10�2) 3.3 (5.2 � 10�5) Ship John Shoal 1.6 (5.6 � 10�2) 4.4 (3.4 � 10�2) Oyster beds 0.38 (<1.0 � 10�5) 2.5 (6.2 � 10�3) A.C. Ross et al. / Estuarine, Coastal and Shelf Science 157 (2015) 79e92 87 more gradual transition to riverine salinity (Ketchum, 1952). Under low-flow conditions, the region of gradual transition could shift farther upstream, increasing the axial salinity gradient in this re- gion and possibly the effect of sea level, wind, and other factors. The model performance could be improved if these variables had a predictable change in effect under low-flow conditions or if there was some other variable that is important only under low-flow conditions and could explain this variability. However, we did not find any variable that did. Finally, we note that the poor perfor- mance at Chester may also be caused by the gauge's proximity to Chester Creek or by possible freshwater release or removal by nearby manufacturing facilities. At all locations, the models often, but not always, underpredict when salinity is high. This can be seen in Fig. 6, which plots the relationship between the observed and modeled salinity. One concern was that this bias could have been caused by the use of smoothing splines to approximate the sharply nonlinear response to streamflow under low-flow conditions. However, experimental results described in the discussion indicated that the modeling methods were reasonable even when approximating exponentials. Furthermore, other smoothing methods did not perform better at correctly predicting high salinities. Long-term trends are also present in the model residuals (Fig. 7). This suggests that after accounting for streamflow, salinity in the estuary has been increasing with time. To test the significance of these trends, a parametric time term was added to Equations (1) and (2). The resulting trends and p-values are provided in Fig. 6. The difference between modeled and observed salinity versus modeled salinity. Mod line marks where there is no difference between observed and modeled values. The sign has higher than the observed salinity and negative values indicate that the model-predicted sa Table 2. Significant upward trends in streamflow-adjusted salinity are found at the oyster beds and Reedy Island. Upward trends are also found at all of the remaining locations; however, none of these trends are significantly different from zero at the 95% confidence level. The lack of statistical significance may reflect the short observational records, particularly at Ship John Shoal and Fort Mifflin, and the possibility of higher variability due to the proximity of the upstream locations (Chester, Fort Mifflin, and Ben Franklin Bridge) to the salt front. It should also be noted that each location covers slightly different time periods. 3.2. Effect of sea level Over the last century, sea level along the East Coast of the United States has risen significantly, and levels have been increasing faster than the global average along much of the coast (Sallenger et al., 2012). Locally, since measurements began in 1911, sea level at Atlantic City, NJ has risen at a rate of 0.41 m per century. The long- term trend has been nearly linear, although some short-term var- iations are present. Is sea-level rise responsible for the increasing trends in streamflow-adjusted salinity? This hypothesis was tested by including a parametric term for sea level in Equations (1) and (2) (and removing the time term previously included). Note that high and low tide water level, without sea-level rise or subtidal fluctuations, is already included in the oyster bed model to allow observations from both high and low tide to be combined. The resulting coefficients for the eled salinity values are from Equations (1) and (2) without the optional term. The black been configured so that positive values indicate that the model-predicted salinity was linity was less than the observed salinity. Fig. 7. Timeseries of model residuals. The gray dots are the raw residuals (observed minus modeled salinity) for the models in Equations (1) and (2) without the optional term. The black lines show the trends that result when a term for time is added to these models. Table 3 Responses of salinity to alongshore wind stress, along-estuary wind stress, and wind stress magnitude. p-values are indicated in parenthesis. Bold indicates locations where the results are significant at the 95% confidence level. Alongshore wind stress is defined as positive when it has a south-southwest to north-northeast component. Negative slopes indicate that salinity is lowered when the alongshore wind stress is from this direction. Along-estuary wind stress is defined as positive when it has a component pointing upstream perpendicular to the alongshore direction. Units for all values are psu per N m�2. Location Alongshore Along-estuary Magnitude Ben Franklin Bridge 7.1 � 10�2 (0.41) 0.11 (0.14) 0.11 (0.11) Fort Mifflin �7.3 � 10�2 (0.67) 0.65 (8.6 � 10�5) 0.21 (0.18) Chester �0.33 (0.31) 0.51 (9.2 � 10�2) 0.13 (0.55) Reedy Island ¡23 (2.0 � 10�4) 6.5 (0.25) 11.0 (4.8 � 10�2) Ship John Shoal ¡49 (1.2 � 10�2) �3.6 (0.72) 4.4 (0.74) A.C. Ross et al. / Estuarine, Coastal and Shelf Science 157 (2015) 79e9288 sensitivity of salinity to sea level are provided in Table 2. These results indicate that there is a significant correlation between sea level and estuarine salinity, with salinity increasing by as much as 4.4 per meter of sea-level rise. Due to the significant upward trends in streamflow-adjusted salinity already detected, any time series with a significant trend would likely be modeled as having a significant effect on salinity. However, the results here are reasonable, as there is a physical reason to expect increasing salinity with increasing sea level. Adjusted R2 values for the models with the time term were similar to those with the sea-level term (0.54 vs. 0.55, 0.60 vs. 0.61, 0.17 vs. 0.17, 0.75 vs. 0.75, 0.62 vs. 0.65, 0.84 vs. 0.83 in order of distance downstream). The model likelihoods (as well as Akaike information criterion values) were slightly better for all models with the sea- level term except at the oyster beds. Finally, several studies examined in the discussion using idealized and dynamical models have arrived at similar sensitivities. Idealized models with exponentially decreasing estuary widths predict that the response of salinity to sea level is not a simple linear function but rather a power function (Savenije, 1993; Hilton et al., 2008). A simple way of testing this non-parametrically is to replace the parametric relationship between salinity and sea level with a smooth term. Likelihood ratios indicated that the null model with a linear response to sea level was never rejected. This does not necessarily mean that the response of salinity to sea level is linear; however, it does indicate that any other form of response is not detected in the current data. 3.3. Effect of winds In addition to streamflow and sea level, wind stress may affect the salinity of the estuary by generating turbulent mixing and estuary-shelf exchanges. Because of the interaction between wind stress and sea level, the sea-level terms in the previous models were dropped and replaced with terms for wind stress. The indi- vidual wind stress components (alongshore, along-estuary, and magnitude) were also tested separately to avoid any problems caused by correlations between the components. Parametric terms were used to represent all of the wind stress terms because, like sea level, smooth relationships between salinity and wind stress were never significantly better. The results are presented in Table 3. Alongshore wind stress from the south-southwest to the north- northeast should induce Ekman transport away from the estuary and lower sea level and salinity. Wind stress from the opposite direction should have the opposite effect. The effect of low- frequency alongshore wind stress variability on subtidal sea level and circulation has been observed in the Delaware Estuary (Wong and Moses-Hall, 1998). Using generalized least squares (Pinheiro and Bates, 2000) with an AR(2) error covariance, there is also a significant negative correlation between sea level and alongshore wind stress in the data used in this study (p < 1 � 10�5). The resulting effects of wind stress on salinity were detected at Reedy Island and Ship John Shoal (Table 3). The values at the remaining locations were not significantly different from zero. Along-estuary wind stress may also affect salinity and sea level by directly inducing set-up inside the estuary. However, with the alongshore component in the previous regression models replaced by the along-estuary component, the effects of along-estuary wind stress were only significant at Fort Mifflin. This result is supported by Wong and Moses-Hall (1998), who found that although local wind effects have some influence on subtidal currents in the upper estuary, there is little coherence between subtidal current and surface salinity in this region. Similarly, Garvine (1985) determined that subtidal sea level and barotropic current fluctuations should primarily be produced by remote wind effects. The magnitude of wind stress may also have an effect on salinity, for example by increasing or decreasing vertical mixing in the es- tuary. The Delaware Estuary has strong tides and is traditionally A.C. Ross et al. / Estuarine, Coastal and Shelf Science 157 (2015) 79e92 89 considered weakly to partially stratified (Wong and Münchow, 1995; Wong, 1995), meaning any effect of wind mixing should be relatively small. However, strong vertical salinity stratification can occur, particularly during spring (Sharp et al., 1986). To test this influence, the along-estuary component in the previous regression model was replaced with a component for wind stress magnitude. A significant relationship was found at Reedy Island. The slopes at Reedy Island and the remaining insignificant locations were all positive, consistent with the hypothesis that higher winds will in- crease surface salinity by increasing vertical mixing. 3.4. Effect of oceanic salinity Finally, variability in oceanic salinity may also influence the salinity of the estuary. When a parametric term for the Gulf Stream Index is included in Equation (1), the term is significant only at Reedy Island (b1 ¼0.10,p ¼ 1.8 � 10�2). The sign of the slope b is positive, consistent with the hypothesis that a northward displacement of the Gulf Stream will increase oceanic salinity offshore on the Mid-Atlantic Bight and drive more saline water into the estuary. The terms are also positive at all of the remaining USGS locations except Chester; however, the terms are not statistically significant. This does not necessarily disprove the influence of oceanic salinity on estuarine salinity. Because the Gulf Stream Index is a highly variable time series, the statistical models may have difficulty identifying any relationship between oceanic and estua- rine salinity. In addition, there could be a lag between movement of the Gulf Stream and any associated changes in estuarine salinity, or the influence of oceanic salinity could occur over longer time scales that would not necessarily be detected by examining monthly averages. 4. Discussion 4.1. Sea level impact on estuarine salinity Assuming that the statistical models are reliable and the effects of unmodeled influences are negligible, the results show that a long-term upward trend in salinity is present after accounting for the effects of streamflow. This trend provides evidence that rising sea levels may be causing salinity to rise. If rising sea levels are in fact causing salinity to increase, salinity is likely to increase significantly in the future as a result of continued and accelerated sea-level rise. The Intergovernmental Panel on Climate Changes's (IPCC's) Fifth Assessment Report (AR5) predicts that global mean sea level will rise by 0.38e0.73 m during the twenty-first century under the RCP6.0 emissions scenario (Church et al., 2013). Other studies suggest that the AR5 may un- derestimate sea-level rise, with global mean sea level possibly rising 1 m or more above the 1990 mean by 2100 (Rahmstorf, 2007; Vermeer and Rahmstorf, 2009). Using statistical models, Vermeer and Rahmstorf (2009) project that global mean sea level will in- crease 1.24 m above the 1990 mean level by the end of the twenty- first century under the A2 emissions scenario (the models in this study range from 0.98 to 1.55 m). Using the modeled sensitivity to sea level of 3.3 m�1 at Reedy Island, this amounts to a 4.1 increase in mean salinity by the end of the twenty-first century with a range of 3.2e5.1. If streamflow is unchanged, this will raise the mean salinity at Reedy Island to 8.5 (7.6e9.5). Using the possibly conservative IPCC estimates of sea-level rise still results in an increase of 1.3e2.4, bringing the mean salinity at Reedy Island to 5.7e6.8. Statistically identifying sea-level rise as the cause of increasing salinity is difficult, since statistical models would correlate any time series with a large upward trend to the increasing salinity in the estuary. Other influences that are difficult to measure, such as dredging, could have also increased salinity in the estuary. How- ever, results from other models and studies of the effect of sea-level rise on salinity in the Mid-Atlantic region are generally similar to the results obtained from the statistical models in this study. Using numerical models, Hull and Tortoriello (1979) determined that a sea-level rise of 0.13 m resulted in a maximum increase in salinity of 0.38 near Reedy Island. This translates to a 2.9 m�1 sensitivity of salinity to sea level, which is slightly less than the 3.3 m�1 determined using the statistical models in this study. At a location 23 km upstream of Reedy Island, the U.S. Army Corps of Engineers Philadelphia District (1997) found that a 0.3-m increase in sea level resulted in a 0.3 increase in surface salinity, which translates to a sensitivity of 1.0 m�1. However, they also found that bottom salinity in the lower oyster bed area would actually decrease by 0.2 for a 0.3-m sea-level rise (a sensitivity of �0.7 m�1), which contradicts the 2.5 m�1 sensitivity determined by the sta- tistical models in this study. The Army Corps hypothesized that the negative sensitivity may be a result of flow diversions such as the C&D canal linking the Chesapeake and Delaware Bays; flow di- versions may have been introduced in the model by their approx- imation of sea-level rise in the Chesapeake. In the area upstream of Reedy Island, Kim and Johnson (2007) used numerical models to simulate the response of salinity under 1965 flow conditions and 1996 consumptive use (a worst-case scenario) to a 0.17-m sea-level rise. They found that chlorinity would increase by 0.14 ppt at Chester and 7 ppm at Ben Franklin Bridge. When converted to salinity in parts per thousand (which is nearly equivalent to the practical salinity unit used in this paper), this results in a sensitivity of 1.5 ppt m�1 at Chester and 7.1 � 10�2 ppt m�1 at Ben Franklin Bridge. This is much larger than the sensitivities identified by the statistical models in this study, which may be a result of the imposed low-flow conditions. Savenije (1993) developed a one-dimensional model for ideal- ized estuaries. The model assumes an estuary width and cross- sectional area that decrease exponentially upstream, assumes steady state at high water slack, and models dispersion using an empirical relation. Since the model is one-dimensional, the model salinity does not vary with depth. This may still be a reasonable first-order assumption, particularly since the Delaware Estuary is generally well-mixed. To compare this model with the statistical models in this study, we ran the calculations with the mean value of streamflow at Trenton during this study (349 m3 s�1), a mean estuary depth of 7.7 m approximated from the National Geophys- ical Data Center's Coastal Relief Model, and with the remaining model parameters as defined by Savenije (1993). The results indicate that the Savenije (1993) model-predicted salinity closely matches the observed salinity values (Fig. 8a). The response of salinity to 1 m of sea-level rise is similar to that predicted by the statistical models (Fig. 8b). Both models project the largest in- crease in the middle of the estuary and only minor changes farther upstream. Recent studies have also identified the implications of sea-level rise for salinity in other estuaries. In the nearby Chesapeake Bay, Hilton et al. (2008) applied statistical and numerical models and found that salinity has a sensitivity to sea level of 2e7 m�1, which is similar to the 2.5e4.4 m�1 sensitivity identified in this study. Similarly, using a 3D numerical model, Hong and Shen (2012) determined that the mean salinity of the Chesapeake Bay would increase by 1.2e2.0 if sea level rises by 1 m. With the same nu- merical model, Rice et al. (2012) found that, during typical flow conditions, a 1-m sea-level rise would increase salinity by nearly 10 in the James River at the mouth of the Chickahominy River. During a simulated dry year, a 1-m sea-level rise would result in a salinity increase of slightly more than 4, which is more in line with the estimates by Hilton et al. (2008). Fig. 8. a: Comparison between observed salinity and salinity predicted by the Savenije (1993) 1D model. The shapes are the observed salinity, and the lines denote the 1D model predictions. The distance for the oyster beds is determined by the mean dis- tance weighted by number of observations. b: Projections of salinity change in response to 1 m of sea-level rise under current mean streamflow. Shapes are from the statistical models developed in this work. Solid lines are from the Savenije (1993) 1D model. c: Projections of salinity change in response to a 35% increase or decrease in streamflow under current mean sea levels. As before, solid lines are from the Savenije (1993) 1D model, and shapes are from the statistical models developed in this work. A.C. Ross et al. / Estuarine, Coastal and Shelf Science 157 (2015) 79e9290 4.2. Delaware Estuary salinity response to streamflow Changes in streamflow are more difficult to project. In the Eastern United States, precipitation primarily determines streamflow. Many studies have projected that precipitation amounts will increase, both globally (Collins et al., 2013) and over the Eastern United States (Najjar et al., 2009). However, land use changes, reduced plant transpiration as a result of increased CO2, and increased evaporation as a result of higher temperatures (Krakauer and Fung, 2008) will also influence streamflow change. Using global climate models, Najjar et al. (2009) projected that precipitation changes will cause streamflows in the Mid-Atlantic region to change by 15 ± 20% by the end of the twenty-first cen- tury under the A2 emissions scenario. However, Najjar et al. (2009) also found that warming-induced evapotranspiration changes could cause a 15e40% decrease in streamflows in the region. If sea level is held constant at the 1964e2012 mean, the statistical model for Reedy Island predicts that salinity would increase from the 1964e2012 mean value of 4.4e5.5 for a 35% decrease in mean streamflow or decrease to 3.7 for a 35% increase in mean streamflow. Thus, only an extreme streamflow increase will be able to offset the salinity change caused by a moderate 1 m sea-level rise. The 1D model of Savenije (1993) also predicts a similar response of salinity to streamflow change (Fig. 8c). As a result of the quasi- exponential shape of the salinity-streamflow curve, salinity changes more in response to an increase of streamflow than it does for an equal decrease of streamflow. The 1D model predicts a larger response to streamflow change than the statistical models predict at Reedy Island and downstream; however, 1D models are known to overpredict the sensitivity of salinity to streamflow in the Delaware Estuary (Garvine et al., 1992). 4.3. Limitations of statistical models Although the statistical models generally performed well and produced results that make sense from a physical standpoint and are reasonably close to other studies, the models have some shortcomings. One issue is that the models often underpredicted extremely high salinities. It was thought that this issue may have reflected an inability of the smoothing splines to fit the roughly exponential salinity-streamflow relationship. This hypothesis was tested in three steps. First, simulations were conducted in which exponential curves were created from the observed salinity- streamflow relationship. For example, a curve was created with shape S ¼ a � exp(Q/b), where the coefficients a and b were deter- mined using nonlinear least squares and observed S and Q values. Next, various levels of random noise (including no noise) were added to these exponential curves. Finally, the ability of the AMM smoothing splines to fit the exponential curves with noise was tested. In all cases, the splines were remarkably close to the actual exponentials. In addition, using other smoothing splines, increasing the maximum degrees of freedom in the splines, and manually setting the spline knots did not improve the model fits. Specifying exponential or power law relationships between salinity and streamflow, rather than using smoothing functions, did not improve the model results either. Since it otherwise appears to work well, fitting splines to the raw quasi-exponential relationship is advantageous over other methods for fitting quasi-exponential curves, such as applying a log transform to salinity, because it preserves the additive nature of the model. Because there are seasonal patterns in many of the variables, there is also the potential for some concurvity issues when including multiple variables as predictors in the model (concurvity, or approximate concurvity, refers to the presence of nonlinear re- lationships between predictor variables (Buja et al., 1989; Ramsay et al., 2003)). This issue may also arise from the time-of-year term that accounts for residual seasonal variation. However, the month term significantly improves the model fits, and excluding it from the models results in seasonal patterns in the residuals. The cause of this residual seasonal variation is unknown. One possibility is seasonal patterns in evaporation and precipitation. In the reanalysis data used to obtain wind speed and direction, the evaporation rate is generally lowest in AprileJune (roughly 10 mm per month) and highest in OctobereDecember (roughly 100 mm per month). The seasonal pattern in precipitation is fairly small, so the difference between precipitation and evaporation follows a similar pattern, peaking at a net evaporation of 56 mm per month in October and a net precipitation of 42 mm per month in March. A.C. Ross et al. / Estuarine, Coastal and Shelf Science 157 (2015) 79e92 91 This seasonal pattern in E � P should cause decreasing salinity in MarcheJune and increasing salinity in OctobereDecember. Indeed, this very closely matches the seasonal pattern in the model re- siduals (Fig. 5). However, the actual impact of evaporation is likely negligible. As a simple example, consider a 1 m2 by 1-m deep volume of water at the surface. After one month of evaporation and precipitation, the salinity of this water is Sf ¼ Si � 1 þ mEPmi�mEP � where Si and mi are the salinity and total mass of water at the start of the month respec- tively and mEP is the E � P accumulated mass loss during the month. Applying this equation to each month of the NARR annual cycle of E � P results in an approximate annual cycle of salinity resulting solely from evaporation and precipitation. Using Reedy Island as an example, with a starting January salinity of 4.0, results in an annual cycle with a range of 0.66 compared to the range of the residual cycle of roughly 3. Furthermore, this example assumes that there is no mixing below 1 m throughout the year, which is extremely unrealistic. Using the approximate mean depth of the estuary (7.7 m) results in an annual cycle range of only 9.4 � 10�2 at Reedy Island. Other factors that are difficult to model, such as changes in the width and depth of the Delaware River navigation channel, may also influence the salinity of the Delaware Estuary. However, the width and depth of the navigation channel remained relatively stable for the majority of the time period of this study. The project to deepen the channel downstream of Philadelphia to a depth of 40 feet (approximately 12.2 m) was completed in 1942 (U.S. Army Corps of Engineers Philadelphia District (2009)), and the autho- rized depth remained at 40 feet until 1992. Walsh (2004) found minor changes after the channel deepening, with primarily small- scale dredging projects such as an extension of the Marcus Hook anchorage in 1964 and minor accretion in the area surrounding the channel from 1980e1987 to 2001. Although the authorized channel depth was increased to 45 feet (approximately 13.7 m) in 1992, work to deepen the channel to this depth did not begin until 2010, and as of 2014 this work has not been completed. Numerical model simulations performed by Kim and Johnson (2007) predicted that the increase in channel depth from 40 feet to 45 feet would cause a 6.33% increase at Chester under 1965 drought conditions. The sta- tistical models in this study did not detect a significant trend in salinity at Chester, suggesting that the effect of the recent deep- ening has not influenced long-term salinity trends. By 2040, Kim and Johnson (2007) projects that sea-level rise will have increased salinity more than channel deepening and increased water consumption combined. 5. Conclusion After accounting for the effects of streamflow and seasonal variations, salinity in many areas of the Delaware Estuary is increasing. This increase may be caused by sea-level rise. If the future response of salinity to sea level matches the modeled past response, salinity will increase significantly in the future as sea level continues to rise. Any increase in streamflow caused by warming will likely be unable to balance the increase in salinity caused by sea-level rise. Although the statistical models used in this study appear to have worked well, additional investigation into the ability of the methods to handle the salinity-streamflow relationship may be beneficial. In addition, the comparison of the statistical model re- sults with results from numerical models would benefit from modern numerical model simulations forced with the full range of possible streamflow conditions (rather than only low-flow conditions). Acknowledgments Mike Loewen provided assistance in reading the oyster bed salinity data from punch cards. We thank the three anonymous re- viewers whose comments improved this paper. Support for this research was provided by the National Science Foundation Physical Oceanography Program (award #0961423) and Pennsylvania Sea Grant (4796-TPSU-NOAA-0061). North American Regional Rean- alysis data was provided by the NOAA/OAR/ESRL PSD, Boulder, Col- orado, USA, from their Web site at http://www.esrl.noaa.gov/psd/. References Aristiz�abal, M., Chant, R., 2013. A numerical study of salt fluxes in Delaware Bay estuary. J. Phys. Oceanogr. 43, 1572e1588. Autin, M., Edwards, D., 2010. Nonparametric harmonic regression for estuarine water quality data. Environmetrics 21, 588e605. Banas, N., Hickey, B., MacCready, P., Newton, J., 2004. Dynamics of Willapa Bay, Washington: a highly unsteady, partially mixed estuary. J. Phys. Oceanogr. 34, 2413e2427. Bernhard, A.E., Donn, T., Giblin, A.E., Stahl, D.A., 2005. Loss of diversity of ammonia- oxidizing bacteria correlates with increasing salinity in an estuary system. Environ. Microbiol. 7, 1289e1297. Breunig, M., Kriegel, H.-P., Ng, R.T., Sander, J., 2000. LOF: identifying density-based local outliers. In: Proceedings of the 2000 ACM SIGMOD International Confer- ence on Management of Data, pp. 93e104. Buja, A., Hastie, T., Tibshirani, R., 1989. Linear smoothers and additive models. Ann. Statistics 17, 453e510. Bushek, D., Ford, S.E., Burt, I., 2012. Long-term patterns of an estuarine pathogen along a salinity gradient. J. Mar. Res. 70, 225e251. Church, J., Clark, P., Cazenave, A., Gregory, J., Jevrejeva, S., Levermann, A., Merrifield, M., Milne, G., Nerem, R., Nunn, P., Payne, A., Pfeffer, W., Stammer, D., Unnikrishnan, A., 2013. Sea level change. In: Stocker, T., Qin, D., Plattner, G.-K., Tignor, M., Allen, S., Boschung, J., Nauels, A., Xia, Y., Bex, V., Midgley, P. (Eds.), Climate Change 2013: the Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change. Cambridge University Press, New York and Cambridge. Cloern, J.E., Knowles, N., Brown, L.R., Cayan, D., Dettinger, M.D., Morgan, T.L., Schoellhamer, D.H., Stacey, M.T., van der Wegen, M., Wagner, R.W., Jassby, A.D., 2011. Projected evolution of California's San Francisco Bay-Delta-River system in a century of climate change. PLoS ONE 6. Collins, M., Knutti, R., Arblaster, J., Dufresne, J.L., Fichefet, T., Friedlingstein, P., Gao, X., Gutowski, W.J., Johns, T., Krinner, G., Shongwe, M., Tebaldi, C., Weaver, A.J., Wehner, M., 2013. Long-term climate change: projections, com- mitments and irreversibility. In: Stocker, T., Qin, D., Plattner, G.-K., Tignor, M., Allen, S., Boschung, J., Nauels, A., Xia, Y., Bex, V., Midgley, P. (Eds.), Climate Change 2013: the Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change. Cambridge University Press, New York and Cambridge. Corbeil, R.R., Searle, S.R., 1976. Restricted maximum likelihood (REML) estimation of variance components in the mixed model. Technometrics 18, 31e38. Coull, B.A., Schwartz, J., Wand, M., 2001. Respiratory health and air pollution: ad- ditive mixed model analyses. Biostatistics 2, 337e349. Cronin, T., Sanford, L., Langland, M., Willard, D., Saenger, C., 2003. Estuarine sedi- ment transport, deposition, and sedimentation. In: Langland, M., Cronin, T. (Eds.), A Summary Report of Sediment Processes in Chesapeake Bay and Watershed. U.S. Geological Survey, pp. 61e79. Ch. 6. Ford, S.E., Scarpa, E., Bushek, D., 2012. Spatial and temporal variability of disease refuges in an estuary: implications for the development of resistance. J. Mar. Res. 70, 253e277. Gallegos, C., Jordan, T., 2002. Impact of the Spring 2000 phytoplankton bloom in Chesapeake Bay on optical properties and light penetration in the Rhode River, Maryland. Estuaries Coasts 25, 508e518. Garvine, R.W., 1985. A simple model of estuarine subtidal fluctuations forced by local and remote wind stress. J. Geophys. Res. Oceans 90, 11945e11948. Garvine, R.W., McCarthy, R.K., Wong, K.-C., 1992. The axial salinity distribution in the Delaware Estuary and its weak response to river discharge. Estuar. Coast. Shelf Sci. 35, 157e165. Gibson, J., Najjar, R.G., 2000. The response of Chesapeake Bay salinity to climate- induced changes in streamflow. Limnol. Oceanogr. 45, 1764e1772. Gordon, A.B., Katzenbach, M., 1983. Guidelines for Use of Water Quality Monitors. Tech. Rep. Open-File Report 83e681. U.S. Geological Survey. Hartmann, D., Klein Tank, A., Rusticucci, M., Alexander, L., Bronnimann, S., Charabi, Y., Dentener, F., Dlugokencky, E., Easterling, D., Kaplan, A., Soden, B., Thorne, P., Wild, M., Zhai, P., 2013. Observations: atmosphere and surface. In: Stocker, T., Qin, D., Plattner, G.-K., Tignor, M., Allen, S., Boschung, J., Nauels, A., Xia, Y., Bex, V., Midgley, P. (Eds.), Climate Change 2013: the Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change. Cambridge University Press, New York and Cambridge. http://www.esrl.noaa.gov/psd/ http://refhub.elsevier.com/S0272-7714(15)00068-2/sref1 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref1 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref1 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref1 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref2 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref2 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref2 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref3 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref3 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref3 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref3 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref4 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref4 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref4 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref4 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref5 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref5 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref5 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref5 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref6 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref6 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref6 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref7 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref7 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref7 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref8 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref8 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref8 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref8 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref8 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref8 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref8 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref9 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref9 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref9 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref9 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref10 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref10 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref10 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref10 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref10 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref10 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref10 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref10 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref11 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref11 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref11 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref12 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref12 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref12 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref13 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref13 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref13 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref13 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref13 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref14 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref14 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref14 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref14 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref15 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref15 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref15 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref15 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref16 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref16 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref16 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref17 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref17 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref17 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref17 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref18 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref18 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref18 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref19 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref19 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref19 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref20 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref20 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref20 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref20 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref20 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref20 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref20 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref20 A.C. Ross et al. / Estuarine, Coastal and Shelf Science 157 (2015) 79e9292 Haskin, H., 1972. Delaware River Flow-bay Salinity Relationships. Tech. rep. Dela- ware River Basin Commission. Haskin, H.H., Ford, S.E., 1982. Haplosporidium nelsoni (MSX) on Delaware Bay seed oyster beds: a host-parasite relationship along a salinity gradient. J. Invertebr. Pathol. 40, 388e405. Hilton, T.W., Najjar, R.G., Zhong, L., Li, M., 2008. Is there a signal of sea-level rise in Chesapeake Bay salinity? J. Geophys. Res. 113. Hirsch, R.M., Slack, J.R., 1984. A nonparametric trend test for seasonal data with serial dependence. Water Resour. Res. 20, 727e732. Holgate, S.J., Matthews, A., Woodworth, P.L., Rickards, L.J., Tamisiea, M.E., Bradshaw, E., Foden, P.R., Gordon, K.M., Jevrejeva, S., Pugh, J., 2013. New data systems and products at the permanent service for mean sea level. J. Coast. Res. 29, 493e504. Hong, B., Shen, J., 2012. Responses of estuarine salinity and transport processes to potential future sea-level rise in the Chesapeake Bay. Estuar. Coast. Shelf Sci. 104e105, 33e45. Hull, C.H., Titus, J., 1986. Greenhouse Effect, Sea Level Rise, and Salinity in the Delaware Estuary. United States Environmental Protection Agency and the Delaware River Basin Commission. Hull, C.H., Tortoriello, R.C., 1979. Sea-level Trend and Salinity in the Delaware Es- tuary. Tech. rep.. Delaware River Basin Commission. Jiang, Y., Luo, Y., Zhao, Z., Tao, S., 2009. Changes in wind speed over China during 1956e2004. Theor. Appl. Climatol. 99, 421e430. Jolly, I.D., Williamson, D.R., Gilfedder, M., Walker, G.R., Morton, R., Robinson, G., Jones, H., Zhang, L., Dowling, T.I., Dyce, P., Nathan, R.J., Nandakumar, N., Clarke, R., McNeill, V., 2001. Historical stream salinity trends and catchment salt balances in the Murray-Darling Basin, Australia. Mar. Freshw. Res. 52, 53e63. Katzenbach, M.S., 1990. Comparison of Accuracy and Completeness of Data Ob- tained from Three Types of Automatic Water-quality Monitors. Tech. Rep. Water-Resources Investigations Report 89-4198. U.S. Geological Survey. Kauffman, G., Homsey, A., Belden, A., Sanchez, J., 2011. Water quality trends in the Delaware River Basin (USA) from 1980 to 2005. Environ. Monit. Assess. 177, 193e225. Kennish, M., 2002. Environmental threats and environmental future of estuaries. Environ. Conserv. 29, 78e107. Ketchum, B., 1952. Circulation in estuaries. In: Proceedings of Third Conference on Coastal Engineering, pp. 65e76. Kim, K.W., Johnson, B.H., 2007. Salinity Re-validation of the Delaware Bay and River 3-D Hydrodynamic Model with Applications to Assess the Impact of Channel Deepening, Consumptive Water Use, and Sea Level Change. Tech. rep.. U.S. Army Research and Development Center, Vicksburg, MS. Krakauer, N., Fung, I., 2008. Mapping and attribution of change in streamflow in the coterminous United States. Hydrology Earth Syst. Sci. 12, 1111e1120. Lee, Y.J., Lwiza, K.M., 2008. Factors driving bottom salinity variability in the Ches- apeake Bay. Cont. Shelf Res. 28, 1352e1362. Lewis, E., Perkin, R., 1981. The practical salinity scale 1978: conversion of existing data. Deep Sea Res. Part A. Oceanogr. Res. Pap. 28, 307e328. Lin, X., Zhang, D., 1999. Inference in generalized additive mixed models by using smoothing splines. J. R. Stat. Soc. Ser. B Stat. Methodol. 61, 381e400. Marshall, F., Smith, D., Nickerson, D., 2011. Empirical tools for simulating salinity in the estuaries in Everglades National Park, Florida. Estuar. Coast. Shelf Sci. 95, 377e387. Mesinger, F., DiMego, G., Kalnay, E., Mitchell, K., Shafran, P., Ebisuzaki, W., Jovic, D., Woollen, J., Rogers, E., Berbery, E., et al., 2006. North American Regional Reanalysis. Bull. Am. Meteorol. Soc. 87, 343e360. Morton, R., Henderson, B.L., 2008. Estimation of nonlinear trends in water quality: an improved approach using generalized additive models. Water Resour. Res. 44. Najjar, R., Patterson, L., Graham, S., 2009. Climate simulations of major estuarine watersheds in the Mid-Atlantic region of the US. Clim. Change 95, 139e168. Pinheiro, J., Bates, D., 2000. Mixed Effects Models in S and S-plus. In: Statistics and Computing Series. Springer-Verlag, New York. Powell, E.N., Gauthier, J.D., Wilson, E.A., Nelson, A., Fay, R.R., Brooks, J.M., 1992. Oyster disease and climate change. Are yearly changes in Perkinsus marinus parasitism in oysters (Crassostrea virginica) controlled by climatic cycles in the Gulf of Mexico? Mar. Ecol. 13, 243e270. Rahmstorf, S., 2007. A semi-empirical approach to projecting future sea-level rise. Science 315, 368e370. Ramsay, T.O., Burnett, R.T., Krewski, D., 2003. The effect of concurvity in generalized additive models linking mortality to ambient particulate matter. Epidemiology 14, 18e23. Rice, K.C., Hong, B., Shen, J., 2012. Assessment of salinity intrusion in the James and Chickahominy Rivers as a result of simulated sea-level rise in Chesapeake Bay, East Coast, USA. J. Environ. Manag. 111, 61e69. Saenger, C., Cronin, T., Thunell, R., Vann, C., 2006. Modelling river discharge and precipitation from estuarine salinity in the northern Chesapeake Bay: appli- cation to Holocene palaeoclimate. Holocene 16, 467e477. Sallenger, A.H., Doran, K.S., Howd, P.A., 2012. Hotspot of accelerated sea-level rise on the Atlantic coast of North America. Nat. Clim. Change 2, 1e5. Sanchez, J., Kauffman, G., Reavy, K., Homsey, A., 2012. Watersheds & landscapes. In: Technical Report for the Delaware Estuary & Basin. Partnership for the Dela- ware Estuary, pp. 14e47. Ch. 1. Savenije, H.H., 1993. Predictive model for salt intrusion in estuaries. J. Hydrol. 148, 203e218. Schemel, L., 2001. Simplified Conversion Between Specific Conductance and Salinity Units for Use with Data from Monitoring Stations. Sharp, J.H. (Ed.), 1983. The Delaware Estuary: Research as Background for Estuarine Management and Development. Sharp, J.H., Cifuentes, L.A., Coffin, R.B., Pennock, J.R., Wong, K.-C., 1986. The influence of river variability on the circulation, chemistry, and microbiology of the Delaware Estuary. Estuaries 9, 261e269. Simpson, G.L., Anderson, N., 2009. Deciphering the effect of climate change and separating the influence of confounding factors in sediment core records using additive models. Limnol. Oceanogr. 54, 2529e2541. Taylor, A., 1995. North-South shifts of the Gulf Stream and their climatic connection with the abundance of zooplankton in the UK and its surrounding seas. ICES J. Mar. Sci. 52, 711e721. U.S. Army Corps of Engineers Philadelphia District, 1997. Chapter 5: hydrodynamic and salinity modeling. In: Delaware River Main Channel Deepening Project: Supplemental Environmental Impact Statement, pp. 5-1e5-62. U.S. Army Corps of Engineers Philadelphia District, 2009. Delaware River Main Stem and Channel Deepening Project, Environmental Assessment. Tech. rep.. U.S. Army Corps of Engineers. Vautard, R., Cattiaux, J., Yiou, P., Thepaut, J.-N., Ciais, P., 2010. Northern Hemisphere atmospheric stilling partly attributed to an increase in surface roughness. Nat. Geosci. 3, 756e761. Vermeer, M., Rahmstorf, S., 2009. Global sea level linked to global temperature. Proc. Natl. Acad. Sci. U. S. A. 106. Walsh, D.R., 2004. Anthropogenic Influences on the Morphology of the Tidal Delaware River and Estuary: 1877e1987 (Master's thesis). University of Delaware. Wilks, S.S., 1938. The large-sample distribution of the likelihood ratio for testing composite hypotheses. Ann. Math. Stat. 9. Wong, K., Moses-Hall, J., 1998. On the relative importance of the remote and local wind effects to the subtidal variability in a coastal plain estuary. J. Geophys. Res. 103, 393e404. Wong, K.-C., 1994. On the nature of transverse variability in a coastal plain estuary. J. Geophys. Res. Oceans 99, 14209e14222. Wong, K.-C., 1995. On the relationship between long-term salinity variations and river discharge in the middle reach of the Delaware estuary. J. Geophys. Res. 100, 20705e20713. Wong, K.-C., Garvine, R.W., 1984. Observations of wind-induced, subtidal variability in the Delaware Estuary. J. Geophys. Res. Oceans 89, 10589e10597. Wong, K.-C., Münchow, A., 1995. Buoyancy forced interaction between estuary and inner shelf: observation. Cont. Shelf Res. 15, 59e88. Wood, S.N., 2003. Thin plate regression splines. J. R. Stat. Soc. Ser. B Stat. Methodol. 65, 95e114. Wood, S.N., 2004. Stable and efficient multiple smoothing parameter estimation for generalized additive models. J. Am. Stat. Assoc. 99, 673e686. Wood, S.N., 2006. Generalized Additive Models: an Introduction with R. Chapman and Hall/CRC, Boca Raton. Wood, S.N., 2014. Package 'mgcv' (accessed 17.11.14.). URL. http://cran.r-project.org/ web/packages/mgcv/mgcv.pdf. Wu, J., 1982. Wind-stress coefficients over sea surface from breeze to hurricane. J. Geophys. Res. 87, 9704e9706. http://refhub.elsevier.com/S0272-7714(15)00068-2/sref21 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref21 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref22 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref22 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref22 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref22 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref23 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref23 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref24 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref24 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref24 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref25 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref25 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref25 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref25 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref25 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref26 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref26 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref26 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref26 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref26 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref27 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref27 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref27 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref28 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref28 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref29 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref29 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref29 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref29 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref30 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref30 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref30 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref30 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref30 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref31 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref31 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref31 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref32 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref32 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref32 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref32 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref33 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref33 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref33 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref34 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref34 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref34 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref35 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref35 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref35 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref35 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref36 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref36 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref36 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref37 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref37 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref37 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref38 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref38 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref38 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref39 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref39 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref39 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref40 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref40 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref40 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref40 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref41 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref41 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref41 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref41 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref42 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref42 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref42 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref43 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref43 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref43 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref44 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref44 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref45 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref45 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref45 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref45 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref45 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref46 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref46 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref46 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref47 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref47 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref47 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref47 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref48 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref48 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref48 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref48 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref49 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref49 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref49 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref49 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref50 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref50 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref50 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref51 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref51 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref51 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref51 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref51 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref51 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref52 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref52 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref52 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref53 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref53 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref54 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref54 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref55 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref55 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref55 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref55 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref56 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref56 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref56 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref56 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref57 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref57 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref57 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref57 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref58 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref58 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref58 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref58 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref59 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref59 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref59 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref60 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref60 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref60 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref60 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref61 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref61 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref62 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref62 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref62 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref62 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref63 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref63 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref64 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref64 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref64 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref64 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref65 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref65 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref65 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref66 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref66 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref66 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref66 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref67 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref67 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref67 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref68 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref68 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref68 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref69 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref69 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref69 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref70 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref70 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref70 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref71 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref71 http://cran.r-project.org/web/packages/mgcv/mgcv.pdf http://cran.r-project.org/web/packages/mgcv/mgcv.pdf http://refhub.elsevier.com/S0272-7714(15)00068-2/sref73 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref73 http://refhub.elsevier.com/S0272-7714(15)00068-2/sref73 Sea-level rise and other influences on decadal-scale salinity variability in a coastal plain estuary 1. Introduction 2. Methods 2.1. Study area and data 2.2. Statistical models 2.2.1. Smoothing functions 2.2.2. Distributions 2.2.3. Fitting and testing 3. Results 3.1. Basic models 3.2. Effect of sea level 3.3. Effect of winds 3.4. Effect of oceanic salinity 4. Discussion 4.1. Sea level impact on estuarine salinity 4.2. Delaware Estuary salinity response to streamflow 4.3. Limitations of statistical models 5. Conclusion Acknowledgments References