1 Introduction

Along the U.S. West Coast, high ocean water levels and coastal flooding can occur when energetic winter swell events coincide with high tides. Serafin et al. (2017) estimate that wave runup contributes about half of the maximum total water levels (TWLs) for a study site in California, with the remainder from tides and mean sea-level anomalies. Improving wave runup predictions for the broad range of shorelines that typify the west coast remains a primary concern for coastal resiliency. Moreover, over the past century sea level has risen ~ 0.2 m at the San Diego and San Francisco tide gauges along the Pacific coast, with 0.6 m or more of additional rise probable by 2100 (Griggs et al. 2017; Table 5a, RCP 8.5 scenario). As sea levels continue to rise, more combinations of tides and wave events will exceed critical levels, significantly increasing the duration, intensity, and frequency of flooding (e.g., Vitousek et al. 2017, Sweet et al. 2017, Takerhani et al. 2020). Southern California is particularly vulnerable to sea-level rise and the 100-year event is expected to become annual by 2050 (Tebaldi et al. 2012).

Estimates of vulnerability to storm impacts are needed for nowcasts and forecasts for a wide range of shoreline types. Methods have been developed to identify and rank coastal vulnerability over regional scales and with various climate change scenarios. The RISC-KIT project provides planning tools to coastal managers, including processes to identify and rank erosion and flooding hotspots both now and with projected climate changes (van Dongeren et al. 2017; Ferreira et al. 2018; Armaroli and Duo 2018, and others). Real (and short lead) time TWL forecasts are becoming operational, such as the USGS Water Level Viewer, which spans the US East and Gulf Coasts and features 6-day forecasts of TWL and interaction with dunes (Aretxabaleta et al. 2019; Doran et al. 2019). Stokes et al (2021) provide overtopping warnings for 1000 km of southwest England coast, including sections defended with seawalls. In addition to regional models that span hundreds of km of shoreline and include a range of beach types and shoreline infrastructure, site specific forecast systems have been developed for wave runup (e.g., Atkinson et al. 2017, Luccio et al. 2018) and TWL (Coastal Data Information Program, https://cdip.ucsd.edu/themes/cdip?pb=1&d2=p112; PacIOOS, https://www.pacioos.hawaii.edu/shoreline-category/highsea/, Guiles et al. 2019). Some models neglect incident wave contribution to TWL (Vousdoukas et al. 2016), while others use empirical runup formulas. Among dozens of formulations, the most widely used is Stockdon et al. (2006) (hereafter S06), where the effect of profile shape on runup is included with a foreshore slope. S06 (or similar) model errors can be unacceptably large when applied to beaches on which they have not been calibrated (Atkinson et al. 2017; Fiedler et al. 2020). Here we describe the performance of a flood warning system where generic empirical estimates (S06) have been replaced with numerical modeling of runup with site-specific historical storm waves and bathymetry (Fiedler et al. 2020).

Low-lying zones around estuaries, lagoons, and wetlands are particularly susceptible to wave and tidal flooding. Such is the case in San Diego County where frequent winter flooding occurs at Imperial Beach (IB) adjacent to the floodplain of the Tijuana River Estuary (Fig. 1). Wave-driven overtopping events at IB currently are short-lived around peak high tides (Gallien 2016); however, they damage infrastructure by salt-water exposure and sand deposition, and disrupt traffic (Fig. 2). A recent sea-level rise assessment for the city identified long-term vulnerabilities impacting land use, roads, public transportation, wastewater systems, stormwater runoff, schools, and hazardous materials (Revell Coastal 2016). Estimates of the impact to Imperial Beach’s Gross Domestic Product (GDP) range from $9.4 M for a 2-m rise in sea level accompanied by a single 1-yr return period storm event, to $20.9 M for a 2-m rise in sea level combined with a 100-yr event (Colgan and DePaolis 2018).

Fig. 1
figure 1

Imperial Beach map showing offshore isobaths relative to MSL, names of often flooded streets, and sensor locations during the Resilient Futures observations. The bulge in bathymetry (e.g., 12-m contour north of present Tijuana River mouth) can focus wave energy and drive sediment transport

Fig. 2
figure 2

(Upper) South Imperial Beach during a moderate flood on 18 Jan 2020. LiDAR scans of the swash zone (discussed in F20) were collected from a condominium balcony at the seaward end of Cortez Ave (Arrows). (Middle) Entry points for street flooding include: directly through public access through the berm at Seacoast Dr, over riprap rock at Cortez Ave, and through an access point in the seawall at Palm Plaza. (Lower) Flooding occurs landward of these entry points

Wave-driven, high-tide flooding is common at IB in winter, but accurate flood forecasts, and their uncertainties, are lacking. To fill this knowledge gap, the Resilient Futures project was initiated to develop an early warning flood forecast capability for IB. The three major contributors to TWL—tidal elevation, wave-driven runup, and non-tidal residual water level—are predicted up to 6 days in advance. As part of the Resilient Futures project, additional wave, water level, and beach elevation observations have been collected to test runup models, quantify forecast uncertainties, and examine TWL variation along the IB shoreline. An online forecast tool, described here, is available to help prepare for flooding events.

Forecasting the wave runup component of TWL, the main focus of this study, requires reliable incident-wave forecasts, a model to transform incident waves to swash-zone runup, and specification of the nearshore bathymetry, which can change significantly over a winter, and even during a single storm. Here we employ a surrogate of a dynamical numerical model, optimized for conditions at IB (Fiedler et al. 2020). The predictions of the runup model and the resulting TWL estimate are compared with recent flood observations. IB is well suited to evaluate a TWL system because dynamical runup models have been shown skillful (Gallien 2016; Fiedler et al. 2020), beach elevation measurements are available through ongoing projects (Ludka et al. 2019), and the incident wave field is well characterized by a regional wave modeling system based on offshore buoy observations (O’Reilly et al. 2016). Application and calibration of the warning system approach to other locations requires incident wave forecasts, characterization of the local wave climate and beach morphology, and observations of storm wave runup.

2 Study area and data

Imperial Beach, the southernmost city on the California coast, borders water on three sides with San Diego Bay to the north, the Tijuana River Estuary to the south, and the Pacific Ocean to the west (Fig. 2 upper). While flooding may occur along all three boundaries, the focus here is on the winter swell and extreme TWLs that impact the oceanfront. The coastal strip features a relatively steep sand beach (0.07–0.14 slope in winter) intersected by the IB pier and a groin a few 100-m north of the pier (Fig. 1). The back beach rises to ~ 3 m above mean sea level (MSL), with the elevation dropping further inland, with especially low points at the northern (near Palm Ave) and southern ends of IB. The southern portion of the city consists of a narrow strip of land along Seacoast Dr that lies between the ocean and estuary. Homes and roadways along the peninsula experience overwash and flooding, notably near Cortez Ave, Descanso Ave, and the south end of Seacoast Dr. Typical pathways for overwash include beach access points in the sand berm or sea wall, and over (through) riprap and other shore protection structures (Fig. 3).

Fig. 3
figure 3

Land elevations above MSL from ac a truck LiDAR survey on 2/7/2019 and d a 2016 airborne a lidar survey (Young et al. 2018). The high back beach is a curved sea wall at Palm Ave, b elevated, discontinuous sand berm at Descanso Ave, and c line of porous riprap at Cortez Ave

Regular beach elevation surveys have been collected at IB since 2008. Ludka et al. (2018, 2019) describe beach evolution and methodology for surveys collected through 2016. Subaerial surveys collected after 2016 use a truck-mounted LiDAR (RIEGL VMZ-2000) scanner, coupled with a geodetic grade dual-frequency Trimble BD982 global navigation satellite system and a Trimble AP20 inertial measurement unit. Point clouds are filtered and edited to remove erroneous points (e.g. birds and people) and processed into 1-m2 resolution digital elevation models used to quantify vertical beach change. LiDAR-derived elevation maps highlight the low-lying topography at persistent IB flood sites (Fig. 3).

Between 2008 and 2012, winter sand levels at the riprap at Cortez Ave were typically between 2.5 and 3.5 m (Fig. 4a). Minimum levels are unknown. A 2012 nourishment with coarser than native sand elevated subaerial sand volumes for several years, including when exposed to the energetic winter waves of the 2015–16 El Niño, followed by a stormy 2016–17 winter (Fig. 4b). Recent winter sand levels (Fig. 4c, d) are again eroded to pre-nourishment levels, with > 1 m seasonal fluctuations (Fig. 4a). Interest here is on conditions in winter 2018/19, 2019/20 when we observed flooding by overtopping.

Fig. 4
figure 4

Winter (December through February) subaerial profiles at Cortez Ave. Vertical sand elevation versus cross-shore distance a winters 2008–12: pre-nourishment, berm crest 2.0–2.8 m, b winters 2012–18: nourishment influenced, berm crest ~ 3.2–3.9 m. Recent profiles c winter 2018/19 and d winter 2019/20 return to the 3 m pre-nourishment crest height in winter, and the summer profile accretes following the seasonal cycle (a single end of summer accreted profile is shown for comparison). Porous riprap (elevated dark mound) spans from sand level to above 5 m. Sand-laden overwash is observed flowing through, as well as over, the porous riprap

To validate incident wave nowcasts and forecasts, the Coastal Data Information Program (CDIP) station 155 (21-m water depth, 3.4-km from shore, Fig. 1) was instrumented with a Datawell Waverider MkIII directional wave buoy. Wave conditions were recorded continuously and reported every 30 min in real time from 24 August 2018 to 20 March 2020. In addition, pressure and near-bottom current (PUV) measurements were made in ~ 6-m water depth offshore of Cortez Ave using Nortek Vector ADVs and Paroscientific Digiquartz pressure sensors sampling at 2 Hz, deployed from 27 November 2018 to 22 April 2019, and 18 November 2019 to 27 May 2020.

Overtopping landward of the berm crest was inferred using 2 Hz pressure time series from RBRsolo3 D pressure sensors (referred to below as street pressure sensors) deployed during both winters near Cortez Ave, Descanso Ave, and the southern end of Seacoast Dr (Fig. 2). The street pressure sensors were deployed from 17 December 2018 to 5 March 2020, and 6 December 2019 to 13 March 2020, with gaps during 10–15 January 2019 and 3–16 January 2020.

Hourly water level and tidal predictions are specified using data from the La Jolla tide gauge, maintained by the NOAA Center for Operational Oceanographic Products and Services (CO-OPS, https://tidesandcurrents.noaa.gov/products.html). To confirm that water levels at La Jolla (located 35 km north of IB) are presentative of IB, a downlooking Xylem Nile502 radar water level sensor was deployed at the end of IB pier in 8-m water depth from 25 October 2019 to 4 February 2020. The CO-OPS datums database is used to define MSL at Imperial Beach as 0.77 m above the North American Vertical Datum of 1988 (NAVD 88) for the 1983–2001 tidal epoch.

3 Components of a flood forecast

Total water level (relative to MSL) is the superposition of offshore ocean water level and the wave effects,

$$ {\rm TWL} = \eta_{A} + \eta_{\rm NTR} + R_{2} , $$
(1)

where \(\eta_{A}\) is the astronomical tide, \(\eta_{\rm NTR}\) is the non-tidal residual (associated primarily with atmospheric forcing), and \(R_{2}\) is the 2% exceedance vertical level of wave runup. The ocean water level, the sum of \(\eta_{A}\) and \(\eta_{\rm NTR}\), is the TWL component unrelated to nearshore waves.

3.1 Ocean water level

Water levels measured with the IB pier radar and the La Jolla tide gauge are very similar (correlation = 0.998, root mean square error (RMSE) = 0.03 m) with no discernible time difference. Although the San Diego Bay tide gauge is closer to IB than to La Jolla, tides measured in San Diego Bay differ noticeably in amplitude and phase from the open coast IB and La Jolla.

The NOAA tidal prediction at La Jolla is \(\eta_{A}\) and \(\eta_{\rm NTR}\) is the difference between measured hourly-averaged water level at La Jolla (referenced to MSL) and \(\eta_{A}\). The tide dominates ocean water level fluctuations at IB, with the standard deviations of \(\eta_{A}\) and \(\eta_{\rm NTR}\) from the pier radar equal to 0.49 m and 0.03 m, respectively. The \(\eta_{\rm NTR} \) at the time of the forecast is assumed to persist for the full 6-day forecast period. The RMSE of ocean water level predicted at IB using the La Jolla tide gauge introduces normally distributed errors ranging from 0.03 m for nowcasts to 0.08 m at 6-days lead time (Fig. 5).

Fig. 5
figure 5

Ocean water level RMSE, computed as the RMS of the difference between the IB pier water level and the La Jolla-based forecast for the sample period Dec. 1, 2019–Feb. 4, 2020, versus forecast lead time

3.2 Runup model

The wave runup contribution to TWL, \(R_{2}\), is defined as the elevation exceeded by 2% of uprushing waves. Assuming Gaussian runup statistics, the 2% exceedance level is given by

$$ R_{2} = \overline{\eta } + \frac{S}{2}\,\,{\rm where}\,\,S = \left( {S_{SS}^{2} + S_{IG}^{2} } \right)^{\frac{1}{2}} $$
(2)

where \(\overline{\eta }\) is the super-elevation of water level due to wave breaking (wave setup), and \({S_{ss}} \) and \(S_{IG}\) are significant sea-swell (4–20 s) and infragravity band (20–250 s) swash heights. Many empirical runup parameterizations have been developed to estimate \(R_{2}\) (see review by Gomes da Silva et al. 2020). The most widely applied is Stockdon et al. (2006), where empirical fits define setup and runup components in terms of the incident wave height (\(H_{o} ), \) the deep-water wavelength (\(L_{o} )\), and the foreshore beach slope (\(\beta_{f} )\)

$$ R_{2} = 1.1\left( {0.35\beta_{f} \left( {H_{o} L_{o} } \right)^{\frac{1}{2}} + \frac{{\left[ {H_{o} L_{o} \left( {0.563\beta_{f}^{2} + 0.004} \right)} \right]^{\frac{1}{2}} }}{2}} \right). $$
(3)

The \(H_{0} L_{0}\) parameterizations follow from Hunt (1959), who proposed that the ratio of sea-swell runup to incident wave height scales with the Iribarren number (Battjes 1974), a dimensionless measure of dynamic beach steepness. Values of free parameters in (3) were fit to measurements from six beaches with varying morphology under a range of wave conditions. Other studies have used surfzone parameters to define runup (e.g. Stockdon et al 2006; Cohn and Ruggiero 2016), and the average period or frequency, rather than the peak (Atkinson et al. 2017, O’Grady et al. 2019; Dodet et al. 2019; Gomes da Silva et al. 2020, and references therein). These formulae usually approximate the beach profile with \(\beta_{f}\) and characterize the incident wave spectrum with \(H_{o}\) and \(L_{o}\).

With shoreward propagating waves specified at an offshore boundary, and with a fixed depth profile, runup in narrow laboratory wave channels can be accurately simulated with numerical models based on the Navier–Stokes equations (Lara et al. 2011; Torres-Freyermouth et al. 2019). However, the computation times are long and may be needlessly complex. Phase-resolving Boussinesq (Lynett et al. 2002) and non-hydrostatic models (Zijlema et al. 2011; Tonelli and Petti 2012; Tissier et al. 2012; Roeber and Cheung 2012; Smit et al. 2014) are a viable compromise between computational effort and accuracy. Bouss1D, BOSZ, and COULWAVE are nonlinear Boussinesq-type models that have been used to investigate wave and runup processes on beaches and coral reefs (Lynett et al. 2002; Yao et al. 2020; Roeber et al. 2010, 2012; Pinault et al. 2020 and many others). The SWASH (Simulating WAves till SHore, Zijlema et al. 2011) model has accurately simulated nonlinear surfzone wave evolution (Smit et al. 2014, Rijnsdorp et al. 2015), wave runup (Ruju et al. 2019) and wave overtopping (Suzuki et al. 2017) in the laboratory. Verification with field observations of runup for this class of models is more limited, but promising. Given the depth profile and in situ observations of shoreward propagating waves at the offshore boundary, the surfzone wave transformation and the resulting runup are well simulated by the 1D SWASH model (Fiedler et al. 2018).

Fiedler et al. (2020) (hereafter F20) showed the utility of runup formulae optimized for a given location based on SWASH runup simulations using typical storm waves and beach bathymetry with variable foreshore slopes. Simulations are used to calibrate a relatively simple, empirical parameterization. A notable difference from the \(H_{o} L_{o}\) dependence of S06 is the use of frequency-weighted integrals of sea-swell wave frequency spectra, \(E\left( f \right),\) which accounts for broad-banded and multi-peaked wave fields. The F20 approach, called the IPA (Integrated Power law Approximation), uses SWASH to simulate extreme storm wave events on a representative eroded beach profile determined from historical surveys. SWASH simulations are computed for many storm wave conditions identified in the regional Monitoring and Prediction (MOP) system (O’Reilly et al. 2016) wave hindcasts (Sect. 3.3). A general empirical power law form for the energy of the wave-driven components of TWL at the shoreline is

$$ \left[ {E_{IG,SS,\eta } } \right]_{\rm shoreline} = \left[ {\alpha \beta_{f}^{l} \mathop \int \limits_{SS}^{{}} E\left( f \right)^{m} f^{n} {\rm d}f} \right]_{\rm deep} $$
(4)

where \(\alpha\) is dimensional, \(E\left( f \right)\) is the incident wave spectrum at frequency \(f\), and (l,m,n) are powers determined for each wave-driven component. F20 used the MOP \(E\left( f \right)\) and SWASH runup to parameterize the components of \(R_{2}\), and found IB conditions at Cortez Ave, yield best-fits of

$$ \overline{\eta } = 0.21\mathop \int \limits_{SS}^{{}} E\left( f \right)^{0.45} f^{ - 1} {\rm d}f $$
(5)
$$ S_{ss}^{2} = 16*0.99\beta_{f}^{2} \mathop \int \limits_{SS}^{{}} E\left( f \right)^{0.45} f^{ - 1.85} {\rm d}f $$
(6)
$$ S_{IG}^{2} = 16*0.15\mathop \int \limits_{SS}^{{}} E\left( f \right)^{0.9} f^{ - 0.65} {\rm d}f $$
(7)

where integrals are evaluated over the incident sea-swell frequency band. Equations 57 are a model of the SWASH model on a typical winter profile with variable foreshore slopes. The left-hand side units are \(m^{2}\) in Eq. 5 and 6, and \(m\) in Eq. 4. The units are the same for the right-hand side, where f is in Hz and \(E\left( f \right)\) is \(m^{2} Hz^{ - 1}\). S06 include a \(\beta_{f}\) dependence for \(\overline{\eta }\) as well as \(S_{SS}\), whereas for the IB conditions, only the IPA \(S_{SS}\) depends on beach slope. The classic Hunt (1959) parametrization uses \(H_{o} L_{o}\), which scales as \(E^{1/2} f^{ - 2}\), most similar to \(S_{SS}\) (Eq. 5). The scaling consistency of IPA with Hunt (1959) is affirming, as Hunt (1959) examined runup on steep slopes where sea-swell may have dominated runup. The negative frequency exponents in \(\overline{\eta }\) and \(S_{IG}\) are smaller than the − 1.85 of \(S_{SS}\), indicating a weaker inverse dependence on wave frequency. (When the exponent is 0, there is no frequency dependence.)

Using SWASH model simulations on historical beach profiles, F20 show that IPA yields a more accurate runup hindcast than S06, even if S06 is tuned to the SWASH database. Note that IPA Eq. 6 is derived with a known foreshore \(\beta_{f}\). Scatter between SWASH and IPA arises from uncertainty in (submerged and subaerial) \(h\left( x \right)\), and in the offshore IG boundary condition for a given sea-swell spectrum. For later uncertainty assessments, note that the IPA-SWASH \(R_{2} \) emulation of F20 for a known beach slope has RMSD scatter ~ 9% of \(R_{2}\). We have increased this to 20% to account for h(x) uncertainty accompanying a larger training dataset (Fig. 6), as errors would be larger if the beach bathymetry were much different than observed in the F20 historical profiles.

Fig. 6
figure 6

IPA error (difference between IPA and SWASH \(R_{2}\)) versus Swash model \(R_{2}\). IPA is a model of the SWASH model. Dots are the 439 storms and 2 boundary conditions used to calibrate IPA on a representative, eroded winter depth profile (Fiedler et al. 2020). Scatter results from two different offshore IG boundary conditions, and the simplistic description of \(h\left( x \right)\) with a single foreshore slope in IPA. Errors in wave nowcasts and forecasts, and ocean water level, are not included. The RMSE noise model (gray, 0.2 \(R_{2}\)) is used in Monte Carlo simulations. Additional simulations with a range of eroded profiles (not shown) demonstrate that forecast model errors can double if the subaqueous profile evolves

3.3 Incident waves

The CDIP MOP system provides the incident wave spectrum \(E\left( f \right) \) in 10-m depth at locations spaced 100-m longshore. MOP hindcasts of \(E\left( f \right)\) are based on observations from a network of wave buoys. MOP forecasts up to 6 days in advance are specified using the NOAA WaveWatch III model forecasts at the buoy locations (O’Reilly et al. (2016)). The relatively strong dependence on wave frequency in \(\overline{\eta }\) and \(S_{IG}\) (Eqs. 5, 7) increases the importance of swell to \(R_{2}\) relative to high frequency seas. The warning system benefits from the relative predictability (a few days in advance) of the swell that drives extreme wave runup. Errors for individual wave events can be substantial. Overall, the buoy-driven MOP model hindcasts have relatively low bias and therefore are better suited for quantifying mean (e.g., monthly or annual) nearshore wave climate conditions rather than extreme or individual wave events. Nevertheless, despite these limitations, the MOP model proved useful for flood warnings at IB.

Errors in runup components caused by errors in MOP incident wave spectrum are evaluated as the difference in \(\overline{\eta },\) \(S_{SS}\), and \(S_{IG}\) (Eqs. 46) estimated using the observed buoy and MOP \(E\left( f \right)\). Typical time series of each are depicted (Fig. 7, a–c). Errors are approximated as normal and increase with amplitude (Fig. 7d–f), with normalized RSMEs of 0.16, 0.12, and 0.17, for \(\overline{\eta },\) \(S_{SS}\), and \(S_{IG} ,\) respectively.

Fig. 7
figure 7

Estimated IPA runup components a setup b \(S_{SS}\)/2 and c \(S_{IG}\)/2 at Cortez Ave for 2019 using incident waves from the offshore wave buoy, the Cortez PUV, and MOP line #45 (see legend). All estimates are linearly shoaled to 21-m depth for comparison. d, e Corresponding scatter plots of difference (MOP-PUV) versus PUV. Note vertical scales differ

3.4 Foreshore beach slope

IPA sea-swell swash energy depends on foreshore beach slope (Eq. 5). Beach slope variability is examined using Cortez Ave elevation profiles above MLLW during Dec-Mar of 2018/19 and 2019/20, when nourishment influences were reduced (Fig. 4c, d). The foreshore profile varies considerably, as expected for an equilibrium-like beach response to energetic winter swell (Yates et al. 2009; Ludka et al. 2015). The time-averaged subaerial beach slope (\(\beta ) \) over these 14 winter surveys increases approximately linearly with elevation z above MSL

$$ \beta \left( z \right) = 0.029z + 0.042 $$
(8)

with RMS slope difference of the 14 profiles compared to the mean profile = 0.03 (Fig. 8). For IPA estimates (Eq. 6), we compute the spatially-averaged foreshore beach slope \(\beta_{f}\) following S06 using Eq. 8, and assume the RMSE for \( \beta_{f} = 0.03\). The increase in slope with elevation (i.e., concave foreshore beach) causes \(\beta_{f}\), and hence \(S_{SS}\), to increase with tidal height for approximately steady wave conditions. The effect of the 0.03 RMSE \( \beta_{f} \) uncertainty on \(S_{ss}^{2}\) (Eq. 5) is substantial, but \(S_{IG}^{2}\) and \(\overline{\eta }\) are unaffected. Extension of IPA to other seasons would require an evaluation of the IPA parameters (Eqs. 57) for non-winter wave and beach conditions, and specification of seasonal changes to \( \beta_{f}\).

Fig. 8
figure 8

Mean Cortez Ave beach slope versus elevation for eroded profiles in Fig. 4c, d. Slope increases approximately linearly with elevation

4 Total water level predictions

4.1 Hindcasts and flood thresholds

Overtopping thresholds, the TWL elevations where overtopping leads to minor to severe flooding, vary with location. At flood sites, access ramps, riprap, and protective berms occasionally bulldozed by the city, create complicated pathways for overtopping flows (Fig. 3). The curved seawall at Palm Ave reflects incident waves and creates oblique and colliding crests. \(R_{2}\) represents a runup exceedance, so statistically an extreme swash uprush with low probability can exceed a given flood threshold, even if the TWL falls below threshold. A post-event survey might identify evidence of overwash associated with only one extreme swash wave. In general, and consistent with the definition of \(R_{2} , \) a threshold for minor flooding represents a high point elevation that is overtopped by 2% of runups, leading to ocean water impacting the otherwise dry landward terrain (Fig. 2, upper panels). Typically, there are several overtopping waves for each high tide flooding event. Streets are wet with some sandy overwash, and water ponding on Seacoast Dr. Moderate to severe flooding causes disruption of traffic on Seacoast Dr, and significant deposition of sand and debris on side streets leading to the beach, requiring bulldozers for removal (Fig. 3, lower panels).

To assess flood thresholds, we compiled anecdotal reports of flooding from the last 10 years of newspaper reports, photos posted on social media, and log reports maintained by the City of Imperial Beach, and computed the daily peak TWL hindcast during the flood event (Fig. 9). In winter 2018/19 and 2019/20, flooding was also documented by bore signals in the street pressure (Fig. 1a) records, and our own visual observations. Thresholds for minor flooding, based on a mean of the four lowest TWL values, are 1.9 m at Cortez Ave, and 2.2 m at Palm Ave and Descanso Ave (Fig. 9). Flooding was also reported in winter 2016/17 and 2017/2018, when a pad of nourishment sand (elevation > 4 m) reduced the direct impact of waves on riprap (Fig. 4b). However, groundwater seeped and pooled behind the pad on streets where the level was 1.6 m below the 4 m elevation of the pad (Ludka et al. 2018). The simplistic approach of a vertical threshold for flooding ignores the important role of infiltration and elevated groundwater on a nourished beach. Based on our observations from 2018 to 2020 and historical flood reports, we estimate the threshold for moderate to severe flooding from overtopping to be 1 m above the mild flooding threshold.

Fig. 9
figure 9

Flooding thresholds. Hindcasts of the daily peak ocean water level versus IPA \(R_{2}\) for historical flood events reported in newspapers, social media, and SIO observations at (left to right) Cortez Ave, Descanso Ave, and Palm Ave. High total water level (TWL) occurs when high \(R_{2}\) and ocean water level combine. TWL elevation thresholds for minor flooding, approximately the minimum TWL with reported flooding, are noted in each panel and shown in Fig. 3. The color of each dot represents the year of observation. Minor floods were not routinely monitored prior to winter 2018/2019. Flooding events between the fall 2012 and winter 2018/19 may be influenced by groundwater effects associated with a perched nourishment berm. Ocean water levels for 2010–2016 are higher than most of 2018–2020, suggesting that only larger flooding events were reported during 2010–2016

Mild flood thresholds from observations (~ 2 m) and elevation maps (Fig. 3) are similar. Lacking historical flood data, an overtopping flood threshold based solely on elevation maps is not unreasonable, unless rain or groundwater effects are important, or the berm elevation changes substantially. Ongoing flooding observations are key to developing accurate thresholds for a range of conditions.

On-site observers and the street pressure sensors are used to identify times of overtopping and flooding for comparison with TWL hindcasts (Fig. 10). The TWL hindcast captures 18 confirmed flood events in 2018–19, and 22 in 2019–20 (Fig. 10). Moreover, the days of TWL > 3 m correspond to moderate flooding at Cortez Ave based on on-site observations. Flooding on 18 Jan. 2019 (Fig. 2) was the highest of the two winters with TWL approaching 4 m, and \(R_{2} \) reaching 2.6 m. Incident waves were more energetic in 2018–19 than 2019–20, resulting in a 51% higher \(R_{2}\) variance.

Fig. 10
figure 10

Hindcasts of ocean water level (dark blue) and total water level at Cortez Ave based on IPA runup (light blue) for the winters of a 2018–2019 and b 2019–2020. The 1.9-m threshold for minor flooding is selected to capture flood events confirmed by an on-site observer, or street pressure sensor (legend). The street pressure sensor was not deployed prior to December 17th, 2018, nor from January 10th, 2019 to January 15th, 2019, and January 3rd, 2020 to January 16th, 2020

The TWL estimate captures most observed flood events; however, false reports do occur. Several events exceed the threshold without apparent overtopping (e.g., early February 2019), and a reported flood (early February 2020) corresponded to a TWL below threshold (Fig. 10). Minor flood levels predicted in early winter were unconfirmed by the street pressure sensor observations in December 2019 (observations were unavailable in December 2018), suggesting that our winter mean beach profile may need to be adjusted for the first sizeable winter wave events when the beach profile is in seasonal transition.

4.2 TWL uncertainties

These discrepancies highlight the need for TWL uncertainties for hindcasts and forecasts, not simply a binary measure of flooding based on threshold exceedance. Monte Carlo simulations are used to estimate the RMSE of TWL. Each TWL value is computed 1000 times, with independent random errors drawn for the ocean water level (Sect. 3.1), the IPA model (Sect. 3.2), the beach slope (Sect. 3.4), and the MOP wave input, accounting for forecast lead time (except for beach slope).

MOP wave input errors are evaluated as the difference in \(R_{2}\) components (\(\overline{\eta },\) \(S_{SS}\), and \(S_{IG}\), Eqs. 57) estimated using the observed buoy and MOP \(E\left( f \right)\). The component errors are used instead of an error estimate of \(R_{2} \) so as to include beach slope uncertainty in \(S_{SS}\) (Eq. 6) in the Monte Carlo analysis. The normalized RMSEs for \(\overline{\eta },\) \(S_{SS}\), and \(S_{IG}\) are listed for increasing forecast lead times (Table 1).

Table 1 Normalized RMSE for runup components computed from the differences in the estimated components using MOP and observed PUV wave spectra

4.3 Flood predictions and probabilities

TWL forecasts with 6-day lead time estimated using MOP and ocean water level forecasts (Sect. 3) are appended to hindcast time series in Fig. 11. Error bars (\(\pm 1\) RMSE) are shown for high tides when flooding is most likely to occur. In this example for Cortez Ave, the forecast was made at 7 pm PST 21 January 2020. Over the 6-day forecast, on 4 occasions TWL is predicted to exceed the minor threshold by at least one RMSE. On 25 January TWL is within one RMSE of the moderate threshold.

Fig. 11
figure 11

Hindcasts and forecasts of TWL at Cortez Ave. a The ocean water level (dark colors) and \(R_{2}\) (light colors) combine to form TWL. A forecast out to 6 days, made at 7:00 pm PST 21 January 2020, is appended to the hindcast elevation above MSL. Error bars at high tide are \(\pm 1\) RMSE. Red dots indicate time of a confirmed flood. MHHW (mean higher high water) is also shown. b The probability that the TWL will exceed the minor and moderate flood thresholds at high tides, computed from the TWL RMSE and assuming normally-distributed errors. (lower) A 60-cm increase in ocean water level more than doubles the probability of moderate events, all else held constant

Instead of tracking the TWL relative to thresholds, we mimic typical meteorological forecasts (e.g., 40% chance of rain) and estimate the probability that TWL will exceed each threshold. The probability of exceedance at each high tide (Fig. 11b) shows that at forecast time of 7 pm on 21 January 2020, the next 5 higher-high tides had at least a 60% chance of reaching minor flood levels, with 2 high tides exceeding 90% probability. Moderate flooding was not expected for this 6-day forecast period, although the high tide on 25 January had a slight chance of reaching that threshold. Measurements confirm that minor flooding did occur at the higher-high tides between 22 and 25 January, and the flooding on 25 January was the most severe of the week. Flooding did occur on 24 January when TWL did not reach the minor flood threshold (i.e., there was still a ~ 40% chance of flooding). Likewise, on 26 January, when the forecast was nearly 5 days out, the predicted TWL was above threshold with a ~ 70% chance of flooding, but flooding was not observed (i.e., the 30% chance that it would not flood was correct). A 0.6-m rise in ocean water level by 2100 is likely (Griggs et al. 2017). Holding all else constant, threshold flood levels will be exceeded more often, for longer time periods, and moderate flooding will become common (Fig. 11b, lower). For example, minor flooding would accompany all high tides during the January 2020 forecast window, with moderate flooding expected to occur during most spring high tides. Assuming fixed conditions, a 1-m rise in sea level would result in minor flooding during all high tides during the Dec-February winter months.

The relative contribution of each input variable to TWL error is assessed by removing the error of each variable separately in the Monte Carlo simulations of the records shown in Fig. 10, and computing the change in TWL error. The largest \(R_{2}\) errors at all forecast times are infragravity wave boundary conditions (included in IPA error estimates, Fig. 12). MOP errors increase significantly over time because 6-days is pushing the limits of swell predictability. Errors in offshore water level increase over time, but are always relatively small, followed by the assumption of a fixed beach profile (Fig. 12).

Fig. 12
figure 12

The estimated contribution to the mean square error (MSE) of hindcast TWL at Cortez Ave, computed over the two winters (Fig. 10). The bars represent the contribution to TWL MSE due to errors in ocean water level (blue), IPA model (red), and MOP wave input and slope (yellow). The component MSEs are computed as the difference in total MSE when all errors are propagated through to TWL, compared to when the error of the indicated component term is excluded. Errors are calculated for \(R_{2} \) > 1 m

5 Discussion

The Imperial Beach forecasts are optimized for probable flood conditions, coincident high tide and high waves on an eroded winter beach profile. However, flooding can occur at other tidal phases. For example, in early January 2020 (Fig. 10b), the TWL reached minor flood levels for an energetic wave event arriving at neap tides. Although flooding was not observed at this time, the event illustrates that tuning the model for only high tide is limiting. Similarly, high runup has occurred at IB on summer accreted profiles in response to long swell arriving from Southern Hemisphere storms. The IPA formalism is being developed for additional sites, with beach, wave, and water level conditions different from Imperial Beach. In the meantime, until the IPA model has been developed and tested for these conditions, we recommend using the S06 parameterization, if possible tuned for the site (F20). Without site-specific tuning, S06 errors can be large, especially with beach bathymetries or wave conditions much different from the original calibration data (Stockdon et al. 2014). Using simulated energetic conditions, F20 showed calibrated S06 outperformed uncalibrated S06 at IB, decreasing the RMSE by half (0.88 m vs. 0.44 m), while the IPA further decreased the RMSE to 0.20 m.

The IPA parameters are only as good as the SWASH simulations from which they were derived. F20 emphasize that uncertainties in the offshore bathymetry and offshore wave boundary condition can lead to large runup errors in SWASH that translate directly to the IPA. IB remains a challenging location to resolve offshore bathymetry and boundary conditions as poor water quality during winter limits in situ observations and surveys. For locations with limited offshore bathymetry, remotely sensed bathymetry is an attractive alternative. The one-dimensional SWASH model neglects alongshore variation in bathymetry and directional spreading in incident waves. Two-dimensional model tests are hampered by both the large spatial domains required to allow infragravity edge waves, and continued uncertainty in offshore infragravity wave boundary conditions. We did not resolve runup at the point of overtopping. The IPA assumes runup on a beach profile, and extrapolates that to heights that exceed physical land levels, at which point flooding is assumed to occur. Additional model testing is needed to establish overtopping characteristics over seawalls, protective riprap, and protective berms.

The IPA model incorporates a spatially variable foreshore beach slope, which imposes a tidal modulation on the sea-swell component of runup at IB. Our Monte Carlo analysis suggests that uncertainties associated with foreshore slope at IB appear to be secondary compared with the MOP forecast error and the IPA model uncertainty (Fig. 12); however, further model skill improvements might be attainable by using a beach morphology model that allows the beach to evolve throughout the winter swell season. F20 document how a changing foreshore beach profile (slope reduction by nearly 40%) during a single storm event influences the swell component of runup.

The specification of flood thresholds with observations is challenging. Street pressure sensors provide useful indicators of flooding occurrence, but a few point sensors cannot fully assess flood impact over a complicated landscape where the path of overwash bores can be blocked and rechanneled, in some cases due to protective barriers put in place by the city. For example, threshold levels at Descanso Ave were observed to vary depending on the height of a protective sand berm built by city bulldozers prior to large flood events, but a static threshold is used in the warning. Sensors that measure the cumulative impact of overtopping, such as a pressure sensor at the low-point in a flood zone, or imagery that can estimate the area, and with detailed topography the volume of a flood zone would be informative, but we found it difficult to make these measurements at IB. Video can capture the spatial extent of overtopping; however, these systems do not work well at night, and the highest tide of the day in winter at IB often occurs in the early morning dark. We have made progress using LiDAR observations, described in a subsequent study. In general, methods to observe overtopping volumes and flood volumes require further development, and future work to quantify the role of groundwater and drainage infrastructure in relation to flood volume will enhance the forecast.

6 Conclusions

An early warning system for wave-driven flooding at Imperial Beach has been developed using regional wave and water level observations, historical beach surveys, and a numerical runup model. Observations are used to calibrate and validate the approach. A crucial element for predicting TWL is the runup model, which transforms incident wave energy to the swash zone. As computing resources improve, real-time SWASH simulations could be incorporated into the warning system, thus reducing error in the forecasts associated with the IPA. Until then, the IPA provides a useful surrogate for SWASH.

Uncertainties associated with the incident sea-swell and infragravity wave energy, foreshore beach slope, offshore morphology, ocean water level, and the IPA model are propagated to TWL uncertainty, from which probabilistic flood forecasts are derived up to 6 days in advance. Flood thresholds, measured as distance above MSL, established using historical reports of flooding at the IB sites, are qualitatively consistent with digital elevation maps. The effect of active mitigation (i.e., time varying thresholds) is understood poorly, but should be incorporated into an operational flood system.

The warning system approach developed for IB is transferable to other sites given reasonable estimates of ocean water level, incident waves, and beach profiles. Observations of flooding severity are required to establish flood thresholds. As discussed by F20, errors in the offshore boundary condition and the offshore bathymetry can lead to large SWASH uncertainties. Error estimates of TWL are critical to prudent use of warnings in practice, and to accommodating a lack of information (e.g., limited offshore bathymetric profiles, uncertain incident wave conditions). As more information is obtained, errors decrease. Observations during two consecutive winters indicate that Cortez Ave was the most flood-prone area of IB, as is well known by IB inhabitants. This flooding “hot spot” occurs because an offshore shoal (Fig. 1) can focus wave energy at Cortez, a relatively narrow beach with low berm top elevations, and close proximity to Seacoast Dr (Fig. 3). Limited adaptation strategies (e.g., enhancement and raising of protective riprap, beach nourishment) will mitigate overwash and flooding of Cortez Ave and Seacoast Dr; however, beachfront structures remain vulnerable.

We initially expected the Datawell buoy and IB pier water level sensor would be needed long-term. Ultimately, the MOP wave predictions and the La Jolla tide gauge were adequate. In situ PUV observations proved valuable to validate the SWASH model, and to characterize offshore infragravity waves. Flood observations on land are the observations most lacking at IB and we anticipate at other sites. Water level sensors in critical locations along the oceanfront can gather time series of flooding, information key for establishing the flood thresholds crucial for graded alerts. The required flood measurements are site specific, reflecting varying degrees of shoreline protection, and are a key consideration for an early warning system in a new location.