MULTI-DIMENSIONAL ERROR ANALYSIS OF NEARSHORE WAVE MODELING TOOLS, WITH APPLICATION TOWARD DATA-DRIVEN BOUNDARY CORRECTION

As the forecasting models become more sophisticated in their physics and possible depictions of the nearshore hydrodynamics, they also become increasingly sensitive to errors in the inputs, such as errors in the specification of boundary information (lateral boundary conditions, initial boundary conditions, etc). Evaluation of the errors on the boundary is less straightforward, and is the subject of this study. The model under investigation herein is the Delft3D modeling suite, developed at Deltares (formerly Delft Hydraulics) in Delft, the Netherlands. Coupling of the wave (SWAN) and hydrodynamic (FLOW) model requires care at the lateral boundaries in order to balance run time and error growth. To this extent, we will use perturbation method and spatio-temporal analysis method such as Empirical Orthogonal Function (EOF) analysis to determine the various scales of motion in the flow field and the extent of their response to imposed boundary errors.


INTRODUCTION
Accurate descriptions of the nearshore wave, hydrodynamic and sediment transport and morphologic processes are necessary for many civilian and military activities in the nearshore.In particular, forecasting environmental conditions at a given area is of great importance for naval exercises and operations.This capability has progressed beyond simple reduced-dimension models (e.g., Navy Standard Surf Model; Earle 1989) to more sophisticated comprehensive three-dimensional hydrodynamic models (e.g., SHORECIRC; van Dongeren and Svendsen 2000; Delft3D; Lesser et al. 2004), with the realization that strong spatial and temporal non-homogeneities will greatly affect the nearshore hydrodynamic environment.
However, sites for operations in coastal and nearshore areas are usually poorly sampled, and input data are often obsolete (bathymetry) or of low resolution (input wave and current conditions).The effect of these characteristics of input data on increasingly sophisticated forecast models is unclear, but is necessary to know so that reasonable confidence limits can be placed on the results.This is especially true of boundary conditions.Additionally, the forecasts must be performed expediently; however many modeling practices developed to mitigate the potential effects of unknown boundary errors are not well suited to speedy forecasts.We propose to develop methods to evaluate the effect of boundary errors and help optimize the typical run configuration to reduce wasteful calculation required to handle open boundaries.In this study we investigate the effect of lateral grid extensions on the model predictions, and will use two test areas (Duck, NC, and La Jolla, CA, both in the USA) for this evaluation.

GENERAL MODEL APPROACH
The model under investigation herein is the Delft3D modeling suite, developed at Deltares (formerly Delft Hydraulics) in Delft, the Netherlands.The model suite consists of a wave (SWAN), hydrodynamic (FLOW) and sediment model chain which can be activated in standalone or coupled modes.FLOW is capable of simulating flow in a wide range of conditions from coastal-scale tidal flows to wave-induced nearshore flows.It is a terrain-following three-dimensional model, and can be used in curvilinear or Cartesian coordinates.Nonhydrostatic flows and shoreline wetting and drying are also selectable options in the model.SWAN is a spectral wave model which simulates the generation of surface waves due to wind and the propagation of wave energy over arbitrarily-varying bathymetry.Wind-wave generation, whitecapping, wave-current interaction, deep and (parameterized) shallow water nonlinearity, bottom friction and wave breaking dissipation are all represented in the model.

Boundary Conditions
In order to allow flow to leave the lateral boundaries with no artificial circulation, Roevink and Walstra (2004) developed a Neumann lateral boundary condition for the FLOW model.This condition assumes that the longshore gradient of the mean sea surface is zero at the lateral boundaries of the grid; the remaining dynamic variables would reach their "natural" values at the boundaries, and the resulting flow would be free of boundary-induced effects.
Since the Neumann boundary conditions require longshore uniform conditions at the lateral boundaries, coupling of the wave and hydrodynamic models requires care at these boundaries.The optimum method of doing this is by laterally extending the grid for the forcing wave model beyond the boundaries of the hydrodynamic model (Figure 1), thus keeping irregularities in the forcing far away from the domain of interest.A major disadvantage of this approach is the increased computational time required for the wave model, which has typically been the bottleneck in the overall calculation time, particularly for stationary (equilibrium) conditions.We wish to investigate the effect of truncating this extraneous lateral extension has on the run time and the errors in the model.
Herein we define the extension ratio A/B, where A is the length of the extension of the SWAN model grid on either side of the FLOW grid, and B is the overall longshore extent of the FLOW grid (Figure 1).The higher the ratio A/B, therefore, the longer the overall extent of the wave model grid is.

Duck94 Field Experiment
Duck94 field experiment was conducted in August and October near the Army Corps of Engineers' Field Research Facility pier located in Duck, North Carolina (shown in Fig. 2).A bathymetry survey for the so called mini-grid, which is marked by the box, was conducted daily during the intensive study period in October.
The SWAN model grid consisted of 85 columns and 160 rows with grid spacing of 10 m in the cross-shore direction (x), and 15 m along the shore (y).The smaller FLOW grid lay inside the SWAN domain with 85 by 80 grid points with grid spacing of 10 m in x and 15 m in y.In order to satisfy the Neumann boundary condition, longshore uniform areas are extended by a few grid points into the FLOW domain on both side boundaries.Directional spectrum from the 8-m array (pressure gauges located on the 8-m contour about 900 m offshore) was used as offshore boundary to the wave model.The dates selected for model computation were based on the availability of mini-grid bathymetry surveys.Model simulation was conducted for every 3 hours and 132 cases were run.

Nearshore Canyon Experiment (NCEX)
The data from the site of the Nearshore Canyon Experiment (NCEX), conducted near Scripps Institution of Oceanography (La Jolla, CA) in 2003, was used in this study.In addition, there are active relevant websites run by the Coastal Data Information Program (CDIP), which serve as the source of input data for our study.The bathymetric data was taken from surveys conducted during the experiment, while the offshore wave spectrum forcing for the model was made available through the CDIP web site for January 2010.The steep topography at the canyons could be expected to cause significant variation in nearshore wave energy; in particular, complex refraction effects are likely to focus/defocus wave energy at various locations alongshore, leading to strong longshore variation of energy.This expected spatial variation in sensitivity of the overall model results to boundary errors was one of the reasons for selecting this as the study area.
Because forcing from the CDIP buoy is located about 8 km offshore, the SWAN model was first run over a larger domain (shown in Fig. 3 left), and the wave spectrum results along the offshore boundary of an approximately 1 km by 2 km grid (shown in Fig. 3, right) was written out to be used as the boundary condition for the SWAN runs over this smaller area of interest.
A computational grid resolution of 5 m in the cross-shore direction (x), and 15 m along the shore (y) was used in SWAN.We adopted higher resolution for FLOW, which is 2 m in x and 5 m in y.This was chosen to be sufficiently high to reduce the effect of numerical artifacts on the model results, and thereby maximize the impact of the input conditions.Model runs were conducted for every hour and 216 cases were run.
All SWAN runs were stationary for both Duck94 and NCEX cases in this study; waves do not change during the one hour flow computation.In addition, all model runs were made without wavecurrent interactions.

Comparison of Field Observations and Model Experiments
The two-dimensional Delft3D model has been shown to be capable of simulating nearshore hydrodynamic processes such as wave, winds, flow, tides and so forth, over complex bathymetries in which case a one-dimensional nearshore model (e.g. the Navy Standard Surf Model) would be inappropriate (Morris 2001).In this study we will use data from the Duck94 experiment (Duck, NC) to establish the accuracy of the basic mode.
As mentioned above, the Delft3D model was run for 132 cases of Duck94 data.The scatter plots for both significant wave height and longshore current for ratio A/B equal to 50% are presented in Fig. 4.This analysis on the 'best' modeled field (A/B = 50%) yields our "Golden Standard" by which we can evaluate further possible degenerative effects of imposed errors defined relative to the standard on the boundaries.
The skill statistics for wave height and longshore current results are summarized in Table 1 and Table 2 respectively, where R represents the linear correlation coefficient, slope is the slope of the linear regression line (solid line), and N is the number of observations.Additionally the root mean square error (RMSE) and mean absolute gross error (MAGE) are listed.It is clear that the model results agree reasonably well with measurement.Obviously, both the RMSE and MAGE grow as the grid extension is decreased.It should be noted that the higher the ratio A/B, the more computational time is required of the model simulation.
With a reasonable degree of confidence in the model skill, we now can investigate the effect of errors on the possible deterioration of model performance.However, we have only confirmed the model performance in pointwide comparisons to data.Is it possible to systematically examine and evaluate the model performance in terms of spatial variation?Are there any certain spatial patterns we can track as the boundary errors increase?To address these questions, we will further introduce and interpret the Empirical Orthogonal Function analysis in the next part.

ERROR ANALYSIS
While it is always possible to compare model results to data to determine the optimum grid configuration, the point-to-point comparison does not offer any information on the spatial response of the model to the forcing, whether correct or erroneous.Is a poor data-model comparison due to a slight spatial shift of a highly-variable field by the model, or is it due to a complete inadequacy of the model physics or numerics?Figure 5 shows the flow field near the boundary for several values of ratio A/B.As the ratio becomes smaller, irregularities in the forcing field from the waves begins to effectively pollute the interior of the hydrodynamic model domain.This is demonstrated by the development of an eddy in the domain, which strengthens as A/B reduces.Despite this, the patterns elsewhere seem reasonably similar, indicating that reliance on a data-model comparison alone may not raise issues of concern elsewhere in the domain.

Empirical Orthogonal Functions
Empirical Orthogonal Function (EOF) analysis offers a glimpse into the spatial and temporal nature of modeled field.As a start, we used EOF analysis to decompose time series of spatiallydistributed current fields into separate, linearly independent modes which cascade in variance with increasing mode number.In general, the analysis herein addresses the two-dimensional behavior of the model response in the face of lateral boundary condition errors.
As outlined by Davis (1976), two advantages of a statistical EOF description of the data are: 1. EOFs provide the most efficient method of compressing the data (a very few empirical modes can be used to describe the fundamental variability in a large data set); and 2. EOFs may be regarded as uncorrelated modes of variability of the data field.
We will perform this analysis on a "best" modeled field, yielding a standard by which we can evaluate the effects of imposed errors on the lateral boundaries, such as seen in Figure 5.The separation into modes will allow a comparison to determine which modes are most affected by these errors, which in turn helps determine the relevant scales of motion most vulnerable to errors.The result will allow value judgments concerning the important modes to be made and allow for the establishment of a balance between the desired level of accuracy and expediency.

Model Output
We used the model results as data for EOF analysis.Eleven stations, marked with stars in Figure 6, were defined in the FLOW domain of NCEX simulations.Water depth increases from right to left on this bathymetry map. Figure 7, 8 and 9 show the time series of cross-shore current, longshore current and wave height respectively with 50% grid extension ratio for eleven stations.Top panels (of Figure 7, 8 and 9) present daily time series from January 1st, 2010 to January 31st, 2010.In order to obtain more data with high variance for EOF analysis, model runs with hourly initial wave spectrum (instead of a single representative daily spectrum) were also conducted, focusing on days with relatively high measured wave energy.The selected hourly time series spanned January 9th -January 17th is marked by the box in the middle panels.We zoomed into the boxes to see the hourly time series features more clearly, as seen in the bottom panels.The top two panels in Figure 10 show hourly time series of longshore current for both 50% ratio (left) and 25% ratio (right) cases.Straightforward comparison between these two time series is not possible; we thus calculate and display a histogram for each variable.These are shown in the lower panels of Figure 10.It is apparent that the introduction of boundary errors affects the model's ability to simulate southward (negative) longshore currents.

EOF Results
Figure 11 depicts the mean velocity, the spatial and temporal behavior for EOF modes 1 -3 based on the NCEX case with the extension ratio A/B equal to 50% (standard) and 25% respectively.First of all, in the mean velocity field, the mean velocity of the standard is mainly in the northwest direction, but it is towards west for 25% ratio case.EOF mode 1 of our standard clearly shows a strong slablike flow structure in the center, compared with the vortical flow structure in 25% case on the right.A shear flow structure can be observed in the area marked with red circle in EOF mode 2 of the standard, while we obtain an eddylike flow structure in the same position of mode 2 of 25% ratio case.Similarly, we can tell the difference in mode 3 between our standard and 25% case as well.In general, the feature of the first few modes of standard does be changed by grid extension reduction, i.e. the imposed boundary errors affect the spatial and temporal nature of modeled field of the standard.

Swirl Strength
We note that it is important to know the characteristics of eddy structures and their distribution in the model output.Many methods can be applied to identify vortical structures two-dimensional data.Herein we choose the local swirl strength criterion, which was shown (Adrian et al. 2000) to be the most effective at identifying the full range of vortices.The local swirl strength can be calculated by building an equivalent two-dimensional velocity gradient tensor D given by where u and v are respectively the instantaneous velocity vectors in the cross-shore x and longshore y coordinate directions.The velocity gradients are calculated from EOF data using a two-point central difference method (Raffel 1998).The local swirl strength, represents the local frequency of rotation, is the imaginary part of the complex conjugate eigenvalues of the tensor D. It is clear in the Figure 12 that the swirl strength identifies the vortexes which are clearly visible in the data.

CONCLUSIONS
Duck94 data were used for evaluating Delft3D performance.In general, Delft3D has been shown to be accurate and robust in predicting nearshore flows.It should be noted that the higher the ratio A/B, the more accurate the representation of the Neumann lateral boundary condition, however, the more computational time is required of the model simulation.
The method of Empirical Orthogonal Functions and Swirl Strength were used to determine the various scales of motion in the flow field and the extent of their response to imposed boundary errors.A limited number of the first few EOFs, those with greatest eigenvalues, can be used to describe the fundamental variability in a very large data set.
Inadequacies in the model setup are reflected in the EOF structures.In the future, it should be possible to further quantify the effect of imposed error by completing inter-mode comparisons among different error cases.

Figure 1 .
Figure 1.Sketch of grid extension (left) and definition of the extension ratio A/B (right).

Figure 3 .
Figure 3. Bathymetry of over the extent of the larger WAVE domain and selected area of interest.

Figure 4 .
Figure 4. Scatter plots for wave height comparison (left) and longshore current comparison (right) with 50% grid extension.

Figure 5 .
Figure 5.The effect of reduction of the lateral extent of SWAN grid on the simulation of wave-driven currents over Duck94 bathymetry.Shoreline is on the left of each plot of velocity vectors.a) A/B=50%; b) A/B=18.8%;c) A/B=5%.See Figure 1 for grid extension ratio definition.

Figure 6 .
Figure 6.Layout of NCEX stations (coordinates display on the right) for Delft3D output.

Figure 7 .
Figure 7. Time series of cross-shore current of NCEX stations with 50% grid extension ratio.

Figure 8 .
Figure 8.Time series of longshore current of NCEX stations with 50% grid extension ratio.

Figure 9 .
Figure 9.Time series of wave height of NCEX stations with 50% grid extension ratio.

Figure 12 .
Figure 12.Velocity field from EOF data (top) and Swirl Strength estimate (bottom).