1 DEVELOPMENT AND FIELD VALIDATION OF A 2 DH CURVILINEAR STORM IMPACT MODEL

The formulations of the 2DH process-based, nearshore morphological model XBeach were extended to allow for curvilinear grids using a finite volume approach. The formulations were tested for schematic cases such as a circular island and a field validation study was carried out for the case of storm erosion during a construction phase of the soft sea defense of the Maasvlakte-2 extension of Rotterdam port. It was found that the 2D curvilinear version performs significantly better than individual 1D runs. The overall skill of the 2D model for this case can be classified as 'reasonable' while the prediction of the erosion can be classified as 'good'.


INTRODUCTION
Modelling of the response of beaches and dunes to storms has received considerable attention over the last decades.Following the successful application of largely empirical approaches for relatively standard cases guidelines and practices have evolved in many countries based on these models (e.g.DUROS, SBEACH).However, for quite a large category of cases the underlying assumptions (uniform beach, no structures, uniform sandy material) do not hold and a more processbased approach is warranted.The recently developed XBeach model (Roelvink et al, 2009) aims to fill this gap, using a process-based approach, allowing for hard layers and graded sediments, resolving to a large extent the swash motions that make it to the dune front, and allowing a full representation of longshore variability (e.g.McCall et al, 2010).
The original implementation of XBeach uses a rectilinear, non-equidistant finite-difference grid, which is adequate for applications with a relatively short longshore extent or with a rather straight coast.However, for longer stretches of a curved coast (typical applications where you would want to use it) it becomes cumbersome to maintain the required fine resolution to resolve dune erosion processes.A curvilinear approach, as we will show, can overcome this problem and allows to include the effects of longshore gradients on coastal evolution in a smooth manner.
The direct motivation for this implementation came from a study of a storm event during the building of the Maasvlakte-II extension of the Port of Rotterdam (see Figure 7).The outer contour of the port extension had a curved, banana-like shape, which was easy to schematize on a curvilinear grid but inefficient on a rectilinear grid.

Grid setup
The new implementation utilises a curvilinear, staggered grid where depths, water levels, wave action and sediment concentrations are given in the cell centres (denoted by subscript z) and velocities and sediment fluxes at the cell interfaces (denoted by subscript u or v).In Figure 1 the z, u, v and c (corner) points with the same numbering are shown.The grid directions are named s and n; grid distances are denoted by s ∆ and n ∆ , with subscripts referring to the point where they are defined.A finite-volume approach is utilized where mass, momentum and wave action are strictly conserved.In the middle panel of Figure 1, the control volume for the mass balance is shown with the corresponding grid distances around the u-and v-points.The right panel explains the numbering of the fluxes Q and the volume V.

Mass balance equation
The mass balance reads as follows: Here, A z is the area of the cell around the cell centre, z s is the surface elevation, u u is the u-velocity in the u-point, h u the water depth in the u-point and v v the v-velocity in the v-point.The indices i,j refer to the grid number in u resp.v direction; the index n refers to the time step.

Momentum balance equation
Second, we will outline the derivation of the u-momentum balance.The control volume is given in Figure 2. It is centered around the u-point.We now consider the rate of change of the momentum in the local u-direction as follows: ( ) where V is the cell volume, u the velocity in local grid direction, Q the fluxes, ρ the density, g acceleration of gravity, , , , , the bed shear stress, wind shear stress and wave force in udirection.We consider that the outgoing fluxes carry the velocity inside the cell, u and that u in is determined at each inflow boundary by interpolation, reconstructing the component in the same direction as u.
The volume balance for the same volume reads: By multiplying the volume balance by u, subtracting it from the momentum balance and dividing the result by V we arrive at the following equation: where A is the cell area and h um is the average depth of the cell around the u-point.The procedure for the second term (the others are straightforward) now boils down to integrating (only) the incoming fluxes over the interfaces and multiplying them with the difference between u in the cell and the component of velocity in the same direction at the upwind cell.In equations (1.6) and (1.7) and (1.8) the procedure for computing the u-momentum balance is outlined.The discharges in the u-points are computed by multiplying the velocity in the u-or v-point by the water depth at that point.These discharges are then interpolated to the borders of the control volume around the u-point.The difference α ∆ in grid orientation between the incoming cell and the u- point is computed and used to compute the component of the incoming velocity in the local u-direction, from the left and right side of the control volume.
q 1 1 q = q +q , q = -q +q 2 2 , cos sin , cos sin The same is done for the top and bottom of the control volume, based on the discharges in v-direction: Finally, the advective term in the momentum balance is given in equation (1.8).

Time integration scheme
The time integration of the mass and momentum balance equations is combined in an explicit leap-frog scheme, as depicted in Figure 3.The velocities (in the '-' points) are updated using the momentum balance, the water levels are updated using the mass balance.The water level gradients influence the momentum balance and the velocities and derived discharges affect the mass balance.Because of the leap-frog scheme these influences are always computed at the half time step level, which makes the scheme second order accurate.

Figure 3 Leap-frog time integration scheme
Using this straightforward finite volume approach, complicated transformations of the equations are avoided and the solution scheme remains transparent.It is also completely compatible with the original rectilinear implementation and is even slightly more efficient.

Wave action balance
The time-varying wave action balance solved in XBeach is as follows: Where E is the wave energy or wave action, Cg is the group velocity, C ϑ the refraction speed in theta- space and Sink refers to effects of wave breaking and bottom friction.Again, the advection terms are the only ones affected by the curvilinear scheme so we will discuss their treatment in detail.the control volume is the same as for the mass balance.In equation(1.10) the procedure to compute the wave energy fluxes across the cell boundaries is outlined.All variables should also have an index itheta referring to the directional grid, but for brevity these are omitted here.
The component of the group velocity normal to the cell boundary, at the cell boundary, is interpolated from the two adjacent cell center points.Depending on the direction of this component, the wave energy at the cell boundary is computed using linear extrapolation based on the two upwind points, taking into account their grid distances.This second order upwind discretization preserves the propagation of wave groups with little numerical diffusion.
, , (1.10) The other three fluxes are computed in a similar way; for brevity we will not present all formulations.The time integration is explicit and the same as in the original implementation.The advection in u-and v-direction is computed simply by adding the four fluxes and dividing by the cell area.This procedure guarantees conservation of wave energy.

Advection-diffusion equation for suspended sediment
The advection-diffusion equation for suspended sediment is the basis for the sediment transport computations in XBeach.The partial differential equation to solve is: Here c is the depth-averaged concentration, c eq the equilibrium concentration, Ts a typical timescale proportional to water depth divided by fall velocity.As is often done to increase robustness, we treat the erosion term explicitly but take an implicit scheme for the sedimentation term: This can be rewritten as: The sediment transport gradient is discretized in a similar way as the mass balance: , The sediment transports in the u-points contain an advective term, a diffusive term and a bed slope term: Here u rep,s is a representative velocity for suspended transport, which contains contributions due to return flow, wave skewness and wave asymmetry; D c is a horizontal diffusion coefficient and f slope a coefficient.In discretized form the expression for the suspended transport in the u-point is: (1.17) The concentrations in the u-points are computed with a θ -method, where 1 θ = means a fully upwind approximation, and 0.5 θ = a central scheme.In practice, we mostly use the upwind approximation for its robustness.

( ) ( )
The erosion and deposition terms, which may also be used in the bed updating, are finally computed from: , , , , i,j,n i,j h eq s , , , 1 i,j,n+1 i,j h s The evaluation of the bedload transport takes place in the same way as in the previous versions of XBeach, except for the fact that the directions are taken in local grid direction, and will not be repeated here.

Bed updating
Two alternative formulations are available for the bed updating: one where the bottom changes are computed based on the gradients of suspended and bed load transport, equation (1.20), and one where the changes due to suspended transport are accounted for through the erosion and deposition terms, equation (1.21).
( ) ( ) In both cases MF is the morphological factor used to accelerate morphological changes.In the first case, the sediment in the bottom is conserved in all cases, but changes in the amount of sediment in the water are not considered; one can also say that the sediment in suspension is added to the bottom sediment.In the second case, the storage of sediment in the water is accounted for, but will be distorted in cases of high MF.Since under most circumstances the real effect of the storage in the water phase is small we prefer the first formulation which guarantees mass conservation in the bottom.
In view of the limitations of this paper we have omitted many details, but we believe that the main ideas have been put on paper with this overview.Details will be found in manual and implementation documents, as well as the source code, available on www.xbeach.org.

Implementation test
To test the correct implementation of the curvilinear scheme we carried out some principle tests.One of these, which also shows a potentially interesting application, is that of a circular island with a reef.In Figure 4 the bathymetry is shown; the offshore depth is a uniform 2m deep shelf, while the reef has a steep slope up to a horizontal reef flat area at 0 m relative to MSL, followed by a more mildly sloping beach.The test case is run with regular wave groups from the south with a wave height varying between 0 and 2 m, and a still water level of 0.5 m.In Figure 5 a snapshot is shown of the wave height, clearly showing that the wave group crests stay parallel on the horizontal part and propagate without losing energy until they dissipate and refract on the reef edge.In Figure 6 we see the time-averaged wave-driven current and the Hrms wave height, which both show a regular and logical pattern correctly aligned with the island orientation.
This test case offers many more fascinating features, such as long wave resonance, but unfortunately we cannot discuss them all here.Suffice it to state that we tested the curvilinear scheme on a number of tests like these and furthermore ran all the standardized tests in our test bed (comprising over 60 cases) to verify that also for rectilinear cases the scheme still performs as before, so that we can have some confidence in the results presented in the next section, pertaining to a real field validation.

Storm event
During a storm event on 3-5 September 2009 considerable scarping of the steep profiles took place, where the eroded sediment was not just deposited seaward but to a large extent was moved alongshore.Wave conditions along the Maasvlakte 2 sea defence has been predicted by PUMA using the measured wave data from Europlatform using the SWAN wave prediction model.The Europlatform station is located approximately 30 km offshore from Hoek van Holland at 32 m water depth.With this constructed wave transformation matrix, the measured Europlatform wave conditions can be transformed to a number of locations along Maasvlakte 2 sea defense.A morphological determining storm condition is defined by PUMA when the significant wave height on the location is above 2 meters.Such a storm occurred on 3-5 September 2009.During the storm, the significant wave height at the defined locations at about 500 m in front of the island varied from 1.8 m to 3.2 m (2.4 m average) while the peak wave period varied from 6 to 7.5s (6.7 s average).Further, the wave direction firstly came from the west direction (259º) and slowly changed its direction to the north-west (305º) at the end of the storm as shown in Figure 10.Water level time series were obtained from a number of observation points near the MV2 site, as given in Figure 9. Clearly, the data at Europlatform show a smaller amplitude and surge level as it is closer to an amphidromic point and misses part of the windinduced setup.

Data analysis
In Figure 10 the main events during the storm are summarized.The storm turns fro WSW to NW, leading to strong longshore transports to the North in the beginning and to the South towards the end.Over most of the length of the site considerable scarping took place, but there is clearly not a COASTAL ENGINEERING 2012 10 cross-shore balance: much of the sediment eroded from the beach/dune is deposited at both ends of the island, and additionally some sand is deposited landward by overwash processes.At the same time, it must be noted that during the storm production continued, and the locations of the production coincide with the areas of deposition; however, this explains only a relatively small part of the deposition, and in comparisons between observations and model results these volumes are corrected for.Typical cross-shore profile behavior is shown in Figure 11, indicating that in the surveys most of the eroded material is lost from the profile; in contrast, in a 1D XBeach simulation, sediment volume in the profile is conserved, and consequently much less erosion is predicted since the deposited sand helps to dissipate the incoming wave energy.

Model simulations
To compare the effects of 1D vs 2DH simulations and to evaluate the performance of the curvilinear grid scheme, two sets of simulations were carried out: one where 1D simulations were performed for each of the transects shown in Figure 8, and one using the new 2D curvilinear version.

COASTAL ENGINEERING 2012
11 A curvilinear grid was set up with resolution of 6 to 20 m in cross-shore and 17 to 50 m in longshore direction; see Figure 12, left panel.Observed timeseries of water levels and wave conditions were imposed and a simulation was carried out over a period of three days, using default sediment transport settings as reported in Roelvink et al., 2009.In Figure 12, right panels, two detailed snapshots of the current pattern are shown, indicating the concentrated longshore current both in northerly direction (in the northern part) and in southerly direction (in the southern part).The effect of these currents on the sedimentation/erosion pattern is clearly visible in Figure 13 where the deposition patches on either side of the island are clearly visible.The overwash patch as observed is hardly present in this simulation, which is based on the Europlatform water level data.As we will see in the sensitivity tests, a more realistic water level time series improves this.There are differences in the size and exact location of the sedimentation patches at either end of the island, which can be ascribed to processes such as wind and tide that were neglected in this simulation, as well as to the sand production for which the observations shown here were not corrected.The pattern of erosion in particular is very well represented, as is also shown in Figure 14.This important result is less sensitive to the problems related to the prediction of the deposition: it is easier to predict the location of erosion since this has to be at the most exposed and/or steepest locations, whereas the sedimentation has a wider range of possible locations, depending on subtle processes such as the mean current, settling lag and such.In Figure 15 the 1D and 2D results are compared for two transects, one in the northern section and one in the middle.Clearly, there is a large improvement in the profile change for the 2D model, obviously related to the longshore transport gradients, but in part also due to the effect of the longshore current on sediment resuspension.In the right panel, the 1D model produces some overwash deposit in contrast with the 2D result.This is because both use the Europlatform water level data, but the infragravity waves are exaggerated somewhat in 1D and therefore still make it to this location.In Figure 16 the Brier Skill Score per survey line is shown, both for the 1D and the 2D run.There is an obvious improvement for the 2D run relative to the 1D results.Results at the ends have no skill as the exact location of the deposition is not captured well.Also, around KP 5200 the skill is less due to an underestimation of the overwash processes.The sensitivity of the overwash deposit to the water level scenario is shown in Figure 17, where a detailed view of this deposit is shown for the base case and the more realistic case with water levels measured close to the site.Clearly, for areas such as this, where it is 'to be or not to be' for the overwashing, peak water level matters a lot.The overall performance of this 2D model is shown for the different water level scenarios.Clearly the more realistic water levels in scenarios 2-4 lead to better performance.The skill of the erosion can be classified as 'good', which is encouraging since this is often the most important issue in cases where storm impact is considered.

CONCLUSIONS
The formulations of the 2DH process-based, nearshore morphological model XBeach were extended to allow for curvilinear grids using a finite volume approach.The formulations were tested for schematic cases such as a circular island and a field validation study was carried out for the case of storm erosion during a construction phase of the soft sea defense of the Maasvlakte-2 extension of Rotterdam port.It was found that the 2D curvilinear version performs significantly better than individual 1D runs.The overall skill of the 2D model for this case can be classified as 'reasonable' while the prediction of the erosion can be classified as 'good'.
Put in a wider context, we believe that longshore currents and non-uniform transport are important in many cases and existing 1D models in general do not capture this.The curvilinear version of XBeach brings efficient simulations for mmore complex, non-uniform cases within easier reach.

Figure 1
Figure 1 Location of staggered grid points (left panel); definition of grid distances (middle) and terms in volume balance (right)

Figure 2
Figure 2 Control volume u-momentum balance and definition of fluxes

Figure 5
Figure 4 Bathymetry of circular island with reef test Figure 6 Hrms wave height and wave-driven current pattern

Figure 7
Figure 7 Final design of Maasvlakte-2 (left panel) and building phase during storm.

Figure 8
Figure 8 Survey lines (left panel) and maximum elevation per line (right)

Figure 9
Figure 9 Water level gauging stations and observed time series

Figure 10
Figure 10 Situation sketch of storm event and consequences (left panel) and areas where sand production continued during storm (right panel)

Figure 11
Figure 11 Typical profile behavior during storm; comparison between surveyed profile change and predictions by 1D XBeach

Figure 12
Figure 12 Curvilinear grid (left panel) and snapshots of the current pattern at northern and southern ends (right)

Figure 15
Figure 15 Comparison of 1D and 2D XBeach profile changes with observations

Figure 16
Figure 16 Brier Skill Score per survey line, comparison between 1D and 2D Xbeach