Spatial distribution of water level impacting back-barrier bays

Water level in semi-enclosed bays, landward of barrier islands, is mainly driven by offshore sea level fluctuations that are modulated by bay geometry and bathymetry, causing spatial variability in the ensuing response (transfer). Local wind setup can have a complementary role that depends on wind speed, fetch, and relative orientation of the wind direction and the bay. Bay area and inlet geometry and bathymetry primarily regulate the magnitude of the transfer between open ocean and bay. Tides and short-period offshore oscillations are more damped in the bays than longer-lasting offshore fluctuations, such as a storm surge and sea level rise. We compare observed and modeled water levels at stations in a mid-Atlantic bay (Barnegat Bay) with offshore water level proxies. Observed water levels in Barnegat Bay are compared and combined with model results from the Coupled Ocean–Atmosphere–Wave–Sediment Transport (COAWST) modeling system to evaluate the spatial structure of the water level transfer. Analytical models based on the dimensional characteristics of the bay are used to combine the observed data and the numerical model results in a physically consistent approach. Model water level transfers match observed values at locations inside the bay in the storm frequency band (transfers ranging from 50 %–100 %) and tidal frequencies (10 %–55 %). The contribution of frequency-dependent local setup caused by wind acting along the bay is also considered. The wind setup effect can be comparable in magnitude to the offshore transfer forcing during intense storms. The approach provides transfer estimates for locations inside the bay where observations were not available, resulting in a complete spatial characterization. An extension of the methodology that takes advantage of the ADCIRC tidal database for the east coast of the United States allows for the expansion of the approach to other bay systems. Detailed spatial estimates of water level transfer can inform decisions on inlet management and contribute to the assessment of current and future flooding hazard in back-barrier bays and along mainland shorelines.

1 Introduction

Back-barrier bays or coastal lagoons are common features along the coast of the United States. Their depths are usually on the order of a few meters and their horizontal extents are on the order of several tens of kilometers. They are often surrounded by highly populated areas and susceptible to intense human and environmental stressors. During storms, a surge and larger-than-normal waves combine to inundate low-elevation areas, resulting in hazards to coastal communities. Both hurricanes and winter storms affect coastal populations, infrastructure, and natural resources along the coastal bays of the United States (Nicholls et al., 2007, 2014; Rahmstorf, 2017; Wahl et al., 2017).

Hazard assessments consist of a characterization of the spatial and temporal extent of damaging physical events and the determination of the specific characteristics of those events (Ludwig et al., 2018). While flooding on the mainland side of back-barrier bays has severe socioeconomic implications, most coastal hazard evaluations (Gornitz et al., 1994; Thieler and Hammar-Klose, 1999; Klein and Nicholls, 1999; Kunreuther et al., 2000; Neumann et al., 2015; Vitousek et al., 2017) have focused on open-coast areas. Vulnerability evaluation of coastal areas around back-barrier bays requires extensive knowledge of the main hazard sources and their physical controls.

Water level in the bays is mainly driven by offshore sea level fluctuations with additional effects from local wind and wave setups. The bay exchange with the ocean usually occurs through narrow inlets. The size of the inlet determines the frictional effects and the amount of dampening offshore fluctuations encounter (Keulegan, 1967). Tides and short-period offshore oscillations tend to be more dampened in the bays than longer-lasting offshore fluctuations, such as a storm surge and sea level rise.

Bay water level fluctuations are linked to offshore forcing, especially at low frequencies, while wind acting directly over the bay is more connected to current fluctuations in the bay (Garvine, 1985). Chuang and Swenson (1981) determined that water level changes at subtidal frequencies in Lake Pontchartrain were controlled by coupled coastal ocean–bay fluctuations. Wong and Wilson (1984) studied subtidal sea level fluctuations in Great South Bay and again found them primarily driven by bay–shelf coupling. In Delaware Bay, a bay–inlet system with a relatively large opening, Wong and DiLorenzo (1988) showed that remote effects dominate over local effects and that fluctuations at both tidal and subtidal frequencies in connected bays of the Delaware Bay system were forced by shelf sea level.

More recently, Aretxabaleta et al. (2014) analyzed water level data in Barnegat Bay and Great South Bay before and after Hurricane Sandy and demonstrated that the offshore–bay transfer was not significantly altered by the geomorphologic changes caused by the storm. Aretxabaleta et al. (2017) described observed changes in both tidal amplitude and bay water level transfer from offshore in Great South Bay and connected bays and related the changes to the dredging of nearby inlets and the changing size of a breach across Fire Island caused by Hurricane Sandy. They also introduced an analytical model, based on the Chuang and Swenson (1981) approach but extended to interconnected bays, that incorporated bay and inlet dimensions and matched the observed transfer of offshore water level fluctuations into the bay system.

In this study, we combine an analysis of observed water levels in Barnegat Bay with the results of numerical models and an analytical description of the system to characterize the spatial characteristics of the bay response to offshore fluctuations. The observations provide detailed information at five locations in the bay, while the numerical simulations can expand the analysis to the entire bay system. The analytical model allows for the evaluation of the importance of the dominant factors affecting water level in bays. The combined approach can be used to provide consistent spatial maps of offshore water level impact into back-barrier bays. The method will be useful for coastal hazard assessment, assisting in the management of nuisance flooding (Moftakhari et al., 2018), providing spatial differences in vulnerability to perigean spring tides (king tides), and planning for flooding in response to storms of different durations. The final hazard estimates will be included as part of the U.S. Geological Survey Coastal Change Hazard portal (USGS, 2018) in an effort to expand the total water level predictions (Aretxabaleta et al., 2019).

2 Regional description

The Barnegat Bay–Little Egg Harbor (BBLEH) estuary is a back-barrier bay along the coast of New Jersey (Fig. 1). It is a shallow (average depth around 1.5 m) bay connected to the ocean through three openings: Little Egg Inlet in the south, Barnegat Inlet in the center, and the Point Pleasant Canal, which is a much smaller connection in the north of the bay. Offshore tidal amplitudes decrease slightly from 0.7 m in northern New York Bay to 0.6 m in central New Jersey. The southern sub-embayment (Little Egg Harbor) is more connected to the open ocean, with tidal amplitudes ranging between 0.2 and 0.5 m, while its northern part (Barnegat Bay) has less exchange and tidal amplitudes are smaller than 0.2 m (Chant, 2001; Defne and Ganju, 2015).

Figure 1 Map of Barnegat Bay and Little Egg Harbor estuary showing the water level stations, bays, and inlets. The water level stations are Tuckerton (TUC), East Thorofare (ETH), Waretown (WAR), Seaside Heights (SEH), and Mantoloking (MAN). Locations of offshore water level proxy stations and wind buoys are indicated in inset. The COAWST model domain boundary is shown in red. Route 72 crosses the bay near ETH station and the breach that occurred during Hurricane Sandy was about 100 m away from MAN station, so they are not indicated in the map.

3 Observational and model data

Water level observations from five stations in the BBLEH system (Table 1) and from two external coastal stations are used to determine transfer from ocean to bay. The bay stations started recording in October 2007, while Sandy Hook and Atlantic City are long-term NOAA water level stations, operational since 1910 and 1911, respectively (Table 1). Wind observations were obtained from the National Data Buoy Center (NDBC) buoy 44065 (New York Harbor Entrance) for the period 2008–2018.

Table 1 Sites used in water level analysis. Check Fig. 1 for locations. Information on instrumentation type, sampling, and quality control methodologies for the USGS stations is available from https://waterdata.usgs.gov/nwis/inventory (last access: 15 August 2019). National Oceanic and Atmospheric Administration (NOAA); U.S. Geological Survey (USGS). MSL: mean sea level; NGVD29: National Geodetic Vertical Datum of 1929; NADV88: North American Vertical Datum of 1988.

We used numerical simulations of Barnegat Bay for the period March–September 2012 (Defne and Ganju, 2015; data available from Defne and Ganju, 2018) and October–December 2012 (USGS, 2019) to obtain the spatial character of the water level response. The simulations used the Coupled Ocean–Atmosphere–Wave–Sediment Transport modeling system (COAWST; Warner et al., 2010). The model resolution ranged from 40 to 200 m, with the higher resolution located near complex geometry and around the inlets. The model is forced at the boundaries with tides from the ADCIRC tidal database for the western North Atlantic Ocean (Szpilka et al., 2016) and open-ocean forcing from subtidal water level and velocity from the ESPreSSO model (Wilkin and Hunter, 2013; http://www.myroms.org/espresso/, last access: 15 August 2019) and COAWST US east coast forecast (Warner et al., 2010; https://catalog.data.gov/dataset/coawst-forecast-system-usgs-us-east-coast-and-gulf-of, last access: 15 August 2019). Defne and Ganju (2015) showed that the numerical model solution had sufficient flow and elevation skill to characterize bay dynamics under normal and storm conditions.

The impact in the bay of offshore forcing can be evaluated spectrally by estimating transfer functions in frequency space between observed water levels offshore (input) and in the bays (output). The transfer functions are calculated using a Hanning window with overlapping (50 %) data segments with lengths of 29 d to provide estimates near the main tidal frequencies (Aretxabaleta et al., 2017). The uncertainty envelopes for the transfer function are estimated using the Bendat and Piersol (1986) formulation.

4 Analytical water level models

4.1 Offshore impact on bay model

The impact of ocean sea level fluctuations in the bay can be explored with an analytical model of a generic bay system (Fig. 2), consisting of multiple interconnected sub-embayments connected to the offshore by three separate inlets: Little Egg Inlet, Barnegat Inlet, and the Point Pleasant Canal. The model assumes that the bay water level responds as a level surface in each sub-embayment to ocean fluctuations, as local forcing in the bay is not included. The formulation is an extension of the approach proposed by Chuang and Swenson (1981) for a single inlet connecting to a bay and expanded by Wong and DiLorenzo (1988) to two connected bays and to multiple bays and inlets by Aretxabaleta et al. (2017). An analytical solution can be found for the entire system, with expressions for all the connections in the system. The model solves the along-channel depth-averaged momentum equation based on the balance between frictional effects and the elevation gradient between the offshore and bay areas and the continuity equation for the bay–channel system based on the changing volume of the bays as water flows through the inlets. The model also allows the estimation of the effect of the breach in Mantoloking during Hurricane Sandy. An analytical solution can be found by dividing the entire system into five sub-embayments (based on constrictions inside the bay system), resulting in a system of equations that includes 13 equations and unknowns (Appendix A).

Figure 2 Schematic diagram of the ocean–inlet–bay system: A_j is the surface area of the bays, η_j the sea level in the bays, η₀ the sea level in the ocean, and u_j is the velocity through channel j. The correspondence with the real bay system includes areas from the bays (Great Bay, A₁; Little Egg Harbor, A₂; Barnegat Bay, A₃; Toms River sub-embayment, A₄; north of Mantoloking, A₅), flow through inlets (Little Egg Inlet, u₁; Barnegat Inlet, u₄; the Point Pleasant Canal, u₈; and Mantoloking breach, u₇), and flow between bays (Tucker Island, u₂; Route 72, u₃; Bayville, u₅; Mantoloking, u₆). The location of the water level stations are indicated with dots, and the names and specifications are in Fig. 1 and Table 1.

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{\displaystyle \frac{\partial }{\partial t}}\left(\begin{array}{l}{u}_{\mathrm{1}}\\ u_{\mathrm{2}}\\ u_{\mathrm{3}}\\ u_{\mathrm{4}}\\ u_{\mathrm{5}}\\ u_{\mathrm{6}}\\ u_{\mathrm{7}}\\ u_{\mathrm{8}}\end{array}\right)& =g\left(\begin{array}{l}\mathrm{1}/L_{\mathrm{1}}\\ \mathrm{1}/L_{\mathrm{2}}\\ \mathrm{1}/L_{\mathrm{3}}\\ \mathrm{1}/L_{\mathrm{4}}\\ \mathrm{1}/L_{\mathrm{5}}\\ \mathrm{1}/L_{\mathrm{6}}\\ \mathrm{1}/L_{\mathrm{7}}\\ \mathrm{1}/L_{\mathrm{8}}\end{array}\right)\left(\begin{array}{l}{\mathit{\eta }}_{\mathrm{o}}{\mathit{\varphi }}_{\mathrm{LEI}}-{\mathit{\eta }}_{\mathrm{1}}\\ {\mathit{\eta }}_{\mathrm{1}}-{\mathit{\eta }}_{\mathrm{2}}\\ {\mathit{\eta }}_{\mathrm{2}}-{\mathit{\eta }}_{\mathrm{3}}\\ {\mathit{\eta }}_{\mathrm{o}}{\mathit{\varphi }}_{\mathrm{BI}}-{\mathit{\eta }}_{\mathrm{3}}\\ {\mathit{\eta }}_{\mathrm{3}}-{\mathit{\eta }}_{\mathrm{4}}\\ {\mathit{\eta }}_{\mathrm{4}}-{\mathit{\eta }}_{\mathrm{5}}\\ {\mathit{\eta }}_{\mathrm{o}}{\mathit{\varphi }}_{\mathrm{breach}}-{\mathit{\eta }}_{\mathrm{5}}\\ {\mathit{\eta }}_{\mathrm{o}}{\mathit{\varphi }}_{\mathrm{PPC}}-{\mathit{\eta }}_{\mathrm{5}}\end{array}\right)\\ & -\left(\begin{array}{l}{u}_{\mathrm{1}}{r}_{\mathrm{1}}/h_{\mathrm{1}}\\ u_{\mathrm{2}}{r}_{\mathrm{2}}/h_{\mathrm{2}}\\ u_{\mathrm{3}}{r}_{\mathrm{3}}/h_{\mathrm{3}}\\ u_{\mathrm{4}}{r}_{\mathrm{4}}/h_{\mathrm{4}}\\ u_{\mathrm{5}}{r}_{\mathrm{5}}/h_{\mathrm{5}}\\ u_{\mathrm{6}}{r}_{\mathrm{6}}/h_{\mathrm{6}}\\ u_{\mathrm{7}}{r}_{\mathrm{7}}/h_{\mathrm{7}}\\ u_{\mathrm{8}}{r}_{\mathrm{8}}/h_{\mathrm{8}}\end{array}\right)\end{array}\end{array}$$

$$\begin{array}{}\text{(2)}& \begin{array}{rl}& \left(\begin{array}{l}{A}_{\mathrm{1}}\\ A_{\mathrm{2}}\\ A_{\mathrm{3}}\\ A_{\mathrm{4}}\\ A_{\mathrm{5}}\end{array}\right){\displaystyle \frac{\partial }{\partial t}}\left(\begin{array}{l}{\mathit{\eta }}_{\mathrm{1}}\\ {\mathit{\eta }}_{\mathrm{2}}\\ {\mathit{\eta }}_{\mathrm{3}}\\ {\mathit{\eta }}_{\mathrm{4}}\\ {\mathit{\eta }}_{\mathrm{5}}\end{array}\right)=\\ & \left(\begin{array}{l}{h}_{\mathrm{1}}{W}_{\mathrm{1}}{u}_{\mathrm{1}}-h_{\mathrm{2}}{W}_{\mathrm{2}}{u}_{\mathrm{2}}\\ h_{\mathrm{2}}{W}_{\mathrm{2}}{u}_{\mathrm{2}}-h_{\mathrm{3}}{W}_{\mathrm{3}}{u}_{\mathrm{3}}\\ h_{\mathrm{3}}{W}_{\mathrm{3}}{u}_{\mathrm{3}}+h_{\mathrm{4}}{W}_{\mathrm{4}}{u}_{\mathrm{4}}-h_{\mathrm{5}}{W}_{\mathrm{5}}{u}_{\mathrm{5}}\\ h_{\mathrm{5}}{W}_{\mathrm{5}}{u}_{\mathrm{5}}-h_{\mathrm{6}}{W}_{\mathrm{6}}{u}_{\mathrm{6}}\\ h_{\mathrm{6}}{W}_{\mathrm{6}}{u}_{\mathrm{6}}+h_{\mathrm{7}}{W}_{\mathrm{7}}{u}_{\mathrm{7}}+h_{\mathrm{8}}{W}_{\mathrm{8}}{u}_{\mathrm{8}}\end{array}\right)\end{array}\end{array}$$

Here g is the gravitational acceleration, A_m is the surface area of sub-embayment m; η_m the sea level in the m sub-embayment; η_o the sea level in the ocean; h_n the water depth; and W_n the width, L_n the length, and r_n the linear drag coefficient of channel n. ϕ_LEI, ϕ_BI, ϕ_breach, and ϕ_PPC are the linear frequency-dependent relationships between the water levels at offshore proxy stations (Sandy Hook or Atlantic City) and the water level just offshore of Little Egg Inlet, Barnegat Inlet, the breach at Mantoloking caused by Hurricane Sandy, and the Point Pleasant Canal.

Assuming $\mathit{\eta }=\stackrel{\mathrm{̃}}{\mathit{\eta }}{e}^{i\mathit{\omega }t}$ and $u=\stackrel{\mathrm{̃}}{u}{e}^{i\mathit{\omega }t}$, where $\stackrel{\mathrm{̃}}{\mathit{\eta }}$ and $\stackrel{\mathrm{̃}}{u}$ represent the magnitude of the water level and velocity oscillations, respectively, to reduce the size of the equations, we can define $K_n=\frac{h_n{W}_{n}g}{L_n\left(\frac{r_n}{h_n}+i\mathit{\omega }\right)}$ for n=1, …, 8, as the relative contribution of each sub-embayment based on its geometric and frictional characteristics. Then, with the proper rearrangement, it yields the following equation.

$$\begin{array}{}\text{(3)}& {\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{3}}={\displaystyle \frac{{\displaystyle \begin{array}{c}{K}_{\mathrm{4}}\left({\mathit{\eta }}_{\mathrm{o}}\right){\mathit{\varphi }}_{\mathrm{BI}}+{\displaystyle \frac{\frac{K_{\mathrm{1}}{K}_{\mathrm{2}}{K}_{\mathrm{3}}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{o}}{\mathit{\varphi }}_{\mathrm{LEI}}}{i\mathit{\omega }{A}_{\mathrm{1}}+K_{\mathrm{1}}+K_{\mathrm{2}}}}{i\mathit{\omega }{A}_{\mathrm{2}}+K_{\mathrm{2}}+K_{\mathrm{3}}-\frac{K_{\mathrm{2}}{K}_{\mathrm{2}}}{i\mathit{\omega }{A}_{\mathrm{1}}+K_{\mathrm{1}}+K_{\mathrm{2}}}}}\\ +\frac{\frac{K_{\mathrm{5}}{K}_{\mathrm{6}}\left(K_{\mathrm{7}}{\mathit{\varphi }}_{\mathrm{breach}}+K_{\mathrm{8}}{\mathit{\varphi }}_{\mathrm{PPC}}\right){\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{o}}}{i\mathit{\omega }{A}_{\mathrm{5}}+K_{\mathrm{6}}+K_{\mathrm{7}}+K_{\mathrm{8}}}}{i\mathit{\omega }{A}_{\mathrm{4}}+K_{\mathrm{5}}+K_{\mathrm{6}}-\frac{K_{\mathrm{6}}{K}_{\mathrm{6}}}{i\mathit{\omega }{A}_{\mathrm{5}}+K_{\mathrm{6}}+K_{\mathrm{7}}+K_{\mathrm{8}}}}\end{array}}}{\begin{array}{c}i\mathit{\omega }{A}_{\mathrm{3}}+K_{\mathrm{3}}+K_{\mathrm{4}}+K_{\mathrm{5}}-\frac{K_{\mathrm{3}}{K}_{\mathrm{3}}}{i\mathit{\omega }{A}_{\mathrm{2}}+K_{\mathrm{2}}+K_{\mathrm{3}}-\frac{K_{\mathrm{2}}{K}_{\mathrm{2}}}{i\mathit{\omega }{A}_{\mathrm{1}}+K_{\mathrm{1}}+K_{\mathrm{2}}}}\\ -\frac{K_{\mathrm{5}}{K}_{\mathrm{5}}}{i\mathit{\omega }{A}_{\mathrm{4}}+K_{\mathrm{5}}+K_{\mathrm{6}}-\frac{K_{\mathrm{6}}{K}_{\mathrm{6}}}{i\mathit{\omega }{A}_{\mathrm{5}}+K_{\mathrm{6}}+K_{\mathrm{7}}+K_{\mathrm{8}}}}\end{array}}}\end{array}$$

The solution for the water level of the central sub-embayment can be used to recursively calculate the solutions for the rest of the sub-embayments as follows.

$$\begin{array}{}\text{(4)}& {\displaystyle }{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{2}}={\displaystyle \frac{K_{\mathrm{3}}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{3}}+\frac{K_{\mathrm{2}}{K}_{\mathrm{1}}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{o}}{\mathit{\varphi }}_{\mathrm{LEI}}}{i\mathit{\omega }{A}_{\mathrm{1}}+K_{\mathrm{1}}+K_{\mathrm{2}}}}{i\mathit{\omega }{A}_{\mathrm{2}}+K_{\mathrm{2}}+K_{\mathrm{3}}-\frac{K_{\mathrm{2}}{K}_{\mathrm{2}}}{i\mathit{\omega }{A}_{\mathrm{1}}+K_{\mathrm{1}}+K_{\mathrm{2}}}}}\text{(5)}& {\displaystyle }{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{4}}={\displaystyle \frac{K_{\mathrm{5}}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{3}}+\frac{K_{\mathrm{6}}\left(K_{\mathrm{7}}{\mathit{\varphi }}_{\mathrm{breach}}+K_{\mathrm{8}}{\mathit{\varphi }}_{\mathrm{PPC}}\right){\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{o}}}{i\mathit{\omega }{A}_{\mathrm{5}}+K_{\mathrm{6}}+K_{\mathrm{7}}+K_{\mathrm{8}}}}{i\mathit{\omega }{A}_{\mathrm{4}}+K_{\mathrm{5}}+K_{\mathrm{6}}-\frac{K_{\mathrm{6}}{K}_{\mathrm{6}}}{i\mathit{\omega }{A}_{\mathrm{5}}+K_{\mathrm{6}}+K_{\mathrm{7}}+K_{\mathrm{8}}}}}\text{(6)}& {\displaystyle }{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{1}}={\displaystyle \frac{K_{\mathrm{1}}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{o}}{\mathit{\varphi }}_{\mathrm{LEI}}+K_{\mathrm{2}}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{2}}}{i\mathit{\omega }{A}_{\mathrm{1}}+K_{\mathrm{1}}+K_{\mathrm{2}}}}\text{(7)}& {\displaystyle }{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{5}}={\displaystyle \frac{\left(K_{\mathrm{7}}{\mathit{\varphi }}_{\mathrm{breach}}+K_{\mathrm{8}}{\mathit{\varphi }}_{\mathrm{PPC}}\right){\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{o}}+K_{\mathrm{6}}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{4}}}{i\mathit{\omega }{A}_{\mathrm{5}}+K_{\mathrm{6}}+K_{\mathrm{7}}+K_{\mathrm{8}}}}\end{array}$$

The resulting expressions include all the sub-embayment and offshore exchanges under the same assumptions of the Chuang and Swenson (1981) model (e.g., no local influences, no overtopping).

4.2 Local wind impact on bay model

The contribution of local wind setup to the spatial distribution of water level inside the bay can be approximated following Wong and Moses-Hall (1998). The bay can be assumed to be a simple, long, and well-mixed embayment for which the cross-bay gradients and vertical stratification can be ignored. The linearized vertically integrated momentum and mass conservation equations are as follows:

$$\begin{array}{}\text{(8)}& {\displaystyle \frac{\partial U}{\partial t}}=-gh{\displaystyle \frac{\partial \mathit{\eta }}{\partial x}}+{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}\left({\mathit{\tau }}_{\mathrm{s}}-{\mathit{\tau }}_{\mathrm{b}}\right)=-gh{\displaystyle \frac{\partial \mathit{\eta }}{\partial x}}+{\mathit{\tau }}_w-r{\displaystyle \frac{U}{h}},\text{(9)}& {\displaystyle }\mathrm{and}{\displaystyle \frac{\partial U}{\partial x}}=-{\displaystyle \frac{\partial \mathit{\eta }}{\partial t}},\end{array}$$

where η is the water level in the bay, U is the depth-integrated along-bay velocity, h is the water depth, τ_s and τ_b are the surface and bottom dynamic stresses, respectively, and ρ₀ is the water density. ${\mathit{\tau }}_w={\mathit{\tau }}_{\mathrm{s}}/{\mathit{\rho }}_{\mathrm{0}}$ is the spatially invariant kinematic wind stress and r is linearized bottom friction.

Under the assumption of $\mathit{\eta }=\stackrel{\mathrm{̃}}{\mathit{\eta }}{e}^{i\mathit{\omega }t}$and $u=\stackrel{\mathrm{̃}}{u}{e}^{i\mathit{\omega }t}$, where ω is the cyclic frequency, the resulting equation is as follows:

$$\begin{array}{}\text{(10)}& {\displaystyle \frac{{\partial }^{\mathrm{2}}\stackrel{\mathrm{̃}}{\mathit{\eta }}}{\partial x^{\mathrm{2}}}}+\stackrel{\mathrm{̃}}{\mathit{\eta }}\left({\displaystyle \frac{{\mathit{\omega }}^{\mathrm{2}}}{gh}}-{\displaystyle \frac{i\mathit{\omega }r}{g{h}^{\mathrm{2}}}}\right)={\displaystyle \frac{{\partial }^{\mathrm{2}}\stackrel{\mathrm{̃}}{\mathit{\eta }}}{\partial x^{\mathrm{2}}}}+\stackrel{\mathrm{̃}}{\mathit{\eta }}{k}^{\mathrm{2}}=\mathrm{0},\end{array}$$

with boundary conditions $\mathit{\eta }\left(x=\mathrm{0},\mathit{\omega }\right)=\mathrm{0}$ assuming no offshore forcing at the entrance (this assumption will be revisited in the next section) and $\frac{\partial \stackrel{\mathrm{̃}}{\mathit{\eta }}\left(x=L,\mathit{\omega }\right)}{\partial x}=\frac{\stackrel{\mathrm{̃}}{{\mathit{\tau }}_w}\left(\mathit{\omega }\right)}{gh}$ assuming no flux at the head (x=L). $\stackrel{\mathrm{̃}}{{\mathit{\tau }}_w}\left(\mathit{\omega }\right)$ represents the magnitude of the kinematic wind stress that results in water level fluctuations at a specific frequency.

The solution is as follows:

$$\begin{array}{}\text{(11)}& {\displaystyle }\stackrel{\mathrm{̃}}{\mathit{\eta }}\left(\mathit{\omega }\right)={\displaystyle \frac{\stackrel{\mathrm{̃}}{{\mathit{\tau }}_w}\left(\mathit{\omega }\right)\sin \left(kx\right)}{ghk\cos \left(kL\right)}},\text{(12)}& {\displaystyle }\mathrm{where}k={\left({\displaystyle \frac{{\mathit{\omega }}^{\mathrm{2}}}{gh}}-{\displaystyle \frac{i\mathit{\omega }r}{g{h}^{\mathrm{2}}}}\right)}^{\mathrm{1}/\mathrm{2}}.\end{array}$$

The wavenumber k determines the spatial response of the transfer between the wind stress and the bay water level. The imaginary wavenumber part leads to exponential decay based on the frictional characteristics. The magnitude of water level is obtained from Eq. (11), while the ratio of real to imaginary parts provides information about the phase lag between wind stress and water level.

4.3 Combining local and remote effects

The local and remote effects can be combined in following the approach by Wong and Moses-Hall (1998). While bays can exhibit complex spatial responses to wind forcing especially in terms of currents (Csanady, 1973; Hunter and Hearn, 1987; Cioffi et al., 2005), the basic response can be summarized as the sum of local (wind) and remote (surge) forcings. The boundary condition for the local wind effect can be altered to account for the influence of the offshore water level, η_o. The resulting model is a modification of the wind effect model that considers the analytical offshore influence in Sect. 4.1. In a system with a single inlet, the solution can be simply stated, as in Wong and Moses-Hall (1998):

$$\begin{array}{}\text{(13)}& \stackrel{\mathrm{̃}}{\mathit{\eta }}\left(\mathit{\omega }\right)={\displaystyle \frac{\stackrel{\mathrm{̃}}{{\mathit{\tau }}_w}\left(\mathit{\omega }\right)\sin \left(kx\right)}{ghk\cos \left(kL\right)}}+{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{o}}\left(\mathit{\omega }\right){\displaystyle \frac{\cos \left(k\left(L-x\right)\right)}{\cos \left(kL\right)}}.\end{array}$$

In a system with multiple connections with the offshore, the solution can be more complex. One limitation of the approach is that it utilizes a linear friction approximation. To produce a better approximation that takes into account the complex frictional conditions of the bay (e.g., varying geometry, diverse bottom conditions, enhanced attenuation over submerged vegetation), we can take a numerical solution of the bay that resolves the tidal and subtidal water level conditions under realistic friction and adjust the spatial distribution of the transfer from offshore accordingly. As most of the water level variability in the bay is associated with the M₂ semidiurnal tidal constituent (Fig. 3) and the distribution of the tide has been properly validated in the numerical simulations (Defne and Ganju, 2015), we can take the spatial distribution of the M₂ tidal amplitude as a proxy for the internal frictional effects in the bay. Bottom friction caused by both wind-driven and tidal effects is considered in the numerical simulations. By adjusting the water level based on the numerical M₂ spatial distribution, we are approximating the complete frictional characteristics of the bay. The adjustment is applied to each of the sub-embayments using the following expression:

$$\begin{array}{}\text{(14)}& {\displaystyle \frac{{\widehat{\mathit{\eta }}}_j}{{\mathit{\eta }}_{\mathrm{o}}}}\left(x,\mathit{\omega }\right)=\mathrm{1}-{\displaystyle \frac{\left(\mathrm{1}-\frac{{\mathit{\eta }}_j}{{\mathit{\eta }}_{\mathrm{o}}}\left(\mathit{\omega }\right)\right)\left(\mathrm{1}-\frac{\mathit{\eta }\left(x,{\mathit{\omega }}_{M_{\mathrm{2}}}\right)}{\mathit{\eta }\left(\mathrm{offshore},{\mathit{\omega }}_{M_{\mathrm{2}}}\right)}\right)}{\left(\mathrm{1}-\frac{{\mathit{\eta }}_j}{{\mathit{\eta }}_{\mathrm{o}}}\left({\mathit{\omega }}_{M_{\mathrm{2}}}\right)\right)}},\end{array}$$

where η_j∕η_o(ω) is the transfer coefficient of the sub-embayment i (single value) at frequency ω, $\mathit{\eta }\left(x,{\mathit{\omega }}_{M_{\mathrm{2}}}\right)$ is the amplitude of the M₂ tidal fluctuations from the numerical model solution (spatially variable), and ${\widehat{\mathit{\eta }}}_j/{\mathit{\eta }}_{\mathrm{o}}\left(x,\mathit{\omega }\right)$ is the spatially variable adjusted transfer coefficient for sub-embayment j. The resulting adjusted transfer coefficients provide estimates of the spatial changes not only between adjacent sub-embayments but also inside each of the sub-embayments. The local wind effects on bay water level can be added to the impact from offshore fluctuations to obtain a combined local and remote water level response estimate. In cases with simultaneous presence of wind and offshore level fluctuations, the system can respond in a weakly nonlinear manner and departures from the presented basic addition of the process are expected.

Figure 3 Energy spectra at all stations, computed using a Hanning 29 d window with overlapping (50 %) data segments. O₁, K₁, N₂, M₂, and S₂ label the principal tidal frequencies and f the inertial frequency. The vertical shaded area indicates the frequencies corresponding to the storm band (2–5 d; cpd: cycles per day). See Table 1 for the key to station abbreviations and Fig. 1 for locations.

Table 2 Sum of energy (m²) in the different bands of the spectra computed for the period 2007–2018 (or longest available record) using a 29 d Hanning window with overlapping (50 %) data segments.

5 Results

5.1 Offshore transfer to bay

The maximum energy in water level spectra (Fig. 3, Table 2) was associated with the M₂ semidiurnal tidal constituent for the offshore proxies (SH and AC) and at the stations TUC and ETH in the southern part of the BBLEH area. For the locations in Barnegat Bay (WAR, SEH, MAN), maximum energy was in the low-frequency band. Large spectral energy also occurred in the other semidiurnal tidal frequencies (S₂ and N₂), the diurnal frequencies (O₁ and K₁), the storm band (periods between 2 and 5 d), and the low-frequency band (Table 2). The energy in the remaining bands exhibited average fluctuations less than 0.03 m in size offshore, while in the bay fluctuations were less than 0.01 m.

Figure 4 Observed transfer from Atlantic City to five bay stations: Mantoloking (MAN), Seaside Heights (SEH), Waretown (WAR), East Thorofare (ETH), and Tuckerton (TUC). Solid lines indicate transfers for the entire available record at each station. Dashed lines represent observed transfers for the period March–December 2012, for which numerical model solutions were available. The vertical shaded area indicates the frequencies corresponding to the storm band (2–5 d).

Figure 5 Comparison between observed (blue) and numerical model (black) transfers for the period when both are available (March–December 2012) at five bay stations. Uncertainty envelopes for the transfer coefficient (Bendat and Piersol, 1986) are provided for observed (light blue) and model (gray) estimates.

Transfer functions between Atlantic City (AC) and the five stations inside the bay (Fig. 4) for the longest available length of record showed a north to south gradient. The transfer of the offshore fluctuations was 50 %–80 % at periods between 2 and 5 d (storm band), except at Tuckerton (TUC; over 95 %). The transfers at diurnal periods were about 35 % for the three Barnegat Bay stations (WAR, SEH, MAN), about 45 % in Little Egg Harbor (ETH), and 80 % in Great Bay (TUC). For frequencies associated with the semidiurnal tides, the transfers were even more attenuated, with values of about 15 % (between 14 % and 16 %) inside Barnegat, between 30 % and 35 % at ETH, and between 60 % and 70 % at TUC. As the numerical model solution was only available for the period March–December 2012, the long-term (2007–2018) transfers were compared with shorter-term observations. The transfers were similar within the uncertainty envelopes for each station (not shown) for both datasets at most frequencies, except at Mantoloking (MAN), which showed enhanced transfers for periods between 1 and 5 d in the 2012 record, and at Seaside Heights (SEH), where transfers in the storm band were slightly attenuated during 2012. Transfer estimates using Sandy Hook (SH) as the offshore proxy instead of Atlantic City produced similar results (not shown). The transfer between stations AC and SH on the open coast (proxies for offshore fluctuations) has been shown to be close to one (Wong and Wilson, 1984; Aretxabaleta et al., 2014), confirming that the offshore forcing at all three inlets is about the same.

Transfers estimated from the numerical model solution (Fig. 5) showed similar magnitudes to the observed transfers within uncertainty envelopes provided by the Bendat and Piersol (1986) formulation at most frequencies. The observed and modeled transfer at diurnal and semidiurnal transfers were similar (within a few percentage points) at all stations, except that the model overestimated the semidiurnal transfer at TUC. Differences between modeled and observed estimates at MAN and SEH were only significant at frequencies that contained minimal energy. The model reproduced the enhanced transfer in the storm band at Mantoloking during 2012, suggesting a physical mechanism for the change that the model was able to capture but remains unexplained. The likely explanation is that the location of the Azores–Bermuda high-pressure system over the Atlantic in 2012 (Mattingly et al., 2015), associated with the negative phase of the North Atlantic Oscillation, resulted in average winds that lined up with the axis of the bay and caused enhanced wind setup in the northern part of the bay. The model overestimated the transfer at ETH in the storm band and underestimated the low-frequency transfer at Waretown. The likely cause for some of the discrepancies, especially at low frequencies, is the relatively short length of the available record.

The analytical model of offshore impact that considered five sub-embayments (Sect. 4.1) was fit to the observed transfers to obtain an estimate of linear friction. The fit considered the unevenly distributed energy spectra (Fig. 3) with adjusted weight to the semidiurnal and low-frequency components. The resulting friction was r=0.021 m s⁻¹. The associated frictional adjustment time, $t_{\mathrm{fr}}=h/r$, was about 1–5 min depending on the depth of the inlet. The analytical curves (Fig. 6) matched the observed transfer function shape at most frequencies. The analytical model with five sub-embayment domains captured the north–south spatial differences in transfer. The analytical model for the central Barnegat Bay sub-embayment (A₃) approximated the transfer estimates from observations at Waretown at most frequencies (less than 5 % in the storm band). The analytical model for Great Bay (A₁) adequately matched observed transfers at Tuckerton at diurnal and semidiurnal frequencies (less than 5 % difference) but underestimated the transfer in the storm band (model estimates about 90 %, while observations were above 95 %). Meanwhile, the analytical model for Little Egg Harbor (A₂) matched the observed transfers at ETH within the uncertainty envelope, except for a slight underprediction at diurnal frequencies (less than 5 %). The observed transfers at the northern stations (MAN and SEH) were reproduced by the analytical model (A₄, A₅, respectively) at diurnal and semidiurnal tides but were underpredicted for the higher storm band frequencies (5 %–10 % less transfer in frequencies close to 2 d periods) and overpredicted at low frequencies (about 10 % differences). The analytical model was used to explore the effect on transfer of the breach (U₇) at Mantoloking that opened during Hurricane Sandy. The transfers were so minimally affected that the curves are indistinguishable, with only a negligible enhancement (< 0.2 %) in transfer in the northernmost sub-embayment (A₅). The breach was too small and shallow for any significant volume transport to occur that would affect the large bay.

Figure 6 Observed transfer for the longest available record (solid lines) and best analytical model fit for each of the sub-embayments (dashed lines). The vertical shaded area indicates the frequencies corresponding to the storm band (2–5 d).

5.2 Local wind influence

The spectrum of the along-bay component (rotated 20^∘) of the wind (Fig. 7a) from the offshore buoy NDBC 44065 (10 years, 2008–2018) showed high energies in the storm band (2–5 d) and in low frequencies. The largest single peak of energy was associated with 24 h period oscillations, likely associated with sea breeze and matched energy values at 5 d frequencies. There was a small peak at inertial frequencies.

Figure 7 (a) Wind speed spectra for the along-bay wind component for NDBC 44065 buoy (2008–2018). (b) Kinematic wind stress contribution to local water level in the bay expressed as η∕τ_w (m⁻¹ s²), following the Wong and Moses-Hall (1998) formulation, as a function of distance from the southern edge of the bay ($x/L=\mathrm{1}$ is the northern edge of bay). (c) Kinematic wind stress contribution to water level as a function of frequency. The vertical shaded area indicates the frequencies corresponding to the storm band (2–5 d).

The local wind contribution to water level setup inside the bay was approximated using the Wong and Moses-Hall (1998) approach (Sect. 4.2). The resulting formulation showed the largest setup magnitudes near the head of the bay (e.g., northern part with wind blowing from the south) with a decay as distance from the head increased (Fig. 7b, c). The magnitude of the setup depended on the magnitude of the linear friction, with a smaller setup under stronger friction (Fig. 7b, c). The setup responded exponentially to fetch (distance), except over long durations and under low friction conditions, which were predominantly linear (Fig. 7b). The frictional control was less important at higher frequencies (Fig. 7c). As frequency increased there was less wind energy (Fig. 7a), so the frictional control is mostly important for low frequency and storm band wind fluctuations.

Figure 8 Local wind setup inside the bay based on the Wong and Moses-Hall (1998) formulation for a wind stress of 0.1 Pa during specific periods: (a) wind with a 1 d period (e.g., sea breeze), (b) 2 d wind, and (c) 5 d wind.

The resulting effect of the wind setup (or set-down) was small (less than 0.1 m with an along-bay wind stress of 0.1 Pa) for most of the domain (Fig. 8). The estimate assumed a linear friction of the same magnitude as in Sect. 5.1 (r=0.021 m s⁻¹). Under persistent wind stress of 0.1 Pa (about 8 m s⁻¹ wind speed) in the along-bay direction, the resulting setups varied depending on the frequency considered. Setup magnitudes over 0.2 m were estimated for the 5 d period wind (Fig. 8c), while under half of that magnitude was achieved for the 2 d persistent wind (Fig. 8b), and a much smaller water level setup (peak smaller than 0.1 m) was estimated for the sea breeze (Fig. 8a). During extreme events like Hurricane Sandy, under intense wind stress, two additional effects should be considered: the depth of the bay increases via the transfer of offshore surge, resulting in altered setup response (Sect. 4.2), and the frictional effect is enhanced (a larger linear friction would be needed) by the presence of wave-induced roughness.

6 Discussion

6.1 Spatially variable water level transfer

Following the approach described in Sect. 4.3, estimates of spatially variable water level impact from offshore can be calculated (Fig. 9). The M₂ tidal constituent transfer (Fig. 9a) showed a large north to south gradient, with values going from around 10 % in the north to over 80 % in the vicinity of Little Egg Inlet. The role of Barnegat Inlet in enhancing tidal transfer was greatly reduced, as most of the tide was attenuated in the inlet. The contribution of the Point Pleasant Canal was also small as expected from the tidal amplitudes (Chant, 2001; Defne and Ganju, 2015). The transfer in the storm band of 2 d fluctuations (Fig. 9b) also showed a strong north–south gradient, with values of about 50 % in Barnegat Bay, around 70 %–80 % in Little Egg Harbor, and larger values in Great Bay. The 5 d offshore fluctuations were transferred more efficiently into the bay (Fig. 9c) with values over 70 % in the entire bay, reaching 80 %–90 % in Little Egg Harbor, and over 90 % in Great Bay. Both storm band transfer estimates were controlled by the exchange through Little Egg Inlet, with very localized transfer enhancements in the vicinity of Barnegat Inlet and the Point Pleasant Canal. While the presented estimates used Atlantic City as the offshore proxy, similar results were obtained when Sandy Hook was used as the offshore reference (as expected from Aretxabaleta et al., 2014).

Figure 9 Spatially variable transfer function (percentage) of offshore fluctuations transferred into the bays using Atlantic City as an offshore proxy for three frequencies: (a) M₂ semidiurnal tide, (b) 2 d fluctuation in the storm band, and (c) 5 d fluctuation in the storm band. The filled circles represent the transfer estimate at each of the observed locations. Spatial pattern computed using the COAWST numerical solution.

When the magnitude of the fluctuations associated with a specific storm is available (η_o), then an estimate of the average water level in the bay during the storm can be obtained. For instance, for Hurricane Sandy the offshore surge associated with the storm was of the order of 2–3 m. Considering that the storm lasted for over a day, the water level transfer would have been above 50 % in Barnegat Bay and above 70 % in Little Egg Harbor. The resulting surge estimate in the bay was between 1 and 2 m just considering the exchange through the existing inlets. There was reported overtopping of the barrier island during the storm (McKenna et al., 2016) that might have further increased water level in the bay that the proposed method does not consider.

The wind setup effect inside the bay due to local wind can also be estimated for Hurricane Sandy using the approach in Sect. 4.2. Maximum wind stress during the storm was about 1 Pa. To obtain a maximum effect (worst-case scenario) it was assumed that the wind was persistently in the along-bay direction and that maximum stress was maintained for the duration of the storm. The maximum resulting water level considering the Wong and Moses-Hall method is linear with regard to wind stress magnitude (Fig. 7b) and would have been 10 times larger than the setup in Fig. 8b. The maximum wind setup would have been between 1 and 2 m, which was of the same order of magnitude as the surge produced from offshore sources. The cross-bay contribution to the wind setup during Hurricane Sandy was comparatively small, as wind direction was predominantly along-bay. Surge estimates from simple analytical formulations (State Committee for the Zuiderzee, 1926; Pugh, 1987) that do not consider storm duration produce similar magnitude results and are also dependent on the frictional response of the bay. Nonlinear interactions between local and remote effects may alter the total bay response, but these effects are likely second order.

Figure 10 Transfer estimate based on ADCIRC tidal database for three frequencies: (a) M₂ semidiurnal tide, (b) 2 d fluctuation in the storm band, and (c) 5 d fluctuation in the storm band.

6.2 Transfer based on tidal database

The approach thus far was based on the combination of observations, analytical models, and numerical models. In many systems, long-term observations that allow for the estimation of transfer coefficients might not be available. Also, numerical solutions for back-barrier bay systems tend to be computationally expensive and might not be available for the period of interest. We propose a relatively simple approach for some of these systems based on the availability of high-resolution tidal solutions for the system. The EC2015 ADCIRC tidal database (Szpilka et al., 2016) showed sufficient resolution (down to 13 m in some areas) in many bays along the east coast of the United States to resolve the tidal conditions with skill when compared to NOAA CO-OPS stations and historic International Hydrographic Organization (IHO) data. The EC2015 tidal database provides estimates for 37 tidal constituents. Based on those constituents and assuming that the totality of the offshore fluctuations at zero frequency reach the interior of the bay, an estimate can be provided for the storm band frequencies. A weighted least-squares interpolation in the frequency domain was performed based on the M₄, K₂, S₂, L₂, M₂, N₂, K₁, P₁, O₁, Q₁ tidal amplitudes ratios between each point of the ADCIRC domain inside the bay and a point in the offshore. Higher weight was given to zero frequency to nudge toward 100 % transfer at zero frequency. Estimates were calculated based on multiple locations inside the bay and averaged to achieve a more robust calculation and also obtain an approximation of the uncertainty associated with the estimate.

The resulting transfer estimates (Fig. 10) exhibited the same general spatial patterns shown in the previous estimates (Fig. 9) with slight differences. Some of the smaller features present in the COAWST numerical solution (Defne and Ganju, 2015) were not present in the ADCIRC EC2015 domain. The M₂ transfer estimate based on the tidal database (Fig. 10a) presented approximately the same magnitudes in most areas (average difference less than 3 %). The 5 d transfer (Fig. 10c) was also comparable to the solution described in Sect. 6.1 with values over 70 % in the entire domain and the southern areas exceeding 90 % transfer. The 2 d transfer from ADCIRC (Fig. 10b) was 5 %–10 % higher than the direct estimates (Fig. 9b). One of the benefits of the ADCIRC approach was the possibility of providing an approximation of the uncertainty (Fig. 11). The uncertainty estimate of the M₂ transfer was about 1 %–2 % (Fig. 11a) with higher values in the southern part of the domain. The 2 d transfer uncertainty (Fig. 11b) was above 4 % in Barnegat Bay in areas of larger discrepancy between the ADCIRC and complete approaches. The uncertainty estimates in the 5 d offshore water level transfer (Fig. 11c) in the northern part of the domain did not exceed 2.5 %.

The magnitude of the difference between the ADCIRC tidal database approach and the complete method highlighted in Sect. 4.3 was of the same order of magnitude or even smaller than the difference between observations and the analytical model (Fig. 6) or between observed transfers and numerical solutions (Fig. 5). This result emphasizes the validity of using the tidal database to calculate offshore transfer estimates, especially when water level observations inside the bay or numerical solutions are not available.

The effect of local wind setup will also need to be added to the ADCIRC-based estimate, especially during severe storms. The approach discussed in Sect. 5.2 or an even simpler surge calculation (e.g., from the steady state vertically averaged momentum equations, as in Pugh, 1987, from the traditional report of the State Committee for the Zuiderzee, 1926, or the updated frequency domain equivalent from Reef et al., 2018) could be used, and the resulting elevation could be added to the offshore transfer estimate obtain based on the ADCIRC tides. Thus, the production of bay water level predictions will require accurate wind forecast products and the quantification of the nonlinear interaction between local and remote effects.

Figure 11 Uncertainty in transfer estimate based on ADCIRC tidal database for three frequencies: (a) M₂ semidiurnal tide, (b) 2 d fluctuation in the storm band, and (c) 5 d fluctuation in the storm band.

6.3 Validity for flooding hazard assessments

The method presented offers a new methodology for coastal hazards assessment and risk analysis. While many methodologies are being used for open-coast regions (Thieler and Hammar-Klose, 1999; Kunreuther et al., 2000), vulnerability evaluation to coastal hazards in back-barrier bays remains underdeveloped. Evaluating bay hazards usually requires expensive computational simulations at appropriately high resolutions to characterize the spatial and temporal effects. The method presented here, using existing ADCIRC results, provides a less expensive approach that is able to properly estimate the spatial differences in vulnerability in response to flooding at different timescales (e.g., perigean spring tides, storms of different duration). It provides guidance for planning in response to “nuisance” flooding at a relatively low cost. It can be expanded to all back-barriers without the need to simulate each storm in each embayment while applying a consistent methodology.

Careful consideration needs to be given to the estimation of coastal hazards, especially for the forecast of intense storm effects. The inclusion of meticulously validated methodologies that consider both offshore influences (e.g., using the transfer estimated from ADCIRC tides) and local wind setup (e.g., Wong and Moses-Hall, 1998; Reef et al., 2018) is necessary. Skill assessment of storm hazard estimates using adequate observations is critical to avoid producing underpredictions or overpredictions of flooding and inundation.

As part of the general needs for hazard assessment (Ludwig et al., 2018), the important hazard characteristics that decision makers require include spatial extent, duration, and magnitude. The proposed methodology provides an approximation to both the extent of the area and the magnitude and also variations based on storm duration. Additionally, the fact that uncertainty estimates accompany the vulnerability provided by the present method enhances the potential value to decision makers. The extension to other bays in the United States will be included as part of the U.S. Geological Survey Coastal Change Hazards portal (USGS, 2018).

7 Summary

The results presented here demonstrate a strategy for estimating the impact of offshore sea level and local wind setup in back-barrier bay water levels. The transfer estimates from offshore to the bay water level used a combination of observations, analytical models based on appropriate simplifications of the bay system, and numerical simulations that provide the needed spatial distribution and more realistic frictional control.

The resulting maps of water level response to offshore forcing showed larger attenuation of the relatively higher-frequency fluctuations such as the semidiurnal tides. Smaller transfers were associated with shorter duration storms than longer duration storms, and transfer was most spatially uniform for storms of long duration. The description of the magnitude and spatial dependence of transfer on storm duration will assist planning for flooding in back-barrier bays.

In the specific case of the Barnegat Bay–Little Egg Harbor system, larger transfers were estimated for the southern embayments (Great Bay and Little Egg Harbor) when compared to Barnegat Bay. The reason for the difference was the dominant role of Little Egg Inlet (wider and deeper) in controlling the exchange between the offshore and bay systems. During relatively small storms, the contribution of local wind to the bay water level setup was smaller than the transfer from offshore fluctuations. During intense events, like hurricanes, local wind setup was of the same order of magnitude or even larger than offshore influences, depending on wind magnitude and especially on the relative angle of the wind with respect to bay orientation.

We introduced two approaches that depend on the availability of observations and numerical solutions. The approach requiring fewer data based on the ADCIRC tidal database provides spatial offshore transfer estimates and measures of uncertainty. In both cases, the inclusion of the local wind setup could be achieved based on simple surge analytical estimates. The approach that includes an analytical model allows for a simple tool to study the response of back-barrier bay systems to alternative conditions and forcing (e.g., geomorphic changes, changing duration of storms, sea level rise).

The proposed method represents an effective and inexpensive approach to flooding hazard evaluation in back-barrier bays and inland waters. The method provides detailed spatial estimates of vulnerabilities and uncertainties that could be an intuitive tool for coastal managers.

Appendix A

The offshore impact on the water level of the bay can be approximated with an analytical model that solves the linearized depth-averaged momentum equations. The system of equations for an idealized simplification of Barnegat Bay (Fig. 2) that includes five sub-embayments (based on constrictions inside the bay system) consists of 13 equations and unknowns.

$$\begin{array}{}\text{(A1)}& \begin{array}{rl}{\displaystyle \frac{\partial }{\partial t}}\left(\begin{array}{l}{u}_{\mathrm{1}}\\ u_{\mathrm{2}}\\ u_{\mathrm{3}}\\ u_{\mathrm{4}}\\ u_{\mathrm{5}}\\ u_{\mathrm{6}}\\ u_{\mathrm{7}}\\ u_{\mathrm{8}}\end{array}\right)& =g\left(\begin{array}{l}\mathrm{1}/L_{\mathrm{1}}\\ \mathrm{1}/L_{\mathrm{2}}\\ \mathrm{1}/L_{\mathrm{3}}\\ \mathrm{1}/L_{\mathrm{4}}\\ \mathrm{1}/L_{\mathrm{5}}\\ \mathrm{1}/L_{\mathrm{6}}\\ \mathrm{1}/L_{\mathrm{7}}\\ \mathrm{1}/L_{\mathrm{8}}\end{array}\right)\left(\begin{array}{l}{\mathit{\eta }}_{\mathrm{o}}{\mathit{\varphi }}_{\mathrm{LEI}}-{\mathit{\eta }}_{\mathrm{1}}\\ {\mathit{\eta }}_{\mathrm{1}}-{\mathit{\eta }}_{\mathrm{2}}\\ {\mathit{\eta }}_{\mathrm{2}}-{\mathit{\eta }}_{\mathrm{3}}\\ {\mathit{\eta }}_{\mathrm{o}}{\mathit{\varphi }}_{\mathrm{BI}}-{\mathit{\eta }}_{\mathrm{3}}\\ {\mathit{\eta }}_{\mathrm{3}}-{\mathit{\eta }}_{\mathrm{4}}\\ {\mathit{\eta }}_{\mathrm{4}}-{\mathit{\eta }}_{\mathrm{5}}\\ {\mathit{\eta }}_{\mathrm{o}}{\mathit{\varphi }}_{\mathrm{breach}}-{\mathit{\eta }}_{\mathrm{5}}\\ {\mathit{\eta }}_{\mathrm{o}}{\mathit{\varphi }}_{\mathrm{PPC}}-{\mathit{\eta }}_{\mathrm{5}}\end{array}\right)\\ & -\left(\begin{array}{l}{u}_{\mathrm{1}}{r}_{\mathrm{1}}/h_{\mathrm{1}}\\ u_{\mathrm{2}}{r}_{\mathrm{2}}/h_{\mathrm{2}}\\ u_{\mathrm{3}}{r}_{\mathrm{3}}/h_{\mathrm{3}}\\ u_{\mathrm{4}}{r}_{\mathrm{4}}/h_{\mathrm{4}}\\ u_{\mathrm{5}}{r}_{\mathrm{5}}/h_{\mathrm{5}}\\ u_{\mathrm{6}}{r}_{\mathrm{6}}/h_{\mathrm{6}}\\ u_{\mathrm{7}}{r}_{\mathrm{7}}/h_{\mathrm{7}}\\ u_{\mathrm{8}}{r}_{\mathrm{8}}/h_{\mathrm{8}}\end{array}\right)\end{array}\end{array}$$

$$\begin{array}{}\text{(A2)}& \begin{array}{rl}& \left(\begin{array}{l}{A}_{\mathrm{1}}\\ A_{\mathrm{2}}\\ A_{\mathrm{3}}\\ A_{\mathrm{4}}\\ A_{\mathrm{5}}\end{array}\right){\displaystyle \frac{\partial }{\partial t}}\left(\begin{array}{l}{\mathit{\eta }}_{\mathrm{1}}\\ {\mathit{\eta }}_{\mathrm{2}}\\ {\mathit{\eta }}_{\mathrm{3}}\\ {\mathit{\eta }}_{\mathrm{4}}\\ {\mathit{\eta }}_{\mathrm{5}}\end{array}\right)=\\ & \left(\begin{array}{l}{h}_{\mathrm{1}}{W}_{\mathrm{1}}{u}_{\mathrm{1}}-h_{\mathrm{2}}{W}_{\mathrm{2}}{u}_{\mathrm{2}}\\ h_{\mathrm{2}}{W}_{\mathrm{2}}{u}_{\mathrm{2}}-h_{\mathrm{3}}{W}_{\mathrm{3}}{u}_{\mathrm{3}}\\ h_{\mathrm{3}}{W}_{\mathrm{3}}{u}_{\mathrm{3}}+h_{\mathrm{4}}{W}_{\mathrm{4}}{u}_{\mathrm{4}}-h_{\mathrm{5}}{W}_{\mathrm{5}}{u}_{\mathrm{5}}\\ h_{\mathrm{5}}{W}_{\mathrm{5}}{u}_{\mathrm{5}}-h_{\mathrm{6}}{W}_{\mathrm{6}}{u}_{\mathrm{6}}\\ h_{\mathrm{6}}{W}_{\mathrm{6}}{u}_{\mathrm{6}}+h_{\mathrm{7}}{W}_{\mathrm{7}}{u}_{\mathrm{7}}+h_{\mathrm{8}}{W}_{\mathrm{8}}{u}_{\mathrm{8}}\end{array}\right)\end{array}\end{array}$$

${\mathit{\varphi }}_{\mathrm{LEI}},{\mathit{\varphi }}_{\mathrm{BI}},{\mathit{\varphi }}_{\mathrm{breach}}$, and ϕ_PPC are linear frequency-dependent relationships between the water levels at offshore proxy stations (Sandy Hook or Atlantic City) and the water level just offshore of Little Egg Inlet, Barnegat Inlet, the breach at Mantoloking caused by Hurricane Sandy, and the Point Pleasant Canal.

By performing Fourier transforms on the momentum equations ($\mathit{\eta }=\stackrel{\mathrm{̃}}{\mathit{\eta }}{e}^{i\mathit{\omega }t}$ and $u=\stackrel{\mathrm{̃}}{u}{e}^{i\mathit{\omega }t}$, where $\stackrel{\mathrm{̃}}{\mathit{\eta }}$ and $\stackrel{\mathrm{̃}}{u}$ represent the magnitude of the water level and velocity oscillations, respectively), we obtain the following equations.

$$\begin{array}{}\text{(A3)}& \begin{array}{rl}i\mathit{\omega }\left(\begin{array}{l}\stackrel{\mathrm{̃}}{u_{\mathrm{1}}}\\ \stackrel{\mathrm{̃}}{u_{\mathrm{2}}}\\ \stackrel{\mathrm{̃}}{u_{\mathrm{3}}}\\ \stackrel{\mathrm{̃}}{u_{\mathrm{4}}}\\ \stackrel{\mathrm{̃}}{u_{\mathrm{5}}}\\ \stackrel{\mathrm{̃}}{u_{\mathrm{6}}}\\ \stackrel{\mathrm{̃}}{u_{\mathrm{7}}}\\ \stackrel{\mathrm{̃}}{u_{\mathrm{8}}}\end{array}\right)& =g\left(\begin{array}{l}\mathrm{1}/L_{\mathrm{1}}\\ \mathrm{1}/L_{\mathrm{2}}\\ \mathrm{1}/L_{\mathrm{3}}\\ \mathrm{1}/L_{\mathrm{4}}\\ \mathrm{1}/L_{\mathrm{5}}\\ \mathrm{1}/L_{\mathrm{6}}\\ \mathrm{1}/L_{\mathrm{7}}\\ \mathrm{1}/L_{\mathrm{8}}\end{array}\right)\left(\begin{array}{l}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{o}}{\mathit{\varphi }}_{\mathrm{LEI}}-{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{1}}\\ {\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{1}}-{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{2}}\\ {\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{2}}-{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{3}}\\ {\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{o}}{\mathit{\varphi }}_{\mathrm{BI}}-{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{3}}\\ {\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{3}}-{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{4}}\\ {\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{4}}-{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{5}}\\ {\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{o}}{\mathit{\varphi }}_{\mathrm{breach}}-{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{5}}\\ {\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{o}}{\mathit{\varphi }}_{\mathrm{PPC}}-{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{5}}\end{array}\right)\\ & -\left(\begin{array}{l}\stackrel{\mathrm{̃}}{u_{\mathrm{1}}}{r}_{\mathrm{1}}/h_{\mathrm{1}}\\ \stackrel{\mathrm{̃}}{u_{\mathrm{2}}}{r}_{\mathrm{2}}/h_{\mathrm{2}}\\ \stackrel{\mathrm{̃}}{u_{\mathrm{3}}}{r}_{\mathrm{3}}/h_{\mathrm{3}}\\ \stackrel{\mathrm{̃}}{u_{\mathrm{4}}}{r}_{\mathrm{4}}/h_{\mathrm{4}}\\ \stackrel{\mathrm{̃}}{u_{\mathrm{5}}}{r}_{\mathrm{5}}/h_{\mathrm{5}}\\ \stackrel{\mathrm{̃}}{u_{\mathrm{6}}}{r}_{\mathrm{6}}/h_{\mathrm{6}}\\ \stackrel{\mathrm{̃}}{u_{\mathrm{7}}}{r}_{\mathrm{7}}/h_{\mathrm{7}}\\ \stackrel{\mathrm{̃}}{u_{\mathrm{8}}}{r}_{\mathrm{8}}/h_{\mathrm{8}}\end{array}\right),\end{array}\end{array}$$

and then

$$\begin{array}{}\text{(A4)}& \begin{array}{rl}\left(\begin{array}{l}\stackrel{\mathrm{̃}}{u_{\mathrm{1}}}\\ \stackrel{\mathrm{̃}}{u_{\mathrm{2}}}\\ \stackrel{\mathrm{̃}}{u_{\mathrm{3}}}\\ \stackrel{\mathrm{̃}}{u_{\mathrm{4}}}\\ \stackrel{\mathrm{̃}}{u_{\mathrm{5}}}\\ \stackrel{\mathrm{̃}}{u_{\mathrm{6}}}\\ \stackrel{\mathrm{̃}}{u_{\mathrm{7}}}\\ \stackrel{\mathrm{̃}}{u_{\mathrm{8}}}\end{array}\right)& =g\left(\begin{array}{l}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{o}}{\mathit{\varphi }}_{\mathrm{LEI}}-{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{1}}\\ {\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{1}}-{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{2}}\\ {\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{2}}-{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{3}}\\ {\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{o}}{\mathit{\varphi }}_{\mathrm{BI}}-{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{3}}\\ {\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{3}}-{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{4}}\\ {\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{4}}-{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{5}}\\ {\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{o}}{\mathit{\varphi }}_{\mathrm{breach}}-{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{5}}\\ {\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{o}}{\mathit{\varphi }}_{\mathrm{PPC}}-{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{5}}\end{array}\right)/\\ & \left[\left(\begin{array}{l}{L}_{\mathrm{1}}\\ L_{\mathrm{2}}\\ L_{\mathrm{3}}\\ L_{\mathrm{4}}\\ L_{\mathrm{5}}\\ L_{\mathrm{6}}\\ L_{\mathrm{7}}\\ L_{\mathrm{8}}\end{array}\right)\left(\begin{array}{l}i\mathit{\omega }+r_{\mathrm{1}}/h_{\mathrm{1}}\\ i\mathit{\omega }+r_{\mathrm{2}}/h_{\mathrm{2}}\\ i\mathit{\omega }+r_{\mathrm{3}}/h_{\mathrm{3}}\\ i\mathit{\omega }+r_{\mathrm{4}}/h_{\mathrm{4}}\\ i\mathit{\omega }+r_{\mathrm{5}}/h_{\mathrm{5}}\\ i\mathit{\omega }+r_{\mathrm{6}}/h_{\mathrm{6}}\\ i\mathit{\omega }+r_{\mathrm{7}}/h_{\mathrm{7}}\\ i\mathit{\omega }+r_{\mathrm{8}}{h}_{\mathrm{8}}\end{array}\right)\right].\end{array}\end{array}$$

Performing the Fourier transform on the continuity equations (Eq. A2) and substituting the velocity values from Eq. (A4) we obtain the following equation.

$$\begin{array}{}\text{(A5)}& \begin{array}{rl}& i\mathit{\omega }\left(\begin{array}{l}{A}_{\mathrm{1}}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{1}}\\ A_{\mathrm{2}}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{2}}\\ A_{\mathrm{3}}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{3}}\\ A_{\mathrm{4}}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{4}}\\ A_{\mathrm{5}}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{5}}\end{array}\right)=g\\ & \left(\begin{array}{l}{h}_{\mathrm{1}}{W}_{\mathrm{1}}\left({\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{o}}{\mathit{\varphi }}_{\mathrm{LEI}}-{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{1}}\right)/L_{\mathrm{1}}\left(i\mathit{\omega }+r_{\mathrm{1}}/h_{\mathrm{1}}\right)-h_{\mathrm{2}}{W}_{\mathrm{2}}\left({\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{1}}-{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{2}}\right)/\\ L_{\mathrm{2}}\left(i\mathit{\omega }+r_{\mathrm{2}}/h_{\mathrm{2}}\right)\\ h_{\mathrm{2}}{W}_{\mathrm{2}}\left({\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{1}}-{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{2}}\right)/L_{\mathrm{2}}\left(i\mathit{\omega }+r_{\mathrm{2}}/h_{\mathrm{2}}\right)-h_{\mathrm{3}}{W}_{\mathrm{3}}\left({\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{2}}-{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{3}}\right)/\\ L_{\mathrm{3}}\left(i\mathit{\omega }+r_{\mathrm{3}}/h_{\mathrm{3}}\right)\\ h_{\mathrm{3}}{W}_{\mathrm{3}}\left({\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{2}}-{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{3}}\right)/L_{\mathrm{3}}\left(i\mathit{\omega }+r_{\mathrm{3}}/h_{\mathrm{3}}\right)+h_{\mathrm{4}}{W}_{\mathrm{4}}\left({\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{o}}{\mathit{\varphi }}_{\mathrm{BI}}-{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{3}}\right)/\\ L_{\mathrm{4}}\left(i\mathit{\omega }+r_{\mathrm{4}}/h_{\mathrm{4}}\right)-h_{\mathrm{5}}{W}_{\mathrm{5}}\left({\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{3}}-{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{4}}\right)/L_{\mathrm{5}}\left(i\mathit{\omega }+r_{\mathrm{5}}/h_{\mathrm{5}}\right)\\ h_{\mathrm{5}}{W}_{\mathrm{5}}\left({\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{3}}-{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{4}}\right)/L_{\mathrm{5}}\left(i\mathit{\omega }+r_{\mathrm{5}}/h_{\mathrm{5}}\right)-h_{\mathrm{6}}{W}_{\mathrm{6}}\left({\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{4}}-{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{5}}\right)/\\ L_{\mathrm{6}}\left(i\mathit{\omega }+r_{\mathrm{6}}/h_{\mathrm{6}}\right)\\ h_{\mathrm{6}}{W}_{\mathrm{6}}\left({\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{4}}-{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{5}}\right)/L_{\mathrm{6}}\left(i\mathit{\omega }+r_{\mathrm{6}}/h_{\mathrm{6}}\right)+h_{\mathrm{7}}{W}_{\mathrm{7}}\left({\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{o}}{\mathit{\varphi }}_{\mathrm{breach}}-{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{5}}\right)/\\ L_{\mathrm{7}}\left(i\mathit{\omega }+r_{\mathrm{7}}/h_{\mathrm{7}}\right)+h_{\mathrm{8}}{W}_{\mathrm{8}}\left({\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{0}}{\mathit{\varphi }}_{\mathrm{PPC}}-{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{5}}\right)/L_{\mathrm{8}}\left(i\mathit{\omega }+r_{\mathrm{8}}/h_{\mathrm{8}}\right)\end{array}\right).\end{array}\end{array}$$

Then, we can define $K_n=\frac{h_n{W}_{n}g}{L_n\left(\frac{r_n}{h_n}+i\mathit{\omega }\right)}$ for each channel n=1, …, 8 and with the proper rearrangement we get the following equation:

$$\begin{array}{}\text{(A6)}& \left(\begin{array}{l}i\mathit{\omega }{A}_{\mathrm{1}}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{1}}=K_{\mathrm{1}}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{o}}{\mathit{\varphi }}_{\mathrm{LEI}}-K_{\mathrm{1}}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{1}}-K_{\mathrm{2}}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{1}}+K_{\mathrm{2}}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{2}},\\ i\mathit{\omega }{A}_{\mathrm{2}}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{2}}=K_{\mathrm{2}}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{1}}-K_{\mathrm{2}}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{2}}-K_{\mathrm{3}}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{2}}+K_{\mathrm{3}}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{3}},\\ i\mathit{\omega }{A}_{\mathrm{3}}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{3}}=K_{\mathrm{3}}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{2}}-K_{\mathrm{3}}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{3}}+K_{\mathrm{4}}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{o}}{\mathit{\varphi }}_{\mathrm{BI}}-K_{\mathrm{4}}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{3}},\\ -K_{\mathrm{5}}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{3}}+K_{\mathrm{5}}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{4}},\\ i\mathit{\omega }{A}_{\mathrm{4}}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{4}}=K_{\mathrm{5}}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{3}}-K_{\mathrm{5}}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{4}}-K_{\mathrm{6}}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{4}}+K_{\mathrm{6}}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{5}},\\ i\mathit{\omega }{A}_{\mathrm{5}}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{5}}=K_{\mathrm{6}}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{4}}-K_{\mathrm{6}}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{5}}+K_{\mathrm{7}}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{o}}{\mathit{\varphi }}_{\mathrm{breach}},\\ -K_{\mathrm{7}}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{5}}+K_{\mathrm{8}}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{o}}{\mathit{\varphi }}_{\mathrm{PPC}}-K_{\mathrm{8}}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{5}}\end{array}\right).\end{array}$$

The system of equations can be solved by substitution.

$$\begin{array}{}\text{(A7)}& {\displaystyle }{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{5}}={\displaystyle \frac{\left(K_{\mathrm{7}}{\mathit{\varphi }}_{\mathrm{breach}}+K_{\mathrm{8}}{\mathit{\varphi }}_{\mathrm{PPC}}\right){\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{o}}+K_{\mathrm{6}}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{4}}}{i\mathit{\omega }{A}_{\mathrm{5}}+K_{\mathrm{6}}+K_{\mathrm{7}}+K_{\mathrm{8}}}},\text{(A8)}& {\displaystyle }{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{4}}={\displaystyle \frac{K_{\mathrm{5}}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{3}}+\frac{K_{\mathrm{6}}\left(K_{\mathrm{7}}{\mathit{\varphi }}_{\mathrm{breach}}+K_{\mathrm{8}}{\mathit{\varphi }}_{\mathrm{PPC}}\right){\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{o}}}{i\mathit{\omega }{A}_{\mathrm{5}}+K_{\mathrm{6}}+K_{\mathrm{7}}+K_{\mathrm{8}}}}{i\mathit{\omega }{A}_{\mathrm{4}}+K_{\mathrm{5}}+K_{\mathrm{6}}-\frac{K_{\mathrm{6}}{K}_{\mathrm{6}}}{i\mathit{\omega }{A}_{\mathrm{5}}+K_{\mathrm{6}}+K_{\mathrm{7}}+K_{\mathrm{8}}}}},\text{(A9)}& {\displaystyle }{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{1}}={\displaystyle \frac{K_{\mathrm{1}}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{o}}{\mathit{\varphi }}_{\mathrm{LEI}}+K_{\mathrm{2}}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{2}}}{i\mathit{\omega }{A}_{\mathrm{1}}+K_{\mathrm{1}}+K_{\mathrm{2}}}},\text{(A10)}& {\displaystyle }{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{2}}={\displaystyle \frac{K_{\mathrm{3}}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{3}}+\frac{K_{\mathrm{2}}{K}_{\mathrm{1}}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{o}}{\mathit{\varphi }}_{\mathrm{LEI}}}{i\mathit{\omega }{A}_{\mathrm{1}}+K_{\mathrm{1}}+K_{\mathrm{2}}}}{i\mathit{\omega }{A}_{\mathrm{2}}+K_{\mathrm{2}}+K_{\mathrm{3}}-\frac{K_{\mathrm{2}}{K}_{\mathrm{2}}}{i\mathit{\omega }{A}_{\mathrm{1}}+K_{\mathrm{1}}+K_{\mathrm{2}}}}},\end{array}$$

and finally,

$$\begin{array}{}\text{(A11)}& {\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{3}}={\displaystyle \frac{{\displaystyle \begin{array}{c}{K}_{\mathrm{4}}{\mathit{\eta }}_{\mathrm{o}}{\mathit{\varphi }}_{\mathrm{BI}}+{\displaystyle \frac{\frac{K_{\mathrm{1}}{K}_{\mathrm{2}}{K}_{\mathrm{3}}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{o}}{\mathit{\varphi }}_{\mathrm{LEI}}}{i\mathit{\omega }{A}_{\mathrm{1}}+K_{\mathrm{1}}+K_{\mathrm{2}}}}{i\mathit{\omega }{A}_{\mathrm{2}}+K_{\mathrm{2}}+K_{\mathrm{3}}-\frac{K_{\mathrm{2}}{K}_{\mathrm{2}}}{i\mathit{\omega }{A}_{\mathrm{1}}+K_{\mathrm{1}}+K_{\mathrm{2}}}}}\\ +\frac{\frac{K_{\mathrm{5}}{K}_{\mathrm{6}}\left(K_{\mathrm{7}}{\mathit{\varphi }}_{\mathrm{breach}}+K_{\mathrm{8}}{\mathit{\varphi }}_{\mathrm{PPC}}\right){\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{\mathrm{o}}}{i\mathit{\omega }{A}_{\mathrm{5}}+K_{\mathrm{6}}+K_{\mathrm{7}}+K_{\mathrm{8}}}}{i\mathit{\omega }{A}_{\mathrm{4}}+K_{\mathrm{5}}+K_{\mathrm{6}}-\frac{K_{\mathrm{6}}{K}_{\mathrm{6}}}{i\mathit{\omega }{A}_{\mathrm{5}}+K_{\mathrm{6}}+K_{\mathrm{7}}+K_{\mathrm{8}}}}\end{array}}}{\begin{array}{c}i\mathit{\omega }{A}_{\mathrm{3}}+K_{\mathrm{3}}+K_{\mathrm{4}}+K_{\mathrm{5}}-\frac{K_{\mathrm{3}}{K}_{\mathrm{3}}}{i\mathit{\omega }{A}_{\mathrm{2}}+K_{\mathrm{2}}+K_{\mathrm{3}}-\frac{K_{\mathrm{2}}{K}_{\mathrm{2}}}{i\mathit{\omega }{A}_{\mathrm{1}}+K_{\mathrm{1}}+K_{\mathrm{2}}}}\\ -\frac{K_{\mathrm{5}}{K}_{\mathrm{5}}}{i\mathit{\omega }{A}_{\mathrm{4}}+K_{\mathrm{5}}+K_{\mathrm{6}}-\frac{K_{\mathrm{6}}{K}_{\mathrm{6}}}{i\mathit{\omega }{A}_{\mathrm{5}}+K_{\mathrm{6}}+K_{\mathrm{7}}+K_{\mathrm{8}}}}\end{array}}}.\end{array}$$

As it only depends on the offshore water level, the solution for the water level of the central sub-embayment (Eq. A11) can be used to recursively calculate the solutions for the rest of the sub-embayments by substituting the water level estimates into Eqs. (A7), (A8), (A9), and (A10).