André Doria

R T Guza
email: rguza@ucsd.edu

William C O'reilly

M L Yates

Observations and modeling of San Diego beaches during El Niño

KEY POINTS

1) Subaerial sand levels were observed at 5 southern California beaches for 16 years.

2) Cobbles and bedrock sometimes reduced the mobility of eroded shorelines.

3) Inclusion of site-specific geological boundaries improves the performance of an equilibrium model.

Abstract: Subaerial sand levels were observed at five southern California beaches for 16 years, including notable El Niños in 1997-98 and 2009-10. An existing, empirical shoreline equilibrium model, driven with wave conditions estimated using a regional buoy network, simulates well the seasonal changes in subaerial beach width (e.g. the cross-shore location of the MSL contour) during non-El Niño years, similar to previous results with a 5-year time series lacking an El Nino winter. The existing model correctly identifies the 1997-98 El Niño winter conditions as more erosive than 2009-10, but overestimates shoreline erosion during both El Niños. The good skill of the existing equilibrium models in typical conditions does not necessarily extrapolate to extreme erosion on these beaches where a few meters thick sand layer often overlies more resistant layers. The modest over-prediction of the 2009-10 El Niño is reduced by gradually decreasing the model mobility of highly eroded shorelines (simulating cobbles, kelp wrack, shell hash, or other stabilizing layers). Over prediction during the more severe 1997-98 El Niño is corrected by stopping model erosion when resilient surfaces (identified with aerial imagery) are reached. The trained model provides a computationally simple (e.g. nonlinear first order differential equation) representation of the observed relationship between incident waves and shoreline change.

Introduction

Coastal communities and beaches provide abundant ecological, recreational, and socio-economic wealth [Nicholls et al., 2007;[START_REF] Yang | Geospatial analysis of barrier island beach availability to tourists[END_REF][START_REF] Mclachlan | The Ecology of Sandy Shores[END_REF]. Increasing coastal populations [START_REF] Moore | Coastal erosion hazards in Santa Cruz and San Diego Counties, California[END_REF], long-term climate change [START_REF] Keeling | [END_REF][START_REF] Rahmstorf | Recent Climate Observations Compared to Projections[END_REF], polar ice melt [START_REF] Dyurgerov | Twentieth century climate change: Evidence from small glaciers[END_REF][START_REF] Bamber | Reassessment of the Potential Sea-Level Rise from a Collapse of the West Antarctic Ice Sheet[END_REF], and sea level rise (SLR) forecasts of between 0.8-2 m of SLR by 2100 have raised concerns about the long-term (e.g. centuries) fate of beaches, coastal infrastructure, and coastal cliff retreat [START_REF] Zhang | Global warming and coastal erosion[END_REF][START_REF] Pfeffer | Kinematic constraints on glacier contributions to 21st-century sea-level rise[END_REF][START_REF] Vermeer | Global sea level linked to global temperature[END_REF][START_REF] Gallien | Predicting tidal flooding of urbanized embayments: A modeling framework and data requirements[END_REF]. At shorter time scales, accelerated coastal erosion may be caused by decadal oscillations in the frequency, severity, and tracks of storms [START_REF] Graham | Evidence for intensification of North Pacific winter cyclones since 1948[END_REF]Allan and[START_REF] Allan | Climate controls on US West Coast erosion processes[END_REF]2002;Ruggiero et al., 2010a]. California, Oregon, and Washington beaches suffered severe erosion from the intense and frequent storms during the El Niños of 1997-98 and 2009-10 [Revell et al., 2002, 2011;;[START_REF] Barnard | The impact of the 2009-10 El Niño Modoki on U.S. West Coast beaches[END_REF].

Effectively managing beaches now, and in a future with potentially altered wave climates, requires quantifying the relationship between beach change and waves.

However, testing of shoreline change models on the U.S. West coast has been limited.

Genres of shoreline models include process-based and empirical. Process models [e.g. SBEACH, [START_REF] Larson | SBEACH: Numerical model for simulating storminduced beach change[END_REF]XBeach, Roelvink et al., 2009;and CSHORE, Johnson et al., 2012] necessarily parameterize the complex physics of sediment transport with combined steady and oscillatory flows. Empirical models based on an equilibrium hypothesis tune "bulk response" parameters, and have skill in simulating observations of shoreline change on time scales of months to a few years [START_REF] Miller | A simple new shoreline change model[END_REF]Yates et al., 2009a;[START_REF] Davidson | Forecasting seasonal to multi-year shoreline change[END_REF], Ruggiero et al., 2010b[START_REF] Davidson | A simple equilibrium model for predicting shoreline change[END_REF][START_REF] Splinter | A generalized equilibrium model for predicting daily to interannual shoreline response[END_REF]. Equilibrium beach models quantify the hypotheses [START_REF] Wright | Short-term changes in the morphodynamic states of beaches and surf zones: An empirical predictive model[END_REF] that: (a) for a constant wave field, there is an equilibrium beach morphology (the equilibrium beach) that would remain constant in time, neither eroding or accreting, (b) a beach in disequilibrium with the ambient waves changes towards the equilibrium shape, and (c) the change rate is proportional to the disequilibrium. [START_REF] Miller | A simple new shoreline change model[END_REF] applied equilibrium concepts to derive (1) where S is the shoreline location (defined as the cross-shore position of a shallow depth contour, here Mean Sea Level (MSL)), is the beach disequilibrium, and the empirical k depends on wave energy, grain size, and other local factors. Yates et al.

[2009a] (hereafter Y09) showed that an equilibrium shoreline model had skill at three southern California beaches over five years (2004)(2005)(2006)(2007)(2008)(2009). [START_REF] Ludka | Field Evidence of Beach Profile Evolution Towards Equilibrium[END_REF] recently developed an equilibrium beach profile model using up to 10 years of observations that included the 2010 El Niño. Here, the southern California observations of previous studies [START_REF] Shepard | How much data is enough? The importance of morphological sampling interval and duration for calibration of empirical shoreline models[END_REF][START_REF] Winant | Description of seasonal beach changes using empirical eigenfunctions[END_REF][START_REF] Nordstrom | Sand level changes on Torrey Pines Beach, California[END_REF][START_REF] Flick | Performance documentation of the Longard Tube at Del Mar, California 1980-1983[END_REF], Yates et al., 2009a, 2009b, 2009c] are expanded to include additional sources spanning up to 16 years (1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014), including the more severe 1997-98 El Niño winter.

The Y09 shoreline model is extended by gradually decreasing the model mobility of highly eroded shorelines (coarsely accounting for cobbles and other natural armoring), and stopping erosion when a non-erodible layers (e.g. bedrock) is reached.

First, the beach sites (Section 2), and wave and sand level observations (Section 3) are described. In Section 4, observations of waves and shoreline (MSL contour) location are used to tune an equilibrium-type shoreline model. Results are discussed in Section 5, and summarized in Section 6.

Beach Sites

In southern California, wave conditions and beach sand levels vary seasonally [START_REF] Shepard | How much data is enough? The importance of morphological sampling interval and duration for calibration of empirical shoreline models[END_REF][START_REF] Winant | Description of seasonal beach changes using empirical eigenfunctions[END_REF][START_REF] Nordstrom | Sand level changes on Torrey Pines Beach, California[END_REF][START_REF] Flick | Performance documentation of the Longard Tube at Del Mar, California 1980-1983[END_REF], Yates et al., 2009a, 2009b, 2009c]. Sand elevations were measured at five San Diego County beaches (from south to north, Figure 1): Imperial Beach (4 km alongshore span), Torrey Pines (8 km), Solana Beach (2.6 km), Cardiff (2 km), and Camp Pendleton (2.5 km). Median sand sizes range between 0.15-0.28 mm (Table 1), and beach slope between 0.01-0.08 (Table 1).

Imperial Beach (Figure 1b) contains a recreational pier, two short groynes in the northern 300 m, and the Tijuana River mouth at the southern end. Most of the beach is backed by low-lying urban development and protective riprap, seawalls, and cobble berms (Figure 2). The southern 6.5 km of Torrey Pines State Beach (Figure 1c), is backed by 50-110 m high-relief sandstone cliffs, and the northern 1.5 km is fringed by riprap and the Los Peñasquitos Lagoon inlet [START_REF] Moore | Coastal erosion hazards in Santa Cruz and San Diego Counties, California[END_REF][START_REF] Young | Coarse sediment yields from seacliff erosion in the Oceanside Littoral Cell[END_REF]. Solana Beach (Figure 1d) is backed by 25 m sandstone cliffs [START_REF] Young | Coarse sediment yields from seacliff erosion in the Oceanside Littoral Cell[END_REF]

Observations

Sand Level Surveys

Surveys of subaerial beach sand levels from 1997-2014 at 5 beaches were obtained from several sources (Figure 1) including ( 1) cross-shore transects surveyed biannually from the back beach to ~8-10 m depth beginning in 1997 (San Diego Association of Governments (SANDAG); red transects in Figure 1) and (2) quarterly transects, beginning in 2004 (SIO; dense black, blue, or white transects in Figure 1) [Yates et al., 2009a].

(3) Monthly subaerial shoreline parallel surveys beginning at Torrey

Pines, and subsequently expanded to four additional sites (Imperial, Cardiff, Solana, and Camp Pendleton). ( 4) Airborne lidar in April 1998 (NASA's airborne topographic mapper (ATM); [START_REF] Brock | Basis and Methods of NASA Airborne Topographic Mapper Lidar Surveys for Coastal Studies[END_REF]) and biannually from May 2002 until October 2010 (Univ. of Texas, [START_REF] Yates | A Technique for Eliminating Water Returns from Lidar Beach Elevation Surveys[END_REF]). Lidar returns were removed offshore of the waterline location, estimated using water levels from a nearby tide gauge and runup approximated using local wave conditions [START_REF] Yates | A Technique for Eliminating Water Returns from Lidar Beach Elevation Surveys[END_REF]. Lidar sand levels were gridded onto 4 m 2 cells, using the cell median elevation to reduce the influence of outliers. Point density in the 1998 NASA lidar survey was low (0.57 points m -2 ), compared with the post-2001 biannual lidar surveys (~2 points m -2 ) [START_REF] Brock | Basis and Methods of NASA Airborne Topographic Mapper Lidar Surveys for Coastal Studies[END_REF][START_REF] Yates | A Technique for Eliminating Water Returns from Lidar Beach Elevation Surveys[END_REF]. Grid cells with less than 3 data points were discarded from the post-2001 lidar surveys. All data was necessarily retained in the lower density 1998 survey.

Surveys from different sources at the same approximate time and beach usually agree, with differences owing to variable amounts of spatial averaging (Figure 3).

Responding to seasonal variations in wave energy, the observed shoreline (e.g. MSL contour) locations usually varied seasonally by 25-30 m at all 5 study beaches (Figure 4; [START_REF] Winant | Description of seasonal beach changes using empirical eigenfunctions[END_REF]Yates et al., 2009b]). During the 1998 El Niño, shoreline retreat was maximal, about 25 m landward of the typical (e.g. 2004-2012) winter shoreline (Figure 4). Recovery from 1997-98 took several years, even with nourishments both shortly before (1997, Imperial Beach, 178,000 m 3 ) and after (1999, Solana Beach, 41,000 m 3 ) El Niño; however, during fall 1997, existing beach sand levels at several sites were historically lower than post-summer level observed in most other years. Accordingly, the erosive change during the 1997-98 El Niño was limited because of low sand levels preceding the event. Recovery following the less erosive 2009-10 El Niño was more rapid, effectively one season (Figures 3 and4). Spring-summer 2001 nourishments at Imperial Beach, Torrey Pines, Solana Beach, and Cardiff elevated sand levels to new maxima (Figure 4). The nourishment was detectable for about two years at Torrey Pines, either as a wider subaerial beach, or as an enhanced offshore winter sand bar [Yates et al., 2009c]. SANDAG winter surveys occur in spring and fall. The spring surveys usually occur after the winter erosion maximum in February-March (compare squares and circles in Figure 3a, in 2005-2008 inclusive), so the 1998 survey may not have captured the maximum erosion.

Waves

Waves typically approach the Southern California Bight from N-NW in winter and from S-SW in summer, and vary alongshore owing to sheltering by the Channel Islands and refraction over complex offshore bathymetry [Pawka,1983]. Local (e.g. < 30 m depth) bathymetric variations further refract and focus waves with appreciable alongshore energy variations over several hundreds meters alongshore. Directional wave buoys (CDIP, http://cdip.ucsd.edu; Figure 1a) initialized a spectral refraction model [START_REF] O'reilly | Comparison of spectral refraction and refractiondiffraction wave models[END_REF][START_REF] O'reilly | A comparison of spectral wave models in the Southern California Bight[END_REF][START_REF] O'reilly | Assimilating coastal wave observations into regional swell predictions. Part I: Inverse methods[END_REF]] that provided hourly wave estimates at 10 m depth every 100 m alongshore. Near-shore buoy deployments confirmed reasonably good model accuracy in relatively shallow water (20-30 m depth) at several of the study sites [START_REF] Young | Crossshore decay of cliff top ground motions driven by local ocean swell and infragravity waves[END_REF].

Waves were most energetic during strong El Niño winters (Figure 5a). For example, at the Oceanside buoy (Figure 1b), the hours of significant wave height H s exceeding 3 m were between 0-26 hours during 13 non-El Niño winters, compared with (3)

where a 0 and a 1 are empirically determined equilibrium wave energy coefficients. New here, mean-square error (RMSE) using surrogate management framework (SMF) optimization [START_REF] Booker | A rigorous framework for optimization of expensive functions by surrogates[END_REF][START_REF] Marsden | Optimal aerocoustic shape design using the surrogate management framework[END_REF].

S

Model-Data Comparison

Shorelines were hindcast for up to 16 years using the wave-driven equilibrium model, initialized with the earliest survey data point (typically fall 1997). Model calibration with a period including an El Niño yielded improved model-data agreement during both El Niño and non-El Niño years, and calibration with 2003-2011 is shown (Figure 6). The average model skill at Solana, Imperial and Torrey Pines beaches are between 0.55-0.60 (Table 3). At Cardiff and Camp Pendleton, two shorter beaches with river or lagoon mouths, skill was often less than 0.5. Two of the four modeled sections at

Cardiff have low skill (e.g. 0.22 and 0.41) and are located near a persistent lagoon mouth or a large bedrock platform extending from the subearial beach to wading depths. Camp

Pendleton was observed for the shortest time, and has the lowest R 2 (less than 0.5 at all modeled sections; Figure 6e), possibly resulting from the adjacent river mouth. Skill at all modeled locations was significant at the 95% level.

The model back beach erosive limit S bb (Figure 2; dashed horizontal line in

Discussion

Parameter Values, Response Times, and Initialization

Optimal model free parameters varied within and between sites (Table 3). Model error is weakly sensitive to the free parameter values, with only a 10% increase in model error for factor of two of changes in parameters (comparable to the differences between sites). Free parameter values surely depend on sediment availability, grain size, and possibly other environmental factors, but are only loosely constrained by the observations.

The best-fit shoreline adjustment time scales , averaged over each site, varied between roughly 10-20 days for erosion (with m), and the accretion spanned 29-64 days (with m; Table 3). Hypothetical initial conditions illustrate the rapid return (weeks to several months) of the model to equilibrium from artificially large disequilibria (crosses and triangular markers in Figure 7). Six rather different initial conditions in 1996, 1997, and 1998 all result in the same modeled shoreline by summer 1998 (grey curve in Figure 7). Model shorelines recovered from strong 1997-98 El Niño erosion by the following winter, more rapid than the observed multi-year recovery, demonstrating the model's failure to properly replicate the slow return of sand evidently displaced further offshore during the strong event (Figures 6b-d). Accretion is crudely parameterized in the model and requires future study.

Calibration Period

At Torrey Pines, Y09 found a relative 1.9 m increase in model-data RMSE during predictive model periods compared to the calibration period RMSE. [START_REF] Davidson | A simple equilibrium model for predicting shoreline change[END_REF] provide a more extensive calibration and validation discussion of a similar equilibrium- ) compared to this study). These differences may be partially attributed to the increased calibration period, as longer tuning generally resulted in different free parameters and a reduction in RMSE [Y09].

However, direct comparison to the Y09 results is cautioned, as modeled sections at the same beach are not necessarily identical to this study.

Additionally, the statistical nature of the calibration technique creates inherent variation to the resulting coefficients, as several solutions in parameter-space may produce similar results. The multiplicative nature of the model terms (2) also allows for changes in one coefficient to be compensated for by another coefficient.

Alongshore-averaged model coefficients provide a broad representation of the site-specific free-parameter value for bulk comparison to Y09 (Table 3). Alongshore averaged, had the greatest disparity (more than double in magnitude) relative to Y09

5-year calibrated at Torrey Pines. However, as noted previously, model skill is fairly insensitive to parameter values, with being the least sensitive parameter [Y09].

Fundamentally, model coefficients are weakly constrained by observations and differences between studies, even at similar beaches, are not necessarily remarkable. [START_REF] Davidson | A simple equilibrium model for predicting shoreline change[END_REF] and [START_REF] Splinter | A generalized equilibrium model for predicting daily to interannual shoreline response[END_REF] use an equilibrium model with forcing governed by wave power (rather than wave energy, E , in (2)) and the Dean parameter, which depends on grain size. The range of sand grain sizes is not taken into account here, and is relatively small (4 of the 5 beaches have D 50 between 0.15-0.18mm, (Table 1)). At Torrey Pines, Y09 showed replacing wave energy, E , in their shoreline model with H s or radiation stress S xx resulted in similar model skill, because E , H s , and S xx are strongly mutually correlated. [START_REF] Davidson | A simple equilibrium model for predicting shoreline change[END_REF] and others use an equilibrium condition based on the weighted average of antecedent waves, rather than on the present beach state. However, the present beach state depends on the previous wave conditions, and for the idealized case of a step change in time to a constant wave forcing, the equilibrium conditions of [START_REF] Davidson | A simple equilibrium model for predicting shoreline change[END_REF] and Y09 yield identical results. These different equilibrium models were also shown to yield similar results for the field observations [START_REF] Castelle | Equilibrium shoreline modelling of a high-energy meso-macrotidal multiple-barred beach[END_REF].

Alternative Model Formulations

The basic equilibrium equation of the present model ( 2), with a linear dependence of dS dt on the present wave energy E , and 4 free parameters, is referred to as the linear 4 model (the subscript specifies the number of free parameters). Additional alternative models are linear 3 , exp 4 , and cubic 4 . The linear 3 model reduces the number of free parameters to three by replacing with single valued C in (2a), following Yates et al., [2011]. The exp 4 and cubic 4 alternative models also simplify with C in (2a), but use more complex forms of E eq , for exp 4 (4) and for cubic 4 .

(5)

The model parameters S eq and E eq , and the rate of change dS dt and the response time , are similar in the range of common S and H s , while differing at the extremes (Figure 10). All models use the same erosion limiter S bb .

Overall (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011) the alternative models perform similarly, with typically small (<15%) improvements in model error relative to the Y09 model, which has no erosion limiter (Figure 9c). Model performance varied by site, but explained more than 50% of the variance over 16 years at most of the sandy beaches, similar to Y09 five-year hindcasts. The models differ from the Y09 model most significantly for extreme conditions only briefly encountered. While the Y09 model correctly identifies the 1997-98 waves as more erosive than 2009-10, it overestimates shoreline erosion during both El Niño events.

The cubic 4 model provided the greatest improvements in model skill (relative to Y09), with improved predictions for El Niño 2009-10 at beaches both where the erosion limiter was and was not reached (Solana Beach and Torrey Pines, respectively, Figure 9a,b). The over-prediction of the winter 2009-10 shoreline erosion ( , Figure 9d) was reduced using the cubic 4 model at all sites except Camp Pendleton, where over-prediction persisted. Model-data comparison at Camp Pendleton was generally poor irrespective of which model was used, perhaps owing to the close proximity of a river mouth. Typical reductions are about 5 m (up to 18 m peak reduction) relative to Y09. With large waves (H s = 4 m) and a heavily eroded shoreline (solid curves, S = -40 m, Figure 10d), dS dt for exp 4 and cubic 4 are much smaller in magnitude than for linear 3 (a simplified version of the Y09 model). Physical explanations for the reduced mobility of eroded beach face include the exposure of resistant strata and/or a reduction of the effective wave energy reaching the beach face owing to well-developed offshore sandbars.

Conclusion

Sixteen years of shoreline and wave observations, including two El Niños, 1997- boundaries (identified with aerial images), and using alternative, nonlinear forms of E eq (e.g. cubic 4 ) that gradually decrease the mobility of highly eroded shorelines (simulating cobbles, kelp wrack, enhanced offshore sand bars, and other stabilizing effects).

The shoreline location depends on complex processes occurring over the crossshore beach profile, and in some cases on adjacent profiles. Even significantly different equilibrium shoreline models often have similar skill [START_REF] Castelle | Equilibrium shoreline modelling of a high-energy meso-macrotidal multiple-barred beach[END_REF], which is also true for existing, more computationally demanding, physical process models for shoreline change. Application of any model to extreme conditions on sand-limited beaches with unknown substrates will requires site and condition specific calibration. Once trained, the present model provides a computationally simple (e.g. nonlinear first order differential equation) representation of the observed relationship between incident waves and shoreline change, including the effect of erosion resistant substrates.

  40 and 51 hours in the 1997-98 and 2009-10 El Niño winters, respectively. Total hours of H s between 2-3 m during the 1997-98 El Niño winter (more than 400 hours) dwarfed all other winters, nearly doubling those found in the second most energetic winter (e.g. 2009-10 winter; 220 hours of H s = 2-3 m; Figure5a). In 1997-98, H s exceeded 2 m for nearly 60 continuous hours, with frequent and prolonged sequences of energetic waves in early December 1997 and February 1998 (Figure5b). January 2010 had the longest period (~ 140 hours) of continuous H s exceeding 2 m (Figure5e).4. Shoreline Modeling4.1. Equilibrium Shoreline ModelAn existing equilibrium shoreline model [Y09] was modified to improve predictions during El Niños and other severe erosion conditions by accounting for durable limits (e.g. bedrock, seawalls, hard cliffs). The model assumes these relatively resilient boundaries were not eroded during the modeling period, and neglects cliff erosion, which would both relocate the back beach boundary and supply new sand to the beach. The comparative beach profile effects between armored and exposed back beaches are not included in the present model. With the shoreline location S defined as the cross- shore location of the MSL contour, the shoreline change rate dS dt depends on the present shoreline position S and incident wave energy E , equilibrium wave energy, is the wave energy for a given S that would cause no shoreline change . For the few occasions when highly accreted shoreline positions yielded non-physical negative

  photographic and lidar surveys. Shoreline retreat stops (e.g.) when . A

Figure

  Figure 6a-d) was reached during the 1997-98 El Niño (except Camp Pendleton), and

based 1 -

 1 D shoreline model. Both Y09 and[START_REF] Davidson | A simple equilibrium model for predicting shoreline change[END_REF] showed that approximately two years of monthly observations suffice to calibrate empirical shoreline model parameters on seasonally variable beaches (Torrey Pines in southern California and along the eastern Australian coast). Here, three calibration periods are examined(Figure 8): 1997(Figure 8): -2013 (all data; (all data; 16 years), 2003-2011 (8 years), and 2003-2008 (5 years). The 2003-2008 period lacks an El Niño. Model errors are characterized with the RMSE over 16 years, and with , the difference between the maximum erosion observed and modeled during the 2009-10 El Niño winter. Solana Beach results weakly depended on calibration period (Figure 8, top). At the other sites, longer calibration periods that included an El Niño consistently decreased and RSME over the entire 16-year observation period, which included years of neutral and La Niña conditions (Figure 9a,b). The sparse 1997-2001 data were not well fit, even when 1997-2001 was included in the calibration (not shown). The 2003-2011 calibration period was used. The alongshore variability of the 8-year calibration model coefficients was qualitatively similar to previous work [Y09] based on ~5 years of calibration that did not include El Niño (similar to the 2003-2008 calibration results in this study). Here, the relative magnitudes of the wave energy slope, a 1 , and were reversed compared to Y09 (e.g. Y09 had larger (smaller) magnitude a 1 (

  98 and 2009-10, illustrate seasonal and long-term fluctuations in wave climate and shoreline sand levels at five southern California beaches. An existing, empirical shoreline model driven with hourly wave conditions simulates well the seasonal changes in subaerial beach width (e.g. the cross-shore location of the MSL contour) during non-El Niño years, similar to previous results [Y09]. During El Niño winters the Y09 model overprediction of shoreline erosion is reduced by including the location of erosion resistant

FiguresFigure 1 .

 1 Figures

  Figure 1. (a) Southern California map with wave buoy locations (black squares). (b)

Figure 2 .

 2 Figure 2. (a) Aerial image of Imperial Beach with subaerial substrate and back beach

Figure 3 .

 3 Figure 3. MSL cross-shore position (demeaned and alongshore averaged) versus time

Figure 4 .

 4 Figure 4. MSL cross-shore position (demeaned and alongshore averaged) versus time

Figure 5 .Figure 6 .

 56 Figure 5. (a) Hours of observed H s between 2-3 m, and greater than 3 m (see legend)

Figure 7 .

 7 Figure 7. Modeled MSL position versus time at Torrey Pines (section T8, calibrated

Figure 8 .

 8 Figure 8. Model (linear 4 ) (a) RMSE (all data) and (b) model-data winter (January-

Figure 9 .

 9 Figure 9. Modeled and observed (white circles) MSL versus time at (a) Solana Beach

Figure 10 .

 10 Figure 10. Example model results for Torrey Pines section T8 parameters: (a)

Table 1 .

 1 using their datasets. This publication was prepared by André Doria under NOAA Grant #NA10OAR4170060. California Sea Grant Project R/RCC-01, through NOAA's National Sea Grant College Program and the U.S. Dept. of Commerce. The statements, findings, conclusions and recommendations are those of the author(s) and do not necessarily reflect the views of the aforementioned organizations. Beach Alongshore Distance, Beach Facing Azimuthal Direction, Median Sand Grain Diameter (D 50 ), Beach Slope at MSL, MSL Minimum, Maximum, and Standard Deviation Horizontal Displacement from Average MSL Location, Number of Surveys, and Survey Date Range for Each Site. Collected spring 2006 near the high tide line[Yates et al., 2009b].

	Tables
	Site	Alongshore	Direction	D 50 (mm)	Beach	MSL	Number	Date
		Distance	(deg)		Slope	min/max (	of	Range
		(km)				) (m)	Surveys
	Imperial	4	250-270	0.28 a	0.02-0.05 -23.6/25.8	97	Oct
	Beach					(10.5)		1997-
								Aug 2012
	Torrey	8	260-270	0.15 b	0.01-0.08 -31.5/26.2	226	Oct
	Pines					(9.1)		1997-Jan
								2014
	Solana	2.6	240-265	0.15 b	0.02-0.08 -22.5/20.5	103	Oct
	Beach					(7.6)		1997-
								Aug 2012
	Cardiff	2	260	0.15 b	0.02-0.11 -27.3/22.8	136	Oct
						(9.5)		1997-
								Aug 2012
	Camp	2.5	235	0.18 b	0.02-0.04 -35.8/19.4	72	Oct
	Pendleton					(9.7)		1998-Oct
								2010

a Collected May 2008 in the swash zone

[START_REF] Warrick | Observations of coastal sediment dynamics of the Tijuana Estuary Fine Sediment Fate and Transport Demonstration Project[END_REF]

]. b

Table 3 .

 3 Alongshore averages and standard deviations of optimal model free parameters and R 2 at each site. Average characteristic adjustment timescales a are shown in parenthesis and have units of days. The calibration period is October 2003-October 2011. R 2 is for model runs over all available data.

		a 1				R 2
				3		3
			(mh -1	m	) (mh -1	m	)
	Imperial	-4.5±1.7	-0.92±0.72 -1.06±0.94	3.5±2.5 (21)	3.6±2.4 (64) 0.55±0.13
	Beach
	Torrey	-2.2±0.9	-3.90±2.67 -4.58±2.02	6.6±3.2 (9)	10.3±5.7 (29) 0.58±0.08
	Pines
	Solana	-5.8±2.5	-1.26±0.67 -0.83±0.43	7.6±3.5 (21)	4.0±1.7 (54) 0.60±0.08
	Beach
	Cardiff	-5.4±3.1	-2.16±1.52 -1.49±1.92 14.3±16.9 (14) 4.7±4.1 (57) 0.43±0.15
	Camp	-5.9±2.3	-0.62±0.10 -0.79±0.13	3.5±0.8 (12)	4.8±2.5 (41) 0.38±0.03
	Pendleton
	a	m (1 m) was used for estimating	( ).

Acknowledgements

This research was supported by the United States Army Corps of Engineers and the California Department of Parks and Recreation Division of Boating and Waterways (program manager Reinhard Flick). The paper was substantially improved by reviewer and editor commments. André Doria was supported by Fellowships from the University of California Regents, NDSEG, and the National Science Foundation (GRFP). SANDAG, Coastal Frontiers Corporation, and Fugro are thanked for their assistance in


Cécile Raoult
email: cecile.raoult@edf.fr

Michel Benoit
email: benoit@irphe.univ-mrs.fr

Marissa L Yates
email: marissa.yates-michelin@cerema.fr

Validation of a fully nonlinear and dispersive wave model with laboratory non-breaking experiments

Keywords: nonlinear, dispersive, water waves, potential theory, Zakharov equations

come

Validation of a fully nonlinear and dispersive wave model with laboratory non-breaking experiments

Introduction

Fast and accurate models of wave transformation from deep water to the beach over large spatial scales are needed for coastal and ocean engineering applications, such as the design of shore protection structures or marine renewable energy projects. This objective is not an easy task given the range of spatial and temporal scales to be modelled. To succeed, this type of model requires an accurate representation of nonlinear phenomena playing an important role in wave interactions with structures and bathymetry, an accurate representation of dispersion to propagate waves over a wide range of depths, and the proper simulation of important coastal physical processes like shoaling, refraction, diffraction, breaking, and run-up.

Most fluid flow problems can be described by the Navier-Stokes equations since they account for nonlinearities, vorticity and viscosity. Models based on these equations can be very accurate when studying wave interactions with structures in the surf zone (e.g. [START_REF] Shao | Incompressible SPH simulation of wave breaking and overtopping with turbulence modelling[END_REF]; [START_REF] Lara | RANS modelling applied to random wave interaction with submerged permeable structures[END_REF]). These equations can be solved either with an Eulerian approach as in the widely used code OpenFOAM R � (Higuera et al., 2013a,b) solving the RANS (Reynolds Averaged Navier-Stokes) equations for two incompressible phases, or with a Lagrangian approach computing the trajectories of interacting fluid particles (e.g. [START_REF] Dalrymple | Numerical modeling of water waves with the SPH method[END_REF]). These models are highly accurate when studying local-scale processes, but the domain size and resolution are limited due to the computational time, even with the use of GPU parallelized codes [START_REF] Dalrymple | GPU-Accelerated SPH model for water waves and free surface flows[END_REF]. To model large spatial domains, these codes are usually coupled with more computationally efficient models, such as potential flow models, to simulate the far-field processes [START_REF] Narayanaswamy | SPHysics-FUNWAVE hybrid model for coastal wave propagation[END_REF].

Some assumptions can be made to simplify the problem and thus reduce the computational time. For example, the nonlinear shallow water equations (NLSWE) are derived by depth integrating the Navier-Stokes equations to model waves with a wavelength significantly longer than the water depth (e.g. to model tidal waves), but this set of equations does not take into account wave dispersion. However, by including non-hydrostatic pressure in the NL-SWE and dividing the water depth into a sufficient number of layers, the frequency dispersion of waves can be greatly improved [START_REF] Stelling | An accurate and efficient finite-difference algorithm for non-hydrostatic free-surface flow with application to wave propagation[END_REF][START_REF] Ziljema | Further experiences with computing nonhydrostatic free-surface flows involving water waves[END_REF]Stelling, 2005, 2008;Ziljema et al., 2011). For example, with only two layers, the accuracy of the deep water dispersion relation is similar to that of extended Boussinesq-type models. The dispersion of such non-hydrostatic models can be further improved by optimizing the location of the levels [START_REF] Zhu | Optimization of non-hydrostatic Euler model for water waves[END_REF].

When viscous and turbulent effects are negligible, the flow can be well represented by potential flow theory, which consists of solving the Laplace problem in the fluid domain, supplemented by nonlinear free surface boundary conditions. One way of solving this problem is to use the Boundary Integral Equations Method (BIEM), which projects the problem on the boundary surface of the fluid domain using Green's Function [START_REF] Grilli | An efficient boundary element method for nonlinear water waves[END_REF][START_REF] Wang | An efficient numerical tank for nonlinear water waves, basedon the multi-subdomain approach with BEM[END_REF]. These models enable an accurate description of nonlinear and even overturning waves and are well adapted to simulate wave-structure interactions (e.g [START_REF] Dombre | Simulation of floating structure dynamics in waves by implicit coupling of a fully non-linear potential flow model and rigid body motion approach[END_REF]). This method is mainly used for calculating local-scale interactions owing to the long computational times. However, with the use Fast Fourier Transform (FFT) [START_REF] Newman | Boundary-Element Methods in offshore structure analysis[END_REF][START_REF] Fructus | An explicit method for the nonlinear interaction between water waves and variable and moving bottom topography[END_REF] or Fast Multipole Algorithm methods [START_REF] Fochesato | Numerical modeling of extreme rogue waves generated by directional energy focusing[END_REF], the computational time can be reduced considerably.

Another way of solving the problem is to make additional assumptions about the nonlinear and dispersive properties of waves. By doing a Taylor expansion of the vertical velocity about a specified level and truncating it to a finite number of terms, Boussinesq-type models assume a polynomial variation of the vertical velocity, thus reducing the problem by one dimension. Boussinesq-type models are derived with the assumption that nonlinearity and/or frequency dispersion are weak or moderate [START_REF] Madsen | Higher order Boussinesq-type equations for surface gravity waves: derivation and analysis[END_REF][START_REF] Kirby | Boussinesq models and applications to nearshore wave propagation, surfzone processes and wave-induced currents[END_REF]. Using only a quadratic polynomial approximation of the vertical flow distribution gives poor results for wave propagation in intermediate depths. A lot of work has been done to improve the frequency dispersion following various approaches such as: using higher degree polynomials for the vertical approximation with the Green-Naghdi equations [START_REF] Zhao | On the steady solitary-wave solution of the Green-Naghdi equations of different levels[END_REF], using Padé approximants [START_REF] Agnon | A new approach to highorder Boussinesq models[END_REF] combined with an expansion of the Laplace solution from an arbitrary level [START_REF] Madsen | A new Boussinesq method for fully nonlinear waves from shallow to deep water[END_REF], and resolving in two arbitrary layers to maintain low-order spatial derivatives [START_REF] Lynett | A two-layer approach to wave modelling[END_REF][START_REF] Chazel | A double-layer Boussinesqtype model for highly nonlinear and dispersive waves[END_REF]. Additional modeling approaches include those of [START_REF] Kennedy | Boussinesqtype equations with improved nonlinear performance[END_REF], [START_REF] Fuhrman | Numerical solution of fully non-linear and highly dispersive Boussinesq equations in two horizontal dimensions[END_REF][START_REF] Fuhrman | Numerical solution of fully non-linear and highly dispersive Boussinesq equations in two horizontal dimensions[END_REF][START_REF] Engsig-Karup | Nodal DG-FEM solution of high-order Boussinesq-type equations[END_REF].

In the present work, a fully nonlinear potential flow theory model is developed solving the Zakharov equations [START_REF] Zakharov | Stability of periodic waves of finite amplitude on the surface of a deep fluid[END_REF]. The temporal evolution of the free surface elevation η and the free surface velocity potential Φ are given as a function of these two variables and the vertical velocity at the free surface w. The primary challenge is to express the vertical velocity w as a function of η and Φ, a problem commonly called 'Dirichlet-to-Neumann' or DtN. One possibility is to solve directly the Laplace equation using finite element [START_REF] Wu | Numerical simulation of sloshing waves in a 3D tank based on a finite element method[END_REF][START_REF] Ma | Finite element simulation of fully non-linear interaction between vertical cylinders and steep waves. Part 1: Methodology and numerical procedure[END_REF] or finite difference [START_REF] Li | A three dimensional multigrid model for fully nonlinear water waves[END_REF][START_REF] Engsig-Karup | An efficient flexible-order model for 3D nonlinear water waves[END_REF] methods. When using finite difference methods, [START_REF] Kreiss | Comparison of accurate methods for the integration of hyperbolic equations[END_REF] and [START_REF] Bingham | On the accuracy of finite-difference solutions for nonlinear water waves[END_REF] recommend using fourth-order schemes with a stretched vertical grid (clustering points near the free surface) instead of using second-order schemes with a regular grid. When considering rectangular domains with a flat bottom, a high-order spectral approach (HOS) is optimal [START_REF] Dommermuth | A high-order spectral method for the study of nonlinear gravity waves[END_REF][START_REF] West | A new numerical method for surface hydrodynamics[END_REF][START_REF] Chern | A pseudospectral σ-transformation model of 2-D nonlinear waves[END_REF][START_REF] Ducrozet | A comparative study of two fast nonlinear free-surface water wave models[END_REF]. This method is faster than finite difference methods but less flexible with regard to the specified geometry and bathymetry, even if progress has been made in taking into account variable and moving bottoms [START_REF] Smith | An operator expansion formalism for nonlinear surface waves over variable depth[END_REF][START_REF] Guyenne | A high-order spectral method for nonlinear water waves over moving bottom topography[END_REF][START_REF] Gouïn | Validation of a nonlinear spectral model for water waves over a variable bathymetry[END_REF]. An additional approach is to use a spectral method only in the vertical dimension either by expanding the velocity potential with a local mode series [START_REF] Belibassakis | A coupled-mode system with application to nonlinear water waves propagating in finite water depth and in variable bathymetry regions[END_REF] or by projecting it on a polynomial basis [START_REF] Kennedy | A fully-nonlinear computational method for wave propagation over topography[END_REF][START_REF] Tian | A numerical model on the interaction between nearshore nonlinear waves and strong currents[END_REF]. By using high-order finite difference schemes in the horizontal, these models maintain a flexible approach for variable domain geometries and bathymetry. A comparison between a vertical spectral approach and a finite difference approach in the vertical dimension shows the improved accuracy and efficiency of the spectral method in 1DH [START_REF] Yates | Accuracy and efficiency of two numerical methods of solving the potential flow problem for highly nonlinear and dispersive water waves[END_REF] and 2DH [START_REF] Christiansen | Hybridspectral model for fully nonlinear numerical wave tank[END_REF]. Others approaches can be used, such as the extension of the DtN operator as a sum of global convolution terms and local integrals with kernels decaying quickly in space [START_REF] Clamond | A fast method for fully nonlinear water-wave computations[END_REF][START_REF] Fructus | An efficient model for three-dimensional surface wave simulations: Part I: Free space problems[END_REF]. Here, the combination of high-order finite difference schemes for the horizontal dimension and a spectral approach using Chebyshev polynomials in the vertical dimension is tested with emphasis on the ability of the model to accurately represent nonlinear and dispersive phenomena.

The mathematical model and its numerical implementation are presented in Section 2. In Section 3, the ability of the model to simulate deep water cases is demonstrated before the model is validated by comparing simulation results to experimental data for four non-breaking wave test cases: (1) nonlinear dynamics of free and bound components in a wave train generated by a piston-type wavemaker in constant water depth, (2) shoaling of a regular wave train on beach with constant slope up to the breaking point, (3) propagation of regular waves over a submerged bar [START_REF] Dingemans | Comparison of computations with Boussinesq-like models and laboratory measurements[END_REF], and (4) propagation of nonlinear irregular waves over a barred beach [START_REF] Becq-Girard | Non-linear propagation of unidirectional wave fields over varying topography[END_REF]. Conclusions and ongoing developments are summarized in the final Section.

Mathematical model and numerical implementation

Mathematical model

By assuming an inviscid and homogeneous fluid of constant density, and irrotational flow, potential flow theory can be used. The velocity potential Φ(x, z, t), where x = (x, y), must satisfy the Laplace equation in the fluid domain:

∇ 2 Φ + Φ zz = 0. (1)
Here, ∇f ≡ (f x , f y ) is the horizontal gradient operator, and partial derivatives are denoted with subscripts. At the free surface z = η(x, t), which is assumed to be single-valued in x (no overturning waves), and the (nonlinear) kinematic and dynamic boundary conditions (BC) for z = η(x, t) are:

η t = -∇η. ∇Φ + Φ z , (2)
Φ t = -gη - 1 2 � (∇Φ) 2 + (Φ z ) 2 � . (3)
In (3), the atmospheric pressure at the free surface is assumed to be uniform (and chosen to be 0 by convention), surface tension is neglected, and g is the acceleration of gravity. At the bottom z = -h(x), assumed constant in time, the kinematic (impermeability) BC is:

∇Φ • ∇h + Φ z = 0. ( 4
)
At the lateral boundaries, periodic, Dirichlet or Neumann BCs can be imposed.

By defining the velocity potential at the free surface Φ(x, t) ≡ Φ(x, z = η(x, t), t), the free surface conditions (2) and (3) are reformulated respectively as [START_REF] Zakharov | Stability of periodic waves of finite amplitude on the surface of a deep fluid[END_REF]:

η t = -∇η. ∇ Φ + w � 1 + (∇η) 2 � , ( 5
) Φt = -gη - 1 2 (∇ Φ) 2 + 1 2 w2 � 1 + (∇η) 2 � , (6)
where w(x, t) ≡ Φ z | z=η(x,t) is the vertical velocity at the free surface. Equations (5-6) involve only free surface quantities, which are functions of x and t only. To integrate these equations in time, it is necessary to determine w(x, t) from (η(x, t), Φ(x, t)), the DtN problem.

As described in [START_REF] Yates | Accuracy and efficiency of two numerical methods of solving the potential flow problem for highly nonlinear and dispersive water waves[END_REF], the DtN problem is resolved by solving a Laplace boundary value problem (BVP) for the potential Φ, namely the Laplace equation (1), supplemented with the Dirichlet BC at the free surface:

Φ = Φ at z = η(x, t), (7)
the bottom BC (4), and specified lateral BCs. Following the work of [START_REF] Tian | A numerical model on the interaction between nearshore nonlinear waves and strong currents[END_REF] a spectral method is applied in the vertical dimension. An outline of the method, detailed in [START_REF] Yates | Accuracy and efficiency of two numerical methods of solving the potential flow problem for highly nonlinear and dispersive water waves[END_REF], is given hereafter for the case of a single horizontal dimension (i.e. x = x), and the extension to two horizontal dimensions is straightforward. First, the fluid domain is mapped onto a strip of constant height by introducing the vertical coordinate s, which varies from s = -1 at the bottom to s = 1 at the free surface:

s(x, z, t) = 2z + h -(x, t) h + (x, t) , (8)
where h + (x, t) = h(x)+η(x, t) and h -(x, t) = h(x)-η(x, t). The BVP is then reformulated in this transformed space (x, s) for Φ(x, z, t) ≡ ϕ(x, s(x, z, t), t).

The Laplace equation (1) in the fluid domain, the bottom BC (4) and the free surface Dirichlet BC (7) become:

ϕ xx + 2s x ϕ xs + (s 2 x + s 2 z )ϕ ss + s xx ϕ s = 0 in the fluid domain, (9) h + h x ϕ x + 2(1 + h 2 x )ϕ s = 0 for s = -1, ( 10
) ϕ(x, 1) = Φ(x) for s = 1. ( 11
)
Then a spectral approach is used in the vertical: the vertical variation of the velocity potential is approximated by a linear combination of Chebyshev polynomials of the first kind, T n (x), where n = 0, 1, 2, ... indicates the order of the polynomial:

ϕ(x, s) ≈ N T � n=0 a n (x)T n (s), ( 12
)
where N T is the maximum order of the Chebyshev polynomials. These polynomials are easy to evaluate, converge rapidly, and have a large convergence domain. [START_REF] Tian | A numerical model on the interaction between nearshore nonlinear waves and strong currents[END_REF] and [START_REF] Yates | Accuracy and efficiency of two numerical methods of solving the potential flow problem for highly nonlinear and dispersive water waves[END_REF] have shown that values of N T lower than 10 (typically in the range [5,8]) are sufficient to provide high accuracy for a variety of test cases. Furthermore, with the inner product defined as:

< f, g >≡ � 1 -1 f (s)g(s) √ 1 -s 2 ds, (13)
the T n polynomials form an orthogonal basis over the range [-1, 1], since it can be shown that:

< T n , T p >=    0 if n � = p, π if n = p = 0, π 2 if n = p � = 0. (14)
For an arbitrary function f (s) in [-1, 1], the following operator is defined:

< f > p ≡ 2 πC p < f, T p >, with � C 0 = 2, C p = 1 for p > 0. ( 15
)
The orthogonality relations ( 14) then become: < T n > p = δ np , where δ np is the Kronecker delta. By substituting the approximation (12) in the BVP (9-11), the Laplace equation ( 9) becomes (dropping the x argument of a n (x) coefficients and the s argument of T n (s) polynomials for brevity):

N T � n=0 [a �� n T n + 2s x a � n T � n + a n ((s 2 x + s 2 z )T �� n + s xx T � n )] = 0, (16)
with:

s x = (h - x -sh + x )/h + , s z = 2/h + , s zz = 0, (17) s xx = [(h - xx h + -2h - x h + x ) + s(2(h + x ) 2 -h + xx h + )]/(h + ) 2 , (18)
and

a � n ≡ da n dx , a �� n ≡ d 2 a n dx 2 , T � n ≡ dT n ds and T �� n ≡ d 2 T n ds 2 . ( 19
)
Then the Chebyshev-tau method, a variant of the Galerkin method (see e.g. [START_REF] Boyd | Chebyshev and Fourier Spectral Methods: Second Edition, Revised[END_REF]), is applied in the vertical to ( 16) to remove the s dependence. The operator < . > p is applied to ( 16) for p = 0, 1, ..., N T -2, yielding N T -1 equations:

a �� p + N T � n=0 C pn a � n + N T � n=0 D pn a n = 0, p = 0, 1, ...., N T -2, ( 20
)
with

C pn ≡ (m 011 B p01n +m 111 B p11n )/m 020 and D pn ≡ (m 002 B p02n +m 102 B p12n + m 202 B p22n + m 001 B p01n + m 101 B p11n )/m 020 . The terms m ijk depend only on h + (x, t), h -(x, t)
, their spatial derivatives (see equations ( 21)-( 28) in Tian and Sato ( 2008)), and the coefficients:

B pikn ≡< s i ∂ k T n ∂s k (s) > p , (21)
which are constant and can be calculated analytically once at the beginning of each simulation.

The N T -1 equations ( 20) are supplemented by the two BCs at the bottom (10) and the free surface (11), which become respectively:

N T � n=0 [(-1) n h + h x a � n + (-1) n-1 2n 2 (1 + h 2 x )a n ] = 0, ( 22
) N T � n=0 a n = Φ. (23)
Thus, at each node x, there are N T + 1 equations to solve for the N T + 1 unknown coefficients a n (x) for n = 0, 1, ..., N T . A similar treatment is applied for lateral boundary points (see [START_REF] Yates | Accuracy and efficiency of two numerical methods of solving the potential flow problem for highly nonlinear and dispersive water waves[END_REF]).

Once the a n (x) are determined, the vertical velocity at the free surface w is readily obtained:

w(x) = φ s s z � � s=1 = 2 h + (x) N T � n=0 a n (x)n 2 , ( 24
)
and equations (5-6) can be stepped forward in time.

Numerical implementation

To simulate this mathematical model, a numerical model, named MIS-THYC, has been developed and coded in Fortran in one horizontal dimension (x, z). The horizontal domain is discretized with N X nodes in the x direction, with constant or variable Δx. High-order finite difference approximations (fourth order here) are used to calculate first and second-order derivatives in x, following [START_REF] Fornberg | Generation of finite difference formulas on arbitrarily spaced grids[END_REF]. The explicit four-step fourth-order Runge-Kutta (RK4) scheme with a constant time step is used to integrate ( 5) and ( 6) in time.

At each step of the RK4 scheme, the discretization of the Laplace BVP in (x, s), with N X nodes in x and N T the maximum order of the Chebyshev polynomials, results in a system of N X (N T + 1) linear equations for the coefficients a n (x i ) for n = 0, ...., N T and i = 1, ...., N X . The corresponding matrix is sparse, and the system is currently solved in MISTHYC using the direct solver MUMPS ("MUltifrontal Massively Parallel Solver", v4.10.0) [START_REF] Amestoy | A fully synchronous multifrontal solver using distributed dynamic scheduling[END_REF][START_REF] Amestoy | Hybrid scheduling for the parallel solution of linear systems[END_REF] using the default settings. Iterative solvers could be also used and will be tested in the future.

Validation test cases

The model is validated by comparing the simulation results to stream function solutions of regular nonlinear waves propagating in deep water, and then to measurements from four laboratory experiments of non-breaking waves: (1) nonlinear dynamics of free and bound components in a wave train generated by a piston-type wavemaker in constant water depth, (2) shoaling of a regular wave train on beach with constant slope up to the breaking point, (3) propagation of regular waves over a submerged bar [START_REF] Dingemans | Comparison of computations with Boussinesq-like models and laboratory measurements[END_REF], and (4) propagation of nonlinear irregular waves over a barred beach [START_REF] Becq-Girard | Non-linear propagation of unidirectional wave fields over varying topography[END_REF].

Applicability of the model to large relative water depths

Incident wave conditions are characterized by two non-dimensional numbers: (i) the relative water depth µ ≡ kh (equivalently h/L) and (ii) the wave steepness � ≡ ka = kH/2 (equivalently H/L), where h is the local water depth, k is the wave number of the fundamental wave for regular waves (or the peak frequency for irregular waves), L = 2π/k is the wavelength, H is the wave height (crest-to-trough height for regular waves or significant wave height H m0 for irregular waves), and a = H/2 is the wave amplitude. Given h and T , k and L are computed from the linear dispersion relation.

In the following four laboratory experiments, the relative water depth kh is moderate or small (in the range [0.53, 0.73]) where incident waves are imposed. As waves propagate, higher harmonics are generated with frequencies that are multiples of the fundamental frequency, with consequently larger values of kh. To represent properly the dynamics of the higher harmonics, the model must simulate accurately wave dispersion for values of kh much larger than that of the fundamental frequency. [START_REF] Yates | Accuracy and efficiency of two numerical methods of solving the potential flow problem for highly nonlinear and dispersive water waves[END_REF] showed that deep water µ > π (or h/L > 1/2) test cases can be modeled accurately by increasing the maximum order N T of the Chebyshev polynomials (12) for: (1) the propagation of a stable, regular nonlinear wave in constant depth, and (2) the periodic motion of a nonlinear standing wave in constant depth. In the first case, the initial wave profile, obtained from the stream function method (order 20), corresponds to a relative water depth µ = 2π ≈ 6.28 (or h/L = 1) and steepness � = π/10 ≈ 0.314 (or H/L = 10%). After propagating for 25 wave periods, errors in the free surface position remained small (< 0.1%) for N T > 7 (see e.g. section 3.1.3 and Figures 7 and 8 in [START_REF] Yates | Accuracy and efficiency of two numerical methods of solving the potential flow problem for highly nonlinear and dispersive water waves[END_REF]). In the second case, after 100 periods of motion of a standing wave with µ = 3 (or h/L ≈ 0.5) and � = 0.42, and with N T = 7 and sufficient horizontal resolution, the free surface position errors also remained relatively small (< 3%). Considering the increase in computational time with N T , [START_REF] Yates | Accuracy and efficiency of two numerical methods of solving the potential flow problem for highly nonlinear and dispersive water waves[END_REF] concluded that the optimal range of N T is [7; 15] in such deep water conditions.

To further demonstrate the ability of the model to simulate accurately waves in deep water, new results are presented here for 25 periods of propagation of a stable, regular nonlinear wave with wave steepness � = π/10 ≈ 0.314 (or H/L = 10%, same as in [START_REF] Yates | Accuracy and efficiency of two numerical methods of solving the potential flow problem for highly nonlinear and dispersive water waves[END_REF]) and µ = 4π ≈ 12.56 (or h/L = 2), which is 4 times the deep water threshold (and twice that shown in [START_REF] Yates | Accuracy and efficiency of two numerical methods of solving the potential flow problem for highly nonlinear and dispersive water waves[END_REF]). The computational domain is one wavelength (L = 64 m), and based on the sensitivity tests shown in [START_REF] Yates | Accuracy and efficiency of two numerical methods of solving the potential flow problem for highly nonlinear and dispersive water waves[END_REF], the spatial and temporal discretization are Δx = L/96 and Δt = T /75, respectively, corresponding to CFL = 1.28. The initial wave is computed from the stream function solution (order 20), and after each period the computed wave should be unchanged from the initial wave profile. Final wave profiles for N T = 5, 7, 8, 10, 12, 15, and 20 are presented in Figure 1a and b for µ = 2π and 4π, respectively. Both cases converge visually when N T is increased, but for the deeper water case, the errors decrease more slowly and the simulation with N T = 5 does not remain stable for long time simulations. A quantitative comparison of the phase difference with respect to its theoretical position after 25 periods of wave propagation (Figure 2) shows that the accuracy of the simulations generally increases to an approximately constant value with increasing N T . For N T < 12, the phase errors are larger for the deeper water case (kh = 4π), but for N T ≥ 12, the phase differences decrease and the model simulates accurately long term wave propagation for both cases. These tests show that the model is able to simulate accurately deep water waves by increasing the order of the maximum Chebyshev polynomial (e.g. using N T = 10 for kh = 2π and N T = 12 for kh = 4π). A harmonic analysis of free surface elevation time series (after the steady state is reached) decomposes the signal into a discrete sum:

η(t) = a 0 + N � n=1 a n cos(nωt + ϕ n ), ( 25
)
where t is time, ω = 2π T is the angular frequency of the wavemaker, and a n and ϕ n are the amplitude and phase of the harmonic component n. The phase difference between the first and the second harmonic is defined as Δ ϕ 1,2 ≡ ϕ 2 -2ϕ 1 , following [START_REF] Chapalain | Observed and modeled resonantly interacting progressive water-waves[END_REF]. The simulated spatial evolution of the amplitudes of the first four harmonics agrees well with the experiments (Figure 3a). Overall, the model correctly represents the energy transfers between the different harmonic components, as well as the resultant beat lengths. However, a decrease in the second harmonic amplitude (after x = 19 m) is observed in the experimental data but is not reproduced in the simulations. This could be explained by dissipation in the experiments that is not taken into account in the model. This phenomenon is more noticeable for short waves, hence more visible for the higher harmonics. The spatial evolution of the phase difference between the first and second harmonic is also reproduced well. The phase difference oscillates between -π/2 and +π/2 with the same periodicity as the harmonic amplitudes. Zero phase difference occurs when either the first harmonic is maximum and the second harmonic is minimum, or the contrary.

The variation of the free surface elevation thus depends on the location in the wave channel, as shown in Figure 4 at x = 4 m, 7 m, 10 m and 14 m. The simulated free surface position qualitatively agrees well with the measurements (compared to Figure 3 of [START_REF] Chapalain | Observed and modeled resonantly interacting progressive water-waves[END_REF]). When the first and second harmonics are in phase, the free surface profile is either quasi-sinusoidal when the first harmonic is maximal and the second minimal (e.g. x = 14 m), or cnoidal when the second harmonic is maximal and the first minimal (e.g. x = 7 m). However, when the first and second harmonics are out of phase, the waves are vertically asymmetrical with either a steeper (gentler) wave front and a gentler (steeper) rear slope if the phase difference is positive (negative) (e.g x = 10 m and x = 4 m).

Nonlinear wave shoaling on a gently sloping beach

The second test case simulates the shoaling of regular nonlinear waves on a gently sloping (1/35) beach following the experiments of Ting and Kirby (1994) (TK94). In the spilling breakers experiment, waves were generated in h = 0.4 m water depth with a wave period of T = 2.0 s and wave height of H = 0.125 m, corresponding to long waves with µ = 0.68 (h/L = 0.104) and � = 0.1 (or H/L = 3.3 %). In the numerical model, waves are generated with the same wave characteristics as the experiments by searching for the analytical solution to the Zakharov equations and numerically (iteratively) finding the solution. Dirichlet boundary conditions are applied to generate waves at the left boundary and in a relaxation zone that extends 2 -3L into the domain (-15.0 ≤ x ≤ -4.5 m). In the experiments, waves begin to break around x = 6.0 m. The present numerical model is not able to take into account the effects of wave breaking; therefore, waves are absorbed in an absorption zone of constant depth (h = 0.2286 m) from x = 6.0 m to x = 15.0 m (shaded gray area, Figure 5). The model results will only be compared to the experimental data before the absorption zone. For these nonlinear waves, the model domain was discretized with Δx = 0.03 m (≈ L/100, with L = 3.8434 m, the wavelength calculated by the numerical solution for h = 0.4 m), the time step was Δt = 0.01 s (T /200), and N T = 7. As the waves become highly nonlinear approaching the breaking limit, high frequency instabilities develop in the numerical model. For the selected spatial and temporal discretization, a fourth-order Savitzky-Golay filter was applied over 9 points (twice the stencil size) at the end of each time step to smooth these high frequency instabilities.

The measured crest and trough elevation envelope is reproduced well by the model during the shoaling phase up to the breaking point (x = 6.0 m, Figure 5a). After this point, the absorption zone dissipates wave energy but does not try to reproduce the effects of wave breaking, and the simulated and measured free surface elevation decay at different rates. Simulated free surface elevation time series at four positions in the wave tank (Figure 6) agree well with the wave gauge measurements. At the bottom of the sloping beach profile, waves are nearly symmetric (Figure 6a). As they progress up the slope and shoal, the wave profiles become asymmetric with steeper wave fronts and gentler rear faces, and more skewed with more pointed peaks and wider troughs (Figure 6b-d). Even near the breaking point, the simulated wave profiles agree well with the measurements, with only small differences in the wave troughs (Figure 6c-d).

To quantify the wave nonlinearity, higher order statistical moments were calculated from the free surface elevation time series, including the skewness (or horizontal asymmetry) and kurtosis (or measure of the flatness of the free surface elevation distribution). The skewness is defined as the normalized, centered, third-order moment of the free surface elevation:

S = �(η -�η�) 3 � σ 3 , ( 26
)
and the kurtosis as the normalized, centered, fourth-order moment of the free surface elevation:

K = �(η -�η�) 4 � σ 4 , ( 27
)
where the variance σ 2 of the free surface elevation is:

σ 2 = �(η -�η�) 2 �, (28)
with �-� denoting the time-average operator.

As observed in the free surface time series, the waves progressively develop narrow, peaked crests and wide, flat troughs as they shoal on the slope, causing an increase in skewness (Figure 7a). With these changes in the wave profile, the free surface elevation distribution becomes less Gaussian, thus causing an increase in the kurtosis (Figure 7b). The simulations reproduce well the measured increases in skewness and kurtosis as the waves shoal, with small differences appearing as the waves become highly nonlinear approaching the breaking point. After the breaking point, the model is not able to simulate accurately the wave profile and the dissipation of energy, and a wave breaking dissipation model is currently being developed to take into account these effects.

Regular waves over a submerged bar

Beji and Battjes (1993) and then Dingemans (1994) (D94) performed a series of experiments of regular wave propagation over a submerged trapezoidal bar. These experiments are now a standard test case for wave models since both nonlinear and dispersive effects are important when waves propagate over the top of the bar. The bottom profile is shown in Figure 8 at the scale of the experiments of [START_REF] Beji | Experimental investigation of wave propagation over a bar[END_REF]. The water depth is h = 0.40 m offshore and reduces to a minimum of 0.10 m on top of the bar. The front slope of the bar is 1:20, and the rear slope is 1:10. Eleven wave probes recorded the free surface elevation time series in the experiments (see positions in Figure 8).

Here, case A is studied with an incident wave height H = 2.0 cm and a period T = 2.02 s, corresponding to long waves with µ = 0.67 (or h/L = 0.107) offshore of the bar and relatively small steepness � = 0.017 (or H/L = 0.53%). Under these conditions, the incident wave train is significantly affected during the propagation over the submerged bar, with transfers of energy toward higher harmonics.

The model domain has a regularly spaced grid with Δx = 0.05 m (i.e. about L/75 offshore of the bar) and extends from x = -6 m to 30 m (721 nodes along the x axis). Waves are generated in an 8-m wide relaxation zone at the left boundary of the domain (-6 ≤ x ≤ 2 m) using a secondorder Stokes solution to impose the free surface elevation and the velocity potential. Waves are absorbed in a 5-m wide relaxation zone (25 ≤ x ≤ 30 m) applied in front of the fully reflective right boundary to avoid reflections.

Waves are propagated during 25 wave periods (i.e. 50.5 s) with a time step Δt = T /100 = 0.0202 s.

Times series of free surface profiles computed with N T = 7 are presented at probes 4 to 11 in Figure 9. The choice of this value of N T will be discussed at the end of this section. When waves propagate over the front slope of the bar, the wave height and steepness increase due to shoaling effects (probes 4 to 6, Figure 9(a-c)). The wave profile becomes asymmetric due to nonlinear wave-bottom interactions that create higher frequency bound harmonic components. These harmonics are released in the shallowest region and on the rear slope of the bar and then propagate at their own phase speed (probes 7 and 8, Figure 9(d-e)). After the bar, the measured wave profiles vary significantly between the probes due to the differences in celerity of the free wave components. At the last three probes (probes 9 to 11, Figure 9(f-h)), the model reproduces well the complex wave profiles, including the dispersive (high frequency) components.

In order to examine more closely the energy transfers between harmonics, a Fourier analysis of the measured and computed wave signals was completed. The model accurately reproduces the spatial evolution of the amplitudes of the first six harmonics, corresponding to frequencies f = 1/T (fundamental component) to 6f (Figure 10). The amplitude of the fundamental harmonic increases due to shoaling on the submerged bar up to the bar's crest and then decreases (starting from about x = 12 m) due to transfers of energy toward higher harmonics. The observed oscillations are hypothesized to be caused by reflections in the wave channel. The second harmonic (2f ) increases as the waves shoal on the front slope of the bar and continues increasing until x = 16 m, after which its amplitude fluctuates. The higher harmonics start increasing in amplitude at shallower depths, and mainly develop around the bar crest (12 ≤ x ≤ 15 m). Harmonics 4f to 6f decrease after the bar, which is not the case for harmonics 2f and 3f . In addition, after the bar, the second harmonic has the largest amplitude, and the amplitude of the third harmonic is nearly comparable to (though smaller than) that of the fundamental component. The model results agree well with the data up to the sixth harmonic, with the exception of slight differences at some locations in the amplitude of the second harmonic (e.g. at x = 19 m).

To test the sensitivity of the model to the parameter N T (all other parameters were kept constant), simulations were run with N T = 3, 4, 5, 7 and 10 and compared to the observations (Figure 11). In the offshore part of the flume and up to the submerged bar (i.e. up to probe 8), only the simulation results with N T = 3 differ significantly from those with higher N T , which are superimposed and agree well with the measured time series. At the last three probes, where dispersive effects are more important, the results with N T = 3 clearly deteriorate, and results with N T = 4 and 5 also show increasing differences with the measurements. Results with N T = 7 and 10 remain superimposed and in good agreement with the measurements at all probes. Therefore, N T = 7 was chosen to optimize model's accuracy and efficiency.

Random waves over a barred beach

The last test case simulates the propagation of irregular nonlinear waves over a barred beach, reproducing the wave flume experiments of Becq-Girard et al. (1999) (B99). The bathymetric profile of these experiments (Figure 12) was specifically designed to study nonlinear wave effects in shallow water (triad interactions). Irregular waves were generated with a piston-type random wavemaker using a JONSWAP wave spectrum with a peak-enhancement factor of γ = 3.3. The bottom profile was created with smooth metal sheets to minimize bottom friction dissipation, and a beach absorber was included on the upper part of the beach to reduce wave reflection. Resistive-type wave probes measured the free surface elevation at 16 locations in the wave flume (Figure 12) during the 40-minute experiment with a sampling time step Δt = 0.07 s.

The model computational grid extends from x = -5 m to x = 25 m (with the foot of the bar at x = 0 m). Waves are generated in a 5-m relaxation zone by imposing the velocity potential at the left boundary and correcting the free surface position and velocity potential in the relaxation zone. Nonbreaking irregular waves are simulated with significant wave height of H m0 = 3.4 cm and a peak period of T p = 2.39 s in the deepest part of the domain (h = 0.65 m). Waves are absorbed in a 10-m long relaxation zone at the right boundary. Time series of these variables are reconstructed using linear wave theory to sum the components of the wave spectrum obtained from the free surface measurements at probe 2 (located at the foot of the submerged bar). The computational grid is regularly meshed with Δx = 0.05 m, and N T = 7, as in the previous test case. The total simulation time is 2380 s (approximately 39.7 min), with a time step equaling the sampling time step of the free surface elevation probes, Δt = 0.07 s.

The measured and simulated wave variance spectra agree well (Figure 13, spectra at probes 2, 5, 7, 9, 11, 13, 15, and 16). The main spectral peak increases from probes 2 to 5 due to wave shoaling. In addition, energy is transferred from lower to higher frequencies, particularly from the peak frequency to its super-harmonics. This phenomenon is visible at probe 5 with the appearance of the second harmonic peak (2f p ). When the water depth becomes nearly constant (probes 7, 9, and 11), the second and higher harmonic peaks become more pronounced. A peak at the fifth harmonic (5f p ) becomes visible in the spectra at probes 9 and 11, and its amplitude is reproduced well by the model. On the back side of the bar, the energy transfer reverses back to the lower harmonics (in particular to the second harmonic). At probe 13, the peak of the fifth harmonic disappears, and at probe 15, the peaks of the third and fourth harmonics also nearly disappear. Only the second harmonic and main spectral peaks remain visible in the trough. Finally, at probe 16, the third harmonic (3f p ) peak reappears due to the new decrease in the water depth.

The simulated spatial evolution of the first five harmonic amplitudes agrees well with the experimental observations (Figure 14). In the deep end of the wave tank, the amplitude of the first harmonic is dominant. Between 0 m and 5 m, the first harmonic amplitude increases due to shoaling while the higher harmonics remain stable. After 5 m, the first harmonic amplitude decreases while higher harmonic amplitudes increase due to transfers of energy from lower to higher frequencies. Around 9 m, the water depth increases again, the first harmonic amplitude continues decreasing and the third, fourth, and fifth harmonic amplitudes also begin to decrease, with an energy transfer to the second harmonic. Finally, in the shallowest part of the domain, the energy transfer from the second to the third harmonics begins again. The oscillations visible in the first and second harmonic amplitudes are likely caused by reflections from the lateral boundaries.

To further evaluate the spatial spectral evolution, a set of integral wave parameters are calculated. From the variance density spectrum (E(f )), the significant wave height H m0 = 4 √ m 0 and mean wave period

T e = T m-1,0 = m -1 m 0 or T m0,2 = ( m 0 m 2 )
1 2 can be expressed in terms of the n-th moment (m n ) of the variance density spectrum where m n = � ∞ 0 f n E(f ) df . The 0 th moment corresponds to the free surface elevation variance or σ 2 (28).

Nonlinear effects are also visible in the spatial evolution of these parameters, which are globally estimated well by the model (Figure 15). The significant wave height evolves similarly to the first harmonic amplitude shown in Figure 14. It increases as the waves shoal, decreases in the trough and finally increases again as the water depth decreases approaching the beach.

The simulated H m0 agrees well with the measured values, with only a slight overestimation for x > 7 m and a maximum difference of 8.8%. The evolution of the mean wave period is similar for the two definitions considered (T m-1,0 and T m0,2 ). The mean period initially decreases when a reduction in the energy in the low frequency range of the spectrum is compensated by an increase in the high frequency range. The subsequent release of higher harmonics in the trough leads to an increase in the mean wave period that persists along the tank. The largest differences in mean period occur near the end of the tank, with errors of less than 3.5% and 7% for T m-1,0 and T m0,2 , respectively.

To further analyze the simulated wave nonlinearity, the skewness (S) or horizontal asymmetry coefficient, vertical asymmetry coefficient (A) and kurtosis (K) are compared to those of the measurements. The skewness can be defined from the free surface elevation (26) or as:

S = � +∞ m=-∞ � +∞ n=-∞ Re[B m,n ] m 3/2 0 , (29)
and the vertical asymmetry coefficient as:

A = � +∞ m=-∞ � +∞ n=-∞ Im[B m,n ] m 3/2 0 , (30)
where B m,n is the complex bispectrum. Finally, the kurtosis measures the flattening of the free surface distribution, which is related to the probability of occurrence of high waves (27). For a linear sea state, both the vertical and horizontal asymmetries are zero. Here (Figure 16), the simulated skewness and vertical asymmetry are approximately zero at the left side of the domain and evolve along the bathymetric profile in close agreement with the measurements. The spatial evolution of the kurtosis also begins with a value of approximately 3, typical of a linear (Gaussian) sea state, and then increases in shallower water, reaching a maximum in the shallowest zone. The model reproduces well the spatial evolution of the kurtosis, only slightly underestimating the maximum.

This test case validates the ability of the model to simulate the generation, propagation, and absorption of irregular, non-breaking waves, including nonlinear processes such as wave shoaling and wave interactions causing the transfer of energy between higher and lower harmonics.

Conclusions

Comparisons of simulation results to four laboratory experiments with non-breaking waves demonstrate that the model is able to simulate accurately nonlinear wave propagation, including nonlinear and dispersive effects. In the first set of experiments studying nonlinear resonant interactions in constant depth, the transfer of energy from the principal wave to its second harmonic was reproduced well, including the beat length of the resonant interaction. In the second case modeling the shoaling of highly nonlinear waves on a gentle beach, the model accurately simulated the changes in the wave profile, including the skewness and kurtosis, up to the breaking point. The last two sets of experiments investigated the propagation of regular and irregular nonlinear waves over a submerged bar and a barred beach, respectively, including the generation and propagation of higher-order harmonics. In both cases, the simulations reproduced well the observations, including the amplitudes of the first five harmonics.

Using the proposed spectral approach to solve for the velocity potential, the model converges exponentially with a flexible level of accuracy depending of the choice of N T , the highest order Chebyshev polynomial. A sensitivity test verified previous results showing that the model simulates accurately wave propagation for N T = 7 for real tests cases over variable bathymetry. To simulate well deep-water conditions, N T can be increased to 10-12, as shown by the propagation of regular nonlinear waves over a flat bottom with relative water depths as high as kh = 4π ≈ 12.6.

Ongoing model development includes taking into account the effects of wave breaking and validating the accuracy of simulations with moving bottoms (e.g. tsunami generation). The model is also currently being extended to 2DH using unstructured meshes to reach the objective of modeling the coastal zone. [START_REF] Ting | Observation of undertow and turbulence in a laboratory surf zone[END_REF] wave skewness and kurtosis at wave gauges leading up to and in the beginning of the wave absorption zone (light gray shaded area starting at x = 6.0 m).

Acknowledgments

  3.2. Nonlinear dynamics of a wave train generated by a piston-type wavemaker in constant water depth The first case consists of simulating the propagation of waves generated by the sinusoidal movement of a piston-type wavemaker over a flat bottom, based on the flume experiments of Chapalain et al. (1992) (C92). The results presented here correspond to trial A (piston stroke amplitude e = 7.8 cm and period T = 2.5 s), with a constant water depth h = 0.4 m. The wavelength of the fundamental component is L = 4.74 m from the linear dispersion relation, corresponding to long waves with µ = 0.53. The model domain is regularly meshed with Δx = 0.1 m (≈ L/47) and extends far enough to prevent reflection from the right boundary. The waves were propagated during 16 periods (i.e. 40 s) with a time step Δt = T /40 = 0.0625 s, with maximum order of the Chebyshev polynomial N T = 7. The model is forced by imposing at the left boundary a sinusoidal time varying horizontal velocity that is uniform in the vertical.

Figure 1 :

 1 Figure 1: Simulated free surface position for a range of N T (highest order Chebyshev polynomial) after 25 periods of wave propagation for (a) kh = 2π and (b) kh = 4π.

Figure 2 :Figure 3 Figure 4 :

 234 Figure 2: Phase difference as a function of N T (highest order Chebyshev polynomial) after 25 periods of wave propagation for the simulations shown in Figure 1.

Figure 5 :

 5 Figure 5: (a) Minimum and maximum free surface elevation η from the MISTHYC simulations in comparison to the wave gauge measurements of Ting and Kirby (1994) for a regular nonlinear wave with H = 0.125 m and T = 2.0 s. (b) Beach profile used in the model simulations, showing the location of the wave gauges. In (a) and (b), the light gray shaded area (x ≥ 6 m) indicates the wave absorption zone, and the vertical gray lines indicate the wave gauge measurements shown in Figure 6.

Figure 6 :

 6 Figure 6: Simulated and measured[START_REF] Ting | Observation of undertow and turbulence in a laboratory surf zone[END_REF] free surface elevation time series at waves gauges located at: (a) x = 0.5 m, (b) x = 3.5 m, (c) x = 5.5 m, and (d) x = 6.0 m (indicated by the vertical gray lines in Figure5).

  29

Figure 7 :

 7 Figure 7: Simulated and measured[START_REF] Ting | Observation of undertow and turbulence in a laboratory surf zone[END_REF] wave skewness and kurtosis at wave gauges leading up to and in the beginning of the wave absorption zone (light gray shaded area starting at x = 6.0 m).

Figure 8 :Figure 9 :

 89 Figure 8: Bathymetry and position of wave probes in the Dingemans (1994) experiments.

Figure 10 :Figure 11 :

 1011 Figure10: Spatial evolution of the first six harmonic amplitudes (at frequencies f , 2f ,...,6f ) of the free surface elevation for case A of[START_REF] Dingemans | Comparison of computations with Boussinesq-like models and laboratory measurements[END_REF] experiments: observations (D94, circles) and MISTHYC simulation results (solid lines) .

  Cécile Raoult's PhD thesis was funded partially by the French ANRT (Association Nationale de la Recherche et de la Technologie) with CIFRE grant 2013-1024.

	Ziljema, M., Stelling, G. S., 2008. Efficient computation of surf zone waves
	using the nonlinear shallow water equations with non-hydrostatic pressure.
	Coast. Eng. 55, 780-790.
	Ziljema, M., Stelling, G. S., Smit, P., 2011. SWASH: An operational public
	domain code for simulating wave fields and rapidly varied flows in coastal
	waters. Coast. Eng. 58, 992-1012.


Yue Li

Akram Abouseif

Luigi Iannone

Damien Saucez

LISP-Views: Monitoring LISP at Large Scale

Keywords: LISP, monitoring architecture, experiment, measurement

Le protocole LISP (Locator/Identifier Separation Protocol) découple les rôles d'identifiants et de locations des adresses IP. Afin de pouvoir associer des localisateurs à des identifiants, LISP repose sur un "système de mapping". Deux plateformes expérimentales -LISP Beta Network et LISP-Lab -déploient déjà de tels systèmes de mapping LISP dans l'Internet. Toutefois, à l'heure actuelle, seule la plateforme LISP Beta Network est monitorée en reposant sur LISPmon qui échantillonne le système de mapping de LISP Beta Network une fois par jour. Afin d'accompagner la croissance de LISP, il est nécessaire de déployer un système de monitoring complet. C'est ainsi que nous proposons la solution "LISP-Views", une architecture polyvalente de monitoring à grande échelle pour LISP. LISP-Views permet d'effectuer ces mesures de manière automatique. Après avoir fait fonctionner LISP-Views pendant plusieurs mois et comparé les résultats avec ceux de LISPmon, nous parvenons à la conclusion que les résultats de LISP-Views sont plus précis et détaillés. Entre autres, LISP-Views a permis d'observer différents comportements transitoires au sein du système de mapping qui n'apparaissent pas avec LISPmon. Cette connaissance plus fine de LISP promet de futures améliorations de performances pour LISP.

Introduction

The Locator/Identifier Separation Protocol (LISP) is currently standardized at the IETF to overcome the flexibility issue of current Internet [START_REF] Farinacci | The locator/ID separation protocol (LISP)[END_REF]. It separates the IP address space into identifier locator. Then it uses a map-and-encap mechanism in order to ensure end-to-end communication. Thus, a new network entity called Mapping Distribution System (MDS) is introduced to store and distribute the bindings between EIDs and RLOCs. In this way, LISP is able to improve the Internet routing scalability, enable enhanced inter-domain traffic engineering, lightweight multi-tenant carrier-grade VPN, and VMs mobility in Data-Centers [START_REF] Saucez | Designing a deployable internet : The locator/identifier separation protocol[END_REF].

To promote the development of LISP and boost the related research, two LISP-related platforms : LISP Beta Network [lisa] and LISP-Lab platform [lisb] have been deployed and interconnected so far. Currently, a unique LISP monitoring system called LISPmon [lisc] supervises the global MDS and publishes the mapping information daily. However, it is known that the MDS is not always stable and consistent [START_REF] Li | Stability and consistency of the lisp pull routing architecture[END_REF]. We hence propose a dynamic versatile LISP monitoring architecture, LISP-Views, to overcome these limitations. LISP-Views automatically explores the whole MDS every 2 hours and stores the detailed mapping information, so to facilitate the experimenters to evaluate the LISP comprehensive performance.

We used a one-month long set of traces produced by LISP-Views to evaluate its performance and compare it with LISPmon. We show that LISP-Views is more accurate than LISPmon since it monitors all the MDS elements in parallel and provides more information.

LISP Background

LISP architecture overview

The Locator/Identifier Separation Protocol (LISP) splits the conventional IP addressing spaces into two logical sub-spaces : (i) Endpoint IDentifier (EID) : the identifier of end host and is also used to locally route within the local domain, i.e., within LISP-Site. (ii) Routing LOCator (RLOC) : refers to the interface of edge routers in the Internet topology, and is also used to route through the Internet core, i.e., between LISP-Sites. The hosts produce conventional IP packets by using their EIDs respectively as source and destination address. After the IP packets routing to edge routers, called xTR (the combination of Ingress and Egress Tunnel Router), it first checks in its mapping cache to find out the association between the destination EID and its RLOC. Otherwise, it sends a Map-Request to Map-Resolver (MR) of MDS to obtain the Map-Reply containing mapping information. It encapsulates the packets by adding the queried RLOC as the outer destination address and its own RLOC as the outer source address. Then, it sends the encapsulated packets through Internet. Once the destination xTR receives these packets, it decapsulates and forwards them to the destination host.

Based on LISP architecture, there are three types of Map-Reply can be received : (i) LISP Map-Reply : the queried IP address belongs to a LISP site and the Map-Reply contains the mapping information for this site. (ii) Negative Map-Reply : the prefix covering the queried IP address belongs to a non-LISP site, and the Map-Reply contains no mapping information. We consider these two kinds of replies as a successful query. (iii) No Map-Reply : the xTR does not receive any reply during 3 seconds. In this case, we consider the query failed.

Current LISP Deployment Status

Two experimental LISP testbeds have been deployed and interconnected so far. One is the LISP Beta Network, which is deployed in 2008 and the other one is LISP-Lab platform, which is opened since 2015. According to the latest architecture, LISP Beta Network has 12 MRs and LISP-Lab platform has 2 MRs. Besides, they also have all necessary LISP network entities, such as : xTRs, MSes and PxTRs.

LISPmon is the only monitoring platform supervising the current LISP status. It scans the whole IPv4 addressing space everyday, normally begins at 7 :00 a.m. (UTC), and queries them by sending the Map-Request to 1 specific MR only of the LISP Beta network. If this MR does not function normally, it will choose another one as a replacement, and will start the queries again, to guarantee that the receiving of LISP Map-Replies changes smoothly between two days. The results of the available mapping information are published once per day, ever since beginning of 2010.

3 Proposed Monitoring Architecture

Motivation

LISPmon was the first step towards a systematic LISP monitoring. However, it monitors the MDS just from one vantage point (VP) once per day and only queries one MR. Upon MR issues, LISPmon must be manually re-configured to monitor another MR. Yet the mapping information sometimes dramatically changes within a day, and the Map-Replies from the different MRs may not coincide at a given time [START_REF] Li | Stability and consistency of the lisp pull routing architecture[END_REF]. From such point of view, LISPmon has strong limitations since it is not able to detect the changes of mapping information within one day, and is not able to show the differences among MRs. Thus, we propose a new versatile LISP monitoring architecture, called LISP-Views, to monitor public LISP deployments, as well as to enable further performance evaluation of LISP defined by the users themselves.

Description of LISP-Views

The architecture of the proposed LISP-Views monitoring tools is depicted in Fig. 1. LISP-Views consists of several modules with different functions. The Measurement, Report, and Raw Data modules are deployed on a centralized server while, the Crawler, Sonar, and Controller modules can be deployed on several different VPs. As for now, LISP-Views is just deployed on one VP, where the monitor and centralized module are deployed together. All the modules are implemented in Python and described as follows :

Sonar module : the main module with two functions : 1. sending the Map-Requests and receiving the Map-Replies to/from all the existing MRs based on RFC 6830 [START_REF] Farinacci | The locator/ID separation protocol (LISP)[END_REF] ; 2. storing the received information in Report and Raw Data with different purposes. The IP address that Sonar uses to query to MRs is selected either from the output of Crawler, or from the recorded information in the previously produced Report. The reasons are explained in the corresponding modules. Report module : contains the collected EID-prefixes in a list, as well as a list of MRs that answered. As crawling the address space may take long time but we aim to obtain the status of MRs as frequently as possible, Sonar sends Map-Requests for the EIDs recorded in Report only to previously responded MRs. Thus, it decreases the possibility to receive no Map-Reply, so to get LISP status within a shorter time.

Raw Data module : contains all the detailed information of Map-Replies for each MR. So, it can be used to perform thorough performance analysis of MDS. Moreover, since the Raw Data is stored according to the MR, it is possible to track the performance of the MRs individually.

Measurement module : provides the composition of requested measurements (i.e., select different measurement plug-ins) based on the analysis of Raw Data and Report.

REST API : is connected to Measurement module so that the users can launch a custom experiment by setting the experiment time, the monitored MRs, and obtain the different aspects of LISP status. We are currently implementing the REST API, so for this paper the execution of measurements was performed using command lines.

Controller module : synchronizes all the modules in LISP-Views and also specifies the start and stop time of producing both Report and Raw Data. Since the inconsistency between MRs exits, the interval of producing Report and Raw Data for every MR differs from each other. Thus, the interval set in Control module should cover the slowest MR.

LISP-Views Validation

Methodology

In order to validate and evaluate the LISP-Views monitoring tool, we used raw data and report collected during one month (from 0 :00 September 4 th to midnight of October 4 th 2016), by deploying LISP-Views on one VP, which is an xTR of LISP-Lab platform. The interval to produce reports and raw data was respectively 6 and 2 hours. Unfortunately, both MRs of LISP-Lab platform just responded to the Map-Requests at the first time and then stopped. They are not able to handle large number of queries, because the Map-Requests fills request queue and are dropped. Per-se this is a success, because this bug was unknown, and the OpenLISP coders fixed the issue after our measurements. All the results of evaluation used in this paper only depend on the 6 working MRs of LISP Beta Network. The aim of the work is to validate LISP-Views by assessing if it provides at least the same results as LISPmon.

LISP-Views vs. LISPmon

Fig. 2 shows the number of LISP Map-Replies received from LISPmon and LISP-Views over time. The upper figure makes a comparison between the 2 monitors. The line with point is obtained from the daily LISPmon publications, indicating the number of LISP Map-Replies returned by one MR in an uncertain time monitoring round every day. To facilitate the comparison, we pick up the results from LISP-Views at 8 :00 a.m., which is the nearest experimental round to LISPmon. As presented in the curve with triangle symbols, LISP-Views reflects the fact of a combination of the number of LISP Map-Reply containing the different EID-prefixes from all the MRs at a fixed experiment time. In most cases, LISP-Views receives around 20 LISP Map-Replies more than LISPmon, reflects the fact that the MRs are not coherent at any time due to the existence of convergence time for synchronizing the mapping information among them. The other 5 days, where LISP-Views receives very few LISP Map-Replies, show that there were issues to normally receive LISP status at that time, whereas LISPmon always presents a more smooth trend but hides the reality.

The bottom figure of Fig. 2 depicts the number of LISP Map-Reply (in dashed line) and the overall Mapping (LISP and Negative Map-Reply, in solid line) from LISP-Views for every 2 hours monitoring round during the whole experiment. Although the number of Map-Replies presented in the figure is still the combination of results from 6 MRs when the returned EID-prefixes are different, it shows that the MRs are not quite stable. Sometimes the MRs provide less than 10 (or even 0 LISP Map-Reply), but in the next experiment round, the number is recovered. Besides, the total number of successful queries every 2 hours also presents the behavior of instability, where it generally oscillates between 400 and 500, but it occasionally drops lower than 300. Both two lines have almost same changing trend, indicating that the change of overall Map-Replies is mainly affected by LISP Map-Reply. In some monitoring rounds, where the number of LISP Map-Reply approaches 0, but we can still receive the Negative Map-Replies, means that our proposed monitoring architecture works well so that this behaviour is caused by issues in the MRs, especially for answering the LISP Map-Reply.

Conclusion

Motivated by the only LISP monitor deployed so far, named LISPmon, which records the current LISP status only from a specific MR just once per day, we propose a more dynamic LISP monitoring architecture, namely LISP-Views, so to deepen the understandings on LISP and to ease day-to-day operations and troubleshooting. After one full month of comparison with LISPmon, we find that LISP-Views provides more information by discovering more mapping information from all MRs and more complete mapping information of each MR. LISP-Views is still in the first phase, so the implementation of a REST API, the deployment on multiple VPs and the test of IPv6 behavior are still ongoing and future work.

FIGURE

  FIGURE 1: Monitoring Architecture.

FIGURE 2 :

 2 FIGURE 2: Comparison between LISPmon and LISP-Views.

Acknowledgments : Research work presented in this paper benefited support from NewNet@Paris, Cisco's Chair "NETWORKS FOR THE FUTURE" at Telecom ParisTech (http ://newnet.telecom-paristech.fr). Any opinions, findings or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of partners of the Chair.