Data assimilation within the Advanced Circulation (ADCIRC) modeling framework for the estimation of Manning’s friction coefficient

Ocean Modelling 76 (2014) 43–58

Contents lists available at ScienceDirect

Ocean Modelling journal homepage: www.elsevier.com/locate/ocemod

Data assimilation within the Advanced Circulation (ADCIRC) modeling framework for the estimation of Manning’s friction coefficient Talea Mayo b,⇑, Troy Butler a, Clint Dawson b, Ibrahim Hoteit c a

Department of Mathematical and Statistical Sciences, University of Colorado Denver, Denver, CO, United States Institute for Computational Engineering and Sciences, The University of Texas at Austin, Austin, TX, United States c King Abdullah University of Science and Technology, Thuwal, Saudi Arabia b

a r t i c l e

i n f o

Article history: Received 1 July 2013 Received in revised form 13 November 2013 Accepted 11 January 2014 Available online 27 January 2014 Keywords: Bottom stress Manning’s n Data assimilation Parameter estimation Shallow water equations SEIK filter

a b s t r a c t Coastal ocean models play a major role in forecasting coastal inundation due to extreme events such as hurricanes and tsunamis. Additionally, they are used to model tides and currents under more moderate conditions. The models numerically solve the shallow water equations, which describe conservation of mass and momentum for processes with large horizontal length scales relative to the vertical length scales. The bottom stress terms that arise in the momentum equations can be defined through the Manning’s n formulation, utilizing the Manning’s n coefficient. The Manning’s n coefficient is an empirically derived, spatially varying parameter, and depends on many factors such as the bottom surface roughness. It is critical to the accuracy of coastal ocean models, however, the coefficient is often unknown or highly uncertain. In this work we reformulate a statistical data assimilation method generally used in the estimation of model state variables to estimate this model parameter. We show that low-dimensional representations of Manning’s n coefficients can be recovered by assimilating water elevation data. This is a promising approach to parameter estimation in coastal ocean modeling. Ó 2014 Elsevier Ltd. All rights reserved.

1. Introduction Coastal and estuarine systems are home to more than 50% of the human population as well as many biological species (Yanagi, 1999). As a result, coastal inundation during extreme events such as hurricanes and tsunamis poses a major threat to the livelihood and vitality of these regions. Coastal ocean modeling plays a crucial role in forecasting of such inundation and is also used to model water surface elevations and currents under more moderate conditions (e.g. tides). Thus, costal ocean modeling is critical to the missions of conservation, emergency and/or economic planning for these systems. The northern coastline of the Gulf of Mexico is of particular interest for a variety of reasons, including the conservation of wetlands in Louisiana and the economic activity at its many ports, e.g. in Galveston Bay. Coastal ocean models numerically solve the shallow water equations (SWEs), which model the flow of water for processes where the horizontal length scales are large relative to the vertical length scales. The SWEs are obtained by the depth integration of the incompressible Navier–Stokes equations, assuming hydrostatic pressure. The result is a first-order hyperbolic continuity equation

for water elevation coupled to the momentum equations for horizontal depth-averaged velocities. We summarize these equations in Appendix A. Accurate estimation of water elevations and currents requires accurate estimation of many model parameters, which are either measured directly, defined empirically, or inferred from (sometimes sparsely collected) data. Of particular importance is the Manning’s n coefficient of roughness, which enters into the SWEs through the bottom stress components in the momentum equations (Kennedy et al., 2011). These bottom stress components, sbx and sby , are defined using a linear or quadratic drag law through the coefficient K slip :

sbx K slip Q x ¼ ; q0 H sby K slip Q y ¼ : q0 H In this work we utilize a quadratic drag law, letting K slip ¼ cf juj. The coefficient cf can be determined using the Manning’s n formulation,

cf ¼ ⇑ Corresponding author. E-mail addresses: [email protected] (T. Mayo), [email protected] (T. Butler), [email protected] (C. Dawson), [email protected] (I. Hoteit). http://dx.doi.org/10.1016/j.ocemod.2014.01.001 1463-5003/Ó 2014 Elsevier Ltd. All rights reserved.

gn2 H1=3

;

ð1Þ

where n is the Manning’s n coefficient of roughness. The Manning’s n coefficient is a highly variable spatial parameter dependent on the

44

T. Mayo et al. / Ocean Modelling 76 (2014) 43–58

surface characteristics of the seabed (e.g. Fig. 1). It cannot be measured directly, and there is no exact method for its selection; values are often assigned from tables defining well established empirical estimates of common land classifications (e.g. Table 1). In general, the Manning’s n coefficient is specified at the vertices (nodes) on a given discretization of a physical domain. From these nodal values, a continuous, piecewise linear representation of bottom friction can be constructed. The nodal values are determined from data (e.g. satellite images used to define land classifications) that may be incomplete or contain data from different, inapplicable time periods (e.g. see Medeiros et al. (2012)). Furthermore, data is often collected on a sub-grid scale and an upscaling procedure (e.g. local spatial averaging) must be used. As a result, the Manning’s n coefficients used in coastal ocean models often contain large amounts of uncertainty. The uncertainty in the Manning’s n coefficient on any given domain leads to uncertainties in common quantities of interest computed from the solution of a coastal ocean model (e.g. maximum water elevations and currents). Moreover, the solution of the model and thus the quantities of interest are often highly sensitive to changes in these parameters. In this paper, we consider a specific numerical model used to solve the SWEs and implement a statistical data assimilation methodology for estimating the value of Manning’s n coefficients defined in coastal lagoons and estuaries. Uncertainties in spatially variable parameters can be modeled by representing the parameter as a particular realization of a random field (Wiener, 1938; Ghanem and Spanos, 2003). In this work, we focus on estimating low-dimensional representations of the Manning’s n field that model macro-scale uncertainties/perturbations of the true field. This is particularly appropriate in two cases: the common case where we are provided with a grid containing the fine-scale features of the Manning’s n coefficient determined from prior data collection efforts, and we seek to update the field subject to large-scale piecewise-smooth perturbations (e.g. due to sediment transport or erosion over a significant time period and/or due to an extreme event), and the case where we model Manning’s n coefficients for the ocean floor, which are far less variable than the coefficients over land. Moreover, since only a modest number of random samples are required to accurately estimate second-order statistics in low-dimensions, these representations of perturbations to a parameter field are well-suited for low-rank non-intrusive data assimilation methodologies such as the singular evolutive interpolated Kalman (SEIK) filter (Pham, 2001; Hoteit et al., 2002; Triantafyllou et al., 2003) used in this paper. We explicitly use

Table 1 Values of the Manning’s n coefficient for various surfaces (Bunya et al., 2010). The left column describes the land characteristics and the right column is the associated empirically defined Manning’s n coefficients. Open water Ice/snow Pasture Commercial Bare rock/sand Gravel pit Fallow Transitional Deciduous forest Evergreen forest Mixed forest Mixed forest Shrub land Grassland Low residential High residential Row crops Small grains Recreational grass Woody wetland Herbaceous wetland

0.020 0.022 0.033 0.050 0.040 0.060 0.032 0.100 0.160 0.180 0.170 0.170 0.070 0.035 0.120 0.121 0.040 0.035 0.030 0.140 0.035

predictive model estimates and model sensitivity to analyze the quality of solutions to parameter estimation problems. It is important to note that by model sensitivity, we mean the sensitivity of the data computed from the solution (i.e. state vector) of the model with respect to perturbations in the model parameter. In this work, our focus is water elevation data. In practice, this type of data is fairly straightforward to accurately record using existing observation networks in contrast to other types of data such as fluxes or time of inundation. We explore the limitations of using this data and present cases where the water elevation data proves insufficient for accurate estimation of model parameters despite improvements in water elevation forecasts. This paper is organized as follows. In Section 2, we discuss the specific numerical model used to solve the SWEs. In Section 3, we discuss the basic components of statistical data assimilation, the reduced order SEIK filter, and a method of estimating model parameters sequentially based on discrepancies in model forecasts and observable data. In Section 4, we describe the framework and notation for the numerical observations and parameterized fields of Manning’s n values. In Section 5, we present several case studies

Fig. 1. Values of the Manning’s n coefficient along the upper Texas coast.

T. Mayo et al. / Ocean Modelling 76 (2014) 43–58

including estimating Manning’s n on a computational domain representing Galveston Bay. Concluding remarks follow and the appendices contain brief details pertaining to the equations of the coastal ocean model and data assimilation methodology that are otherwise well-established in the literature. 2. The Advanced Circulation model The Advanced Circulation (ADCIRC) model (Luettich and Westerink, 2004) is a widely used coastal ocean model (Bacopoulos et al., 2012; Blanton et al., 2012; Dietrich et al., 2012; Hagen and Bacopoulos, 2012). In this model, the continuity equation (A.1) is replaced by what is known as the generalized wave continuity equation (GWCE) (Kinnmark, 1986; Lynch and Gray, 1979). The GWCE is a second-order, hyperbolic equation, obtained in several steps. See Appendix B for more details. Use of the GWCE reduces the occurrence of spurious oscillations often associated with the numerical solution of the original form (Dawson and Proft, 2004). Together, the GWCE (B.1) and the momentum equations (A.2) define the modified form of the SWEs solved by ADCIRC. In the ADCIRC model, the GWCE and the momentum equations are discretized in space by a first-order continuous Galerkin finite element method. The elements that make up the computational domain are triangular. The time derivatives are approximated using centered finite differences in the GWCE (B.1) and forward differences in the momentum equations (A.2). The ADCIRC model has been verified and validated for hurricanes in the Gulf of Mexico using hindcast studies, which compare common quantities of interest computed from the model solution (e.g. maximum water elevations at various spatial locations) to data from actual storms. Storms used for hindcasting include Hurricanes Betsy (1965), Ivan (2004), Dennis (2004), Katrina (2005), Rita (2005) (Bunya et al., 2010; Dietrich et al., 2010; Westerink et al., 2008), Gustav (2008) (Dietrich et al., 2011), and Ike (2008) (Kennedy et al., 2011), and these studies have guided decision-makers in coastal development planning and coastal hazard mitigation. The accuracy of the hindcasts depends on highly accurate model inputs and parameters (e.g. bottom topography, coastal elevation, and friction characteristics) and accurate data (Kennedy et al., 2011). The past decade has seen an extensive deployment of sensor networks in coastal regions, allowing the real-time collection of data such as water elevations during a storm. Massive data collection and analysis has also provided snapshots in time of the state of coastal features, such as bottom topography and land characteristics. We propose that data collected under typical weather conditions (i.e. when the predominant forcing is tides) be assimilated with model forecasts to accurately estimate large-scale perturbations to parameter fields of particular importance during an extreme event (e.g. the Manning’s n coefficient). These types of perturbations can result from large-scale erosion or sediment transport. Erosion and sediment transport can either occur over long time periods, on the order of months or years, between extensive data collection in the field, or over just several days, as a result of extreme events. Correctly estimating these types of perturbations to a parameter field may lead to more accurate forecasting of inundation for future extreme events. We focus on the estimation of the bottom friction characteristics as these have proven to be particularly difficult to determine from data and for which small perturbations can lead to dramatically different model forecasts. 3. Statistical data assimilation Solving a physical model numerically introduces errors (i.e. uncertainties) into model solutions that may already contain other significant sources of uncertainty. In the case of coastal ocean

45

models such as ADCIRC, major sources of uncertainty are unknown or poorly resolved model inputs, including boundary conditions, initial conditions, and model parameters such as the Manning’s n coefficient. This results in discrepancies between model output and observable data. Resolving such model/data discrepancies is a difficult task in itself, but this problem is further complicated by the fact that observable data are often contaminated by measurement error. We refer to all types of error and uncertainty as simply uncertainty in the model and data. These uncertainties can be quantified and reduced using statistical data assimilation (DA) methods. We present a brief overview of the Kalman and ensemble-based SEIK filter below and direct the interested reader to Tuan Pham et al. (1998), Hoteit et al. (2002), Evensen (2003), Evensen (2009a) and Kalman (1960) for a thorough pedagogical development of these filters as well as the renown ensemble Kalman filter. 3.1. General DA framework DA methods assume we are given an initial estimate of the model state, xak1 , at time tk1 , and the corresponding error covariance, Pak1 . The forecast model, Mk , is used to predict the state at some later time, tk :

xtk ¼ Mk xak1 þ Mk :

ð2Þ

While the model error, Mk , is assumed unknown, its associated error covariance, Q k , is given. The forecasted model state is defined as xfk :¼ Mk xak1 to differentiate it from the true, error-free state, xtk . In consideration of the unknown model error, Mk , statistical DA (filtering) methods attempt to determine the best estimate of the model state at time tk by utilizing data observed both prior to and at that time. We define the data obtained at time t k as yok . The current model state, xtk , is projected into the observation space using an observation operator, Hk , so that it can be directly compared to the current (and possibly noisy) data, yok :

ytk ¼ Hk xtk þ Hk :

ð3Þ

Here, ytk denotes the true noise-free data of the physical system. Again, it is assumed that the error in the observation operator, Hk , is unknown, but that the associated error covariance, Rk , is available. Often Rk is set to r2 I, where I is the identity matrix (with dimensions defined by those of the observation space) and r2 thus represents the error variance. This models the case of statistical errors for each datum that are independent and identically distributed, which is a standard assumption. We use the observation operator to compute the residual between the forecasted model state, xfk , and observed data, yok , which is then weighted by a matrix Gk . This weighted residual is used to update the forecasted model state. We define the updated state as the analysis, xak :

  xak ¼ xfk þ Gk yok  Hk xfk :

ð4Þ

Optimally determining xak defines the DA problem. 3.2. Ensemble Kalman filtering For a linear model, Mk , the DA problem can be solved by the Kalman filter (Kalman, 1960; Evensen, 2009a). In this method, the weight (known as the Kalman gain matrix), Gk , is computed in such a way that the error covariance of the analyzed model state, Pak , is minimized. First the initial error covariance matrix is evolved by the forecast model:

Pfk ¼ Mk Pak1 MTk þ Q k :

ð5Þ

46

T. Mayo et al. / Ocean Modelling 76 (2014) 43–58

This forecast error covariance matrix is then used to compute Gk : 1

Gk ¼ Pfk HTk ½Hk Pfk HTk þ Rk  : The updated error covariance, Pak , can be expressed as

Pak ¼ ½I  Gk Hk Pfk :

ð6Þ

In the case of a nonlinear model, Mk , a common approach to solving the DA problem is to use the ensemble Kalman filter (EnKF) (Evensen, 2003). Rather than forecasting the analyzed error covariance using the forecast model as in (5), which would require linearization of the model (as in the extended Kalman filter (Evensen, 2009a)), the error statistics of the system are represented by an ensemble of forecasted model states, xfk;i ; the EnKF uses a sampling based approach akin to Monte Carlo for propagating uncertainties forward in time. The initial ensemble xak1;i is determined by adding Gaussian noise to the initial state vector, xak1 , and the forecast operator is applied to each of the ensemble members to create the xfk;i . The error covariances in the Kalman filter equations are computed using the sample covariance of the ensemble members:

  T Pfk ¼ xfk;i  xfk xfk;i  xfk ;   T Pak ¼ xak;i  xak xak;i  xak : An ensemble of analyzed states, xak;i , results from updating the ensemble of forecasted states using yok;i , observations that have been perturbed from yok by the addition of Gaussian noise. This method has become increasingly popular in recent years due to its ease of implementation for computationally complex models for which the linearization required in methods such as the extended Kalman filter would be difficult. 3.3. The SEIK filter for state estimation The singular evolutive interpolated Kalman (SEIK) filter is a variation of the EnKF (Hoteit et al., 2002). It is from the class of squareroot ensemble Kalman filters (Nerger et al., 2012), which update the mean of the ensemble of forecasted states, xfk , along with a square root matrix of the respective error covariance (Tippett et al., 2003). This update of the mean eliminates the addition of Gaussian noise to the observations required in the EnKF. The SEIK filter is generally implemented in three steps: (1) the sampling step, (2) the forecast step, and (3) the analysis step. The sampling step uses second-order exact sampling so that the sample mean and sample covariance computed from the ensemble are exactly the analyzed state and error covariance, xak1 and Pak1 . The forecast step uses the model to evolve the ensemble members forward in time. The analysis step updates the forecasted mean and covariance. See Appendix C for details on computation and implementation specific to the SEIK filter. 3.4. The SEIK filter for parameter estimation Though traditionally used for state estimation, Kalman filters can also be used to estimate model parameters. Recently this has been done using joint estimation. In this approach, the state vector is augmented with uncertain model parameters. Then, both the state and parameters are updated by observed data through the Kalman filter equations (see e.g. Anderson (2001), Aksoy et al. (2006) and Evensen (2009b)). In this work, we use statistical data assimilation to estimate the model parameters alone. This can be done by redefining the state operators, Mk and Hk . It is assumed that the evolution of the model parameters, wak1 , is a static process, thus the identity operator is used as the parameter forecast model. Additive random noise can be included in the

forecast step to represent process noise, though this addition is an artificial inflation of the forecast error covariance and we do not include it in the numerical results shown in this paper, i.e Mwk ¼ 0, and wfk ¼ I wak1 ¼ wtk :

Mw k ¼ I; wtk ¼ I wak1 þ Mwk : A composition of the numerical forecast model and the observation operator used to forecast and observe the estimate of the initial model state xak1 in (2) and (3), respectively, is used to define a map from parameter space to observation space. The numerical forecast model uses the initial estimate of the model state along with the initial estimate of the unknown parameters to determine the model state at time tk . The state observation operator is then applied to this state vector. Since the model and consequently the state vector depend on the parameters, the observations of the state are implicitly a function of the parameters. Thus, the vector of observations is considered to be an observation of the initial parameters.

Hw k ¼ Hk  Mk ; t w ytk ¼ Hw k wk þ Hk :

The superscript w in the equations above denotes the operator defined for the model parameters (as opposed to those defined for the model state). The residual between the observed parameters and the observed data is weighted and used to update the forecasted parameters, wfk . These steps are then repeated. The updated parameters, wak , are first forecasted and then they are used to initialize the ADCIRC model along with the previous model state, projecting the forecasted parameters into the state space. The observation operator is applied to this projection to create an observation that is implicitly a function of the parameters. Note that the state vector is never updated, only the parameter vector. By redefining the operators in this way, the Kalman filters can be implemented as described above. 4. Numerical simulation and parameter estimation framework In this section we describe the basic framework used to apply the SEIK filter for parameter estimation to the ADCIRC model, including the general goal of numerical experiments for parameter estimation, how synthetic data is generated for observation system simulation experiments (OSSEs), and a description of the parameterizations of the Manning’s n coefficient in the domains considered. In Butler et al. (2012) and Altaf et al. (2013), we implemented statistical DA methodologies for state estimation with the goal of improving the accuracy of ADCIRC model output. We demonstrated an improvement in forecasts of water elevations using ADCIRC forced by hurricane wind data. Here, we aim to show improvement in both ADCIRC model input (i.e. parameters) and output by using statistical DA for parameter estimation. The primary goal is to obtain more accurate forecasts of water elevations by using the SEIK filter to estimate the bottom stress parameters, specifically the Manning’s n coefficients. We simulate data by specifying a ‘‘true’’ field of Manning’s n coefficients. This field is used along with tidal forcing in ADCIRC to generate ‘‘true’’ observations of water elevations. We then attempt to recover the true field from an inaccurate initial guess using the SEIK filter for parameter estimation. 4.1. A 1-D parameterized Manning’s n field We first aim to estimate a Manning’s n field parametrically defined by one value on a small computational domain. We let the

47

T. Mayo et al. / Ocean Modelling 76 (2014) 43–58

1000 500

0.3

nodes elements stations station 1 open ocean

0.2 0.1 0 elevation (m)

1500

0

−0.1 −0.2 −0.3

−500

−0.4

−1000 −0.5

−1500 −3000

−0.6

−2000

−1000

0

1000

2000

0

3000

10

20

30

40

50 60 time (h)

70

80

90

100

Fig. 2. Idealized inlet with ebb shoal (left) and water elevations modeled at station 1 for various Manning’s n coefficients, with larger tidal amplitudes corresponding to lower Manning’s n coefficients (right). (For interpretation of the references to colour in this figure caption, the reader is referred to the web version of this article.)

Manning’s n coefficient equal this value at every node of the grid to represent a field that is constant in space. 4.1.1. Classifying Manning’s n coefficients We seek to classify Manning’s n coefficients with respect to the response of modeled water elevation data. To begin, we divide Manning’s n coefficients, ranging from 0.005 to 0.2 and incremented by 0.005, into five classes. The classes are defined by first running the ADCIRC model with the smallest Manning’s n coefficient, 0.005, as smaller friction coefficients generate larger water elevations (i.e. tides with larger amplitudes). ADCIRC is forced using the principal lunar semi-diurnal (M2) tidal constituent on the open ocean boundary of the domains. The largest mean amplitude of the tides generated at several locations throughout the domain is computed, and the Manning’s n coefficients that generate tides with mean amplitudes less than 20%, 40%, 60%, 80%, and 100% of this value are then divided into five classes (i.e. subintervals) to identify Manning’s n coefficients that generate similar water elevation data. 4.1.2. Generating synthetic water elevation data We generate synthetic water elevation data by running the ADCIRC model with the constant Manning’s n coefficient set to a value from the center of each of the predefined classes. We force the ADCIRC model for 5 days, including a 12 h ramp up period, using the same M2 tidal constituent used in defining the parameter classes. The data generated from this model run is considered truth. 4.1.3. Recovering the true class of Manning’s n coefficients We attempt to recover the true class of Manning’s n coefficients by applying the SEIK filter for parameter estimation to the ADCIRC model. We begin with an initial guess of the true parameter from each of the four remaining classes. After the 12 h ramp up period, we assimilate the synthetic water elevation data every hour over the 5 days that data is available. We only assimilate the data at several locations throughout the domain, which we refer to as observation stations. The data assimilated does not contain measurement error, however, we discuss the effects of including measurement error in Appendix D. 4.2. A 2-D parameterized Manning’s n field We also aim to estimate two model parameters, a and b, which we use to parametrically define a 2-D field of Manning’s n coefficients. We denote this field by na;b . The exact type of parameterization will depend on the computational domain, and is described in the numerical results section below for each of our test cases. We

Table 2 Manning’s n coefficients in idealized inlet with ebb shoal classified by amplitude of corresponding tidal data at station 1. Class

Mean of tidal amplitude (m)

Manning’s n

A B C D E

0.267989–0.214391 0.214391–0.160793 0.160793–0.107195 0.107195–0.053598 0.053598–0

0.005–0.025 0.03–0.04 0.045–0.07 0.075–0.135 0.14–0.2

perform OSSEs in a similar manner as above; however, in this case we cannot create classes as we did when estimating one coefficient without significant computation. Instead we seek to obtain Manning’s n coefficients that produce maximum tidal amplitudes that are within 20% of those amplitudes that result when ADCIRC is run using the true bottom stress parameters. It should be noted that in all of our test cases we deliberately use the same model to generate synthetic data that we use to perform data assimilation experiments. This is to isolate the source of error in the water elevations; the parameter estimates are not overcompensating for other sources of model error. Additionally, the coastal ocean model has been verified and validated extensively. It gives reliable data provided the model parameters are accurate and the spatial domain is adequately resolved. Regarding the spatial resolution, refining the model domain used to generate synthetic data introduces only negligible changes to the numerical results presented below. For all of the experiments, we also compute the corresponding baseline cases. The baseline cases provide the results that would be

Table 3 x- and y-coordinates and bathymetry at the 15 observation stations in the idealized inlet. Sta.

x

y

h (m)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

1000 1000 1000 2000 2000 500 0 2000 1500 1500 1500 0 0 1000 1000

0 500 500 0 1000 0 0 1000 700 0 700 500 500 1000 1000

1.000 1.000 1.000 1.000 1.000 1.000 0.995 1.000 1.000 1.000 1.000 1.004 1.004 2.389 2.389

48

T. Mayo et al. / Ocean Modelling 76 (2014) 43–58

obtained under ordinary circumstances where a coastal ocean model is run with uncertain bottom stress parameters. This allows the effects of these incorrect parameters on the estimate of the model state to be directly measured. Furthermore, because our DA method is nonintrusive, we see the potential for improving estimated model states with only modest increases in computational resources. Implementation of the DA method requires a minimal number of additional forward model runs (one per ensemble member), and the matrix operations used to complete the analysis step of the filter. By comparing the results of the baseline cases to that of the cases where data is assimilated, we are able to quantify the improvement in the state estimation that the DA provides. 5. Numerical results and discussion We perform test cases on two computational domains, an idealized inlet with an ebb shoal, and a coarse representation of Galveston Bay. The first domain consists of a bay that is connected to the open ocean on the west via an inlet with twin jetties (Fig. 2). Immediately to the west of the inlet is an ebb shoal. The formation of an ebb shoal is a common occurrence at coastal inlets. An ebb shoal forms due to the deceleration of water as it exits the inlet, which results in the deposition of large amounts of sediment that are being transported by the flow (Kubatko et al., 2006). The second domain is Galveston Bay, a large estuary located on the upper Texas coast in the Gulf of Mexico. It is a densely populated, large metropolitan area, and it is also home to the Port of Houston, one of the world’s busiest ports, making it vital to the US economy. The location of Galveston Bay makes it particularly susceptible to hurricane storm surge. The Great Galveston Hurricane destroyed the region in 1990, and in 2008, Hurricane Ike devastated the area again in what has since been estimated as the third costliest hurricane in US history (Berg, 2008). 5.1. Inlet with ebb shoal The idealized inlet contains 1,518 nodes and 2,828 elements (Fig. 2). The west side of the inlet has a bathymetry that slopes linearly from 3.8 m to 1 m measured positive downward from the

geoid, and the remaining, landlocked portion on the east side has a constant bathymetry of 1 m. At the mouth of the channel there is an ebb shoal which is 750 m in diameter and has a maximum height of 0.4 m in the center. We run ADCIRC on this grid using a 2 s time step. ADCIRC is forced using the M2 tidal constituent with an amplitude of 0.25 m relative to the geoid.

5.1.1. Estimating a 1-D parameterized Manning’s n field The mean amplitude of the tides generated by running the ADCIRC model with a Manning’s n coefficient of 0.005 is 0.267989 m at station 1. We use this value along with the mean amplitudes generated by the other coefficients (Fig. 2)) to determine the five classes of Manning’s n coefficients (Table 2). We choose a Manning’s n coefficient from the approximate center of each of the five classes to generate synthetic water elevation data at 15 observation stations throughout the inlet (Table 3), denoted by the red circles in Fig. 2. (We discuss the effect of increasing the number of observation stations in Appendix E.) In five separate experiments, we let the true Manning’s n coefficient equal 0.015, 0.035, 0.06, 0.105, and 0.17. We then use the SEIK filter to estimate the true Manning’s n coefficient from initial guesses equal to each of the four remaining values. For the SEIK filter, we use an ensemble size of 10 and a value of r ¼ 0:01 (recall r2 is the coefficient in the observation error covariance matrix, R; our choice of small r signifies the confidence we have in our data, i.e. correct model parameters will result in the data we observe, however, we find that variation in the value chosen for this filter parameter does not have a significant effect on the results presented here). As previously mentioned, in addition to the data assimilation experiments, we run a baseline case, where no data is assimilated, to quantify the impact of the SEIK filter. In general, parameter estimation problems are approached with an attempt to exactly recover model parameters that produce model output matching observed data. However, for complex computational models, these problems are often ill-posed, with various parameter values resulting in very similar model output. Thus, the aforementioned approach is not necessarily the most effective. Instead, our efforts may be better spent attempting to recover values that produce the desired model output. In this experiement, we

Table 4 Final estimate errors of various true Manning’s n fields parameterized by two values for idealized inlet with ebb shoal domain. The RMSEs are of forecasted water elevations for the parameter estimation (SEIK filter) and baseline cases. True Manning’s n 0.015

Initial guess Final estimate Relative error RMSE (SEIK) RMSE (baseline)

0.035 0.014948 0.003460 0.014790 0.076792

0.06 0.014840 0.010640 0.022188 0.115541

0.105 0.014622 0.025180 0.011590 0.140996

0.17 0.014353 0.043134 0.017258 0.157038

True Manning’s n 0.035

Initial guess Final estimate Relative error RMSE (SEIK) RMSE (baseline)

0.015 0.034406 0.016967 0.018030 0.076792

0.06 0.034739 0.007469 0.014338 0.042259

0.105 0.034930 0.002011 0.017707 0.072842

0.17 0.034684 0.009015 0.018068 0.094232

True Manning’s n 0.06

Initial guess Final estimate Relative error RMSE (SEIK) RMSE (baseline)

0.015 0.052506 0.124892 0.035843 0.115541

0.035 0.054884 0.085252 0.023266 0.042259

0.105 0.060664 0.011067 0.013747 0.032302

0.17 0.060719 0.011977 0.013204 0.056430

True Manning’s n 0.105

Initial guess Final estimate Relative error RMSE (SEIK) RMSE (baseline)

0.015 0.058388 0.443927 0.058573 0.140996

0.035 0.090433 0.138733 0.018113 0.072843

0.06 0.087582 0.165878 0.024664 0.032302

0.17 0.107812 0.026781 0.007318 0.025872

True Manning’s n 0.17

Initial guess Final estimate Relative error RMSE (SEIK) RMSE (baseline)

0.015 0.061195 0.640027 0.076949 0.157038

0.035 0.136759 0.195536 0.028322 0.094232

0.06 0.094165 0.446086 0.046545 0.056430

0.105 0.151548 0.108542 0.014223 0.025872

49

T. Mayo et al. / Ocean Modelling 76 (2014) 43–58

seek to recover true Manning’s n coefficients with the goal of improving water elevation forecasts. We find that the modeled water elevations are generally not responsive to small perturbations in the value used for the Manning’s n coefficient in some classes, i.e. the forecast of water elevations are sensitive to changes in Manning’s n values between classes but the sensitivity within a class can often be small. Thus, rather than attempting to estimate the coefficient exactly, it suffices to recover a value that produces similar water elevation data. Here, we aim to obtain an estimate that lies within the same class as that of the true value (Table 2). Using the SEIK filter for parameter estimation, we find that we are able to accurately estimate the parameterized Manning’s n field defined by the three smallest classes using any of the four initial guesses, i.e. the final estimate lies in the same class as that of the true value (Table 4). In estimating the Manning’s n parameterization defined by 0.105 (i.e. class D), final estimates belonging to class D were obtained from all but the smallest of the initial guesses. In estimating the Manning’s n parameterization defined by 0.17 (i.e. class E), an estimate belonging to class E was only

0.18

0.16

0.035 0.06 0.105 0.17 truth class A

0.14 0.12 0.1 0.08 0.06 0.04

0.12 0.1 0.08 0.06 0.04 0.02

0.02 0

0 0

5

10

15

20 25 assimilation #

30

0.18

35

40

0.12

Manning’s n estimate

0.14

0

5

10

15

20 25 assimilation #

30

0.2

0.015 0.035 0.105 0.17 truth class C

0.16 Manning’s n estimate

0.015 0.06 0.105 0.17 truth class B

0.14 Manning’s n estimate

0.16 Manning’s n estimate

obtained from the largest initial guess. Because we incremented the Manning’s n coefficients by 0.005 when defining the classes, there are Manning’s n coefficients that do not lie within a class, i.e. they lie between two classes (Table 2). In these cases, we consider the parameter estimation a success if the final estimate of the Manning’s n coefficient lies in one of the intervals adjacent to the true class. With this consideration, a correct final estimate of 0.17 was also obtained from an initial guess from class B. In both of these cases (the estimates of 0.105 and 0.17), ‘‘correct’’ final estimates were recovered from the initial guesses that were from classes closest to that of the true coefficient. In general, the correct coefficient is underestimated, particularly in the case when the initial guess is less than the true value. Also, the relative error in the estimate of the Manning’s n coefficient generally increases with the true value. We discuss this trend below. In addition to estimating the model parameters, we also seek to improve forecasts of the model state. We compute the global root mean square error (RMSE) in the water elevations across all nodes at every assimilation time. As a result, this error is influenced by

0.1 0.08 0.06 0.04

35

40

0.015 0.035 0.06 0.17 truth class D

0.15

0.1

0.05

0.02 0

0 0

5

10

15

20 25 assimilation #

Manning’s n estimate

0.2

30

35

40

0

5

10

15

20 25 assimilation #

30

35

40

0.015 0.035 0.06 0.105 truth class E

0.15

0.1

0.05

0

0

5

10

15

20 25 assimilation #

30

35

40

Fig. 3. Estimates of various true Manning’s n fields parameterized by a single value for the idealized inlet with ebb shoal case. From top to bottom, left to right the true Manning’s n value is 0.015 (class A), 0.035 (class B), 0.06 (class C), 0.105 (class D), 0.17 (class E).

50

T. Mayo et al. / Ocean Modelling 76 (2014) 43–58 0.18

0.2

0.035 0.06 0.105 0.17

0.16 0.14

0.015 0.035 0.06 0.105

0.15

RMSE

RMSE

0.12 0.1 0.08

0.1

0.06 0.05

0.04 0.02 0

0 0

20

40

60 assimilation #

80

100

120

0

20

40

60 assimilation #

80

100

Fig. 4. The time evolution of the RMSE of the forecasted water elevations computed globally over all vertices of the idealized inlet domain for the estimation of 0.015 (left) and 0.17 (right).

the (inaccurate) initial guess of the Manning’s n coefficient (Fig. 4), but we find that it is significantly improved by the SEIK filtering over the baseline case in every instance. We note that the largest RMSE in the SEIK filtering case is 0.076949, and by comparison the baseline case has a corresponding RMSE that is over 100% larger (Table 4). In general, we see that the size of the classes of Manning’s n coefficients (Table 2) increases with the values, indicating that the water elevation data is not as sensitive to small variations in the larger Manning’s n coefficients. This idea is supported by noting that the baseline RMSEs in the water elevations generally decrease with increasing values of the true Manning’s n coefficients. The result is that smaller model/data discrepancies, and thus smaller subsequent filter updates, occur when larger

Manning’s n coefficients are estimated from classes with smaller values. As a result, we expect to experience more difficulty in recovering larger coefficients. For the case of 0.105, the filter only failed to attain an estimate in the desired class (class D) when the initial guess was from the class furthest from that of the true value (class A). We see similar behavior in the case of 0.17. Thus, for larger coefficients, the filter might simply require a more accurate initial guess or more time to attain the correct value. This idea is reflected in Fig. 3, as the variation in the estimates of each Manning’s n coefficient after 40 assimilation cycles increases with the true value; the parameter estimate requires more time to ‘‘converge’’ for larger values. We note that in a statistical DA scheme used for parameter estimation, the estimates never truly converge but rather ‘‘settle’’ near a value where the magnitude of the

Table 5 Final estimate errors of various true Manning’s n fields parameterized by two values for idealized inlet with ebb shoal domain. The RMSEs are of forecasted water elevations for the parameter estimation (SEIK filter) and baseline cases. Initial guess

n0:001;0:001

n0:005;0:005

n0:06;0:06

n0:1;0:1

n0:1;0:005

Estimate of a ¼ 0:005 Estimate of b ¼ 0:1

0.005092 0.085802

0.006633 0.087819

0.005004 0.098679

0.005721 0.100251

0.005325 0.097747

‘2 Error of n0:005;0:1 (SEIK) ‘2 Error of n0:005;0:1 (baseline) Relative error (a) Relative error (b)

0.014198 0.099081 0.018590 0.141977

0.012290 0.095000 0.326543 0.121809

0.001321 0.068007 0.000852 0.013206

0.000764 0.095000 0.144290 0.002512

0.002276 0.134350 0.065053 0.022526

RMSE (SEIK) RMSE (baseline)

0.028814 0.170614

0.027044 0.156683

0.004613 0.068212

0.010573 0.082921

0.008532 0.106257

Mean tidal amplitude (sta. 1) Relative error

0.160958 0.287481

0.155605 0.196570

0.128850 0.029745

0.115380 0.083531

0.132544 0.056787

0.001, 0.001 0.005, 0.005 0.1, 0.1 0.06, 0.06 0.1, 0.005 truth (alpha)

Manning’s n estimate

0.06 0.04 0.02 0

0.2

Manning’s n estimate

0.08

0.001, 0.001 0.005, 0.005 0.1, 0.1 0.06, 0.06 0.1, 0.005 truth (beta)

0.15

0.1

0.05

−0.02 −0.04

0 0

20

40

60 assimilation #

80

100

0

20

40

60 assimilation #

80

100

Fig. 5. Estimates of a (left) and b (right), parameters that define a true 2-D parameterization of a Manning’s n field defined on idealized inlet with ebb shoal.

51

T. Mayo et al. / Ocean Modelling 76 (2014) 43–58 0.3

0.3

0.2

0.2

0.1

0.1 0 elevation (m)

elevation (m)

0 −0.1 −0.2

−0.1 −0.2

−0.3

−0.3

−0.4

−0.4

−0.5

−0.5

0

10

20

30

40

50 60 time (h)

70

80

90

100

0

10

20

30

40

50 60 time (h)

70

80

90

100

Fig. 6. Water elevations at station 1 that result when one of the two parameters defining the Manning’s n field for the idealized inlet domain is fixed while the other is varied (fixed a and varying b (left), fixed b and varying a (right)).

6

3.3

x 10

0.15 0.1

3.28 0.05 elevation (m)

3.26 3.24 nodes elements stations station 1 station 9 station 21 open ocean

3.22 3.2 3.18 2.8

2.9

3

3.1

3.2

3.3

3.4

3.5

3.6

0 −0.05 −0.1 −0.15 −0.2

3.7

0

10

20

30

40

5

x 10

50 60 time (h)

70

80

90

100

Fig. 7. Galveston Bay (left) and water elevations modeled at station 1 for various Manning’s n coefficients, with larger tidal amplitudes corresponding to lower Manning’s n coefficients (right). (For interpretation of the references to colour in this figure caption, the reader is referred to the web version of this article.)

updates can be orders of magnitude smaller compared to initial updates to the parameter estimate. It is in this sense that we use the term converge to loosely describe this behavior of repeatedly small updates relative to the magnitude of initial updates. 5.1.2. Estimating a 2-D parameterized Manning’s n field For the second set of OSSEs, we define a field of Manning’s n coefficients over the domain (Fig. 2) through the two parameters, a and b. We let the Manning’s n coefficient equal a in the deep ocean area of the domain (west side, shown in green), and b in the shallow landlocked area (east side, shown in blue). We let the value of the coefficient increase linearly from a to b within the inlet (purple). We generate the synthetic water elevation data letting a = 0.005 and b = 0.1, i.e. using a ‘‘true’’ Manning’s n field n0:005;0:1 . We then use the SEIK filter to estimate the parameters that define this field from the various initial guesses n0:001;0:001 ; n0:005;0:005 ; n0:1;0:1 ; n0:06;0:06 , and n0:1;0:005 . In addition, we run the

Table 6 Manning’s n coefficients in Galveston Bay classified by amplitude of corresponding tidal data. Class

Mean of tidal amplitude (m)

Manning’s n values

A B C D E

0.099569–0.079655 0.079655–0.059742 0.059742–0.039828 0.039828–0.019914 0.019914–0

0.005–0.01 0.015–0.02 0.025–0.035 0.04–0.075 0.08–0.2

respective baseline cases (i.e. no data is assimilated) to quantify the impact of the SEIK filter. We are able to accurately estimate the true pair of Manning’s n coefficients from each of the five initial guesses (Table 5). In each of the five cases, the ‘2 (i.e., the euclidean norm) error of the final

Table 7 x- and y-coordinates and bathymetry at the 21 observation stations in Galveston Bay. Sta.

x (105 m)

y (106 m)

h (m)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21

3.362 3.308 3.231 3.116 3.199 3.079 3.174 3.205 3.322 3.250 3.237 3.311 3.229 3.545 3.169 3.077 3.249 3.053 2.901 2.922 2.938

3.294 3.281 3.277 3.273 3.267 3.284 3.264 3.262 3.263 3.249 3.249 3.247 3.244 3.265 3.242 3.237 3.241 3.233 3.225 3.218 3.219

1.573 2.524 2.793 2.797 1.946 3.547 1.827 8.086 2.183 7.446 4.876 6.388 10.670 1.875 3.439 1.878 1.405 0.895 0.845 0.963 4.239

52

T. Mayo et al. / Ocean Modelling 76 (2014) 43–58

estimate is less than 0.015. This is a significant reduction (at least 91%) in the ‘2 error of the initial guess, which is as high as 0.134. The final estimates of each distinct Manning’s n coefficient, a and b, are also accurate with respect to the classes defined in the first set of OSSEs (Table 2). We find that smaller ‘2 errors in the initial guess do not necessarily result in smaller ‘2 errors in the final estimate. This is also true for the estimation of a ¼ 0:005. On the other hand, we do see a relationship between the error in the initial guess of the second parameter, b ¼ 0:1, and the resulting error in the final estimate; a more accurate initial guess produces a more accurate final estimate. In all five cases, 0.005 is slightly overestimated and 0.1 is slightly underestimated with the exception of the initial guess n0:1;0:1 (see See Fig. 5). With regard to estimating the model state, the RMSEs of the water elevations obtained using the SEIK filter for parameter estimation decreases by at least 82% compared to the respective baseline cases. The water elevation error also decreases with decreasing ‘2 error in the parameters; we see interdependence between the model input and output. For four of the five cases, we find we are able to accurately forecast the mean of the tidal 0.14

0.14

0.02 0.03 0.06 0.14 truth class A

0.1 0.08 0.06 0.04 0.02

0.1 0.08 0.06 0.04 0.02

0

0 0

20

40

60 assimilation #

80

100

0.14

0

40

60 assimilation #

80

0.08 0.06 0.04

100

0.01 0.02 0.03 0.14 truth class D

0.12 Manning’s n estimate

0.1

20

0.14

0.01 0.02 0.06 0.14 truth class C

0.12 Manning’s n estimate

0.01 0.03 0.06 0.14 truth class B

0.12 Manning’s n estimate

0.12 Manning’s n estimate

amplitude (0.125018 m). In each of these cases, the relative error is less than 20%, as desired. The exception is the case with the initial guess n0:001;0:001 . However, this case is largely effected by the initial elevation errors. The estimated tidal amplitude at the end of this simulation is 0.141806 m, which is within 20% of the truth. In estimating b ¼ 0:1, we see the same behavior observed in the first set of OSSEs; higher Manning’s n coefficients are more difficult to reproduce and are generally underestimated. The results also give us some intuition about the sensitivity of water elevation data computed from the ADCIRC model. The accuracy of the initial guess of the first parameter, a ¼ 0:005, does not directly relate to the accuracy of the final estimate. Furthermore, if we hold the Manning’s n coefficient on the west side of the inlet constant with a ¼ 0:005 and vary b, the coefficient on the east side, we see far more variation in the resulting tidal amplitudes than we do if we hold the coefficient on the east side constant at b ¼ 0:1 and vary the coefficient a on the west side. The variation in the amplitudes in the former case is 0.197915 m whereas the variation in the latter case is only 0.078774 m (Fig. 6). Based on these results, we presume that the water elevation data is less sensitive to the value of the Manning’s n coefficient on the west side of the domain, a.

0.1 0.08 0.06 0.04

0.02

0.02 0

20

40

60 assimilation #

80

Manning’s n estimate

0.2

100

0

20

40

60 assimilation #

80

100

0.01 0.02 0.03 0.06 truth class E

0.15

0.1

0.05

0 0

20

40

60 assimilation #

80

100

Fig. 8. Estimates of various true Manning’s n fields parameterized by a single value for Galveston Bay case. From top to bottom, left to right the true Manning’s n value is 0.01 (class A), 0.02 (class B), 0.03 (class C), 0.06 (class D), and 0.14 (class E).

53

T. Mayo et al. / Ocean Modelling 76 (2014) 43–58

Table 8 Final estimate errors of various true Manning’s n fields parameterized by a single value for Galveston Bay case. The RMSEs are of forecasted water elevations for parameter estimation (SEIK filter) and baseline cases. True Manning’s n 0.01

Initial guess Final estimate Relative error RMSE (SEIK) RMSE (baseline)

0.02 0.009776 0.022383 0.003502 0.016351

0.03 0.009903 0.009720 0.004191 0.024598

0.06 0.009921 0.007877 0.007290 0.033514

0.14 0.009479 0.052067 0.015166 0.038380

True Manning’s n 0.02

Initial guess Final estimate Relative error RMSE (SEIK) RMSE (baseline)

0.01 0.020069 0.003440 0.002411 0.016351

0.03 0.019979 0.001044 0.001935 0.008966

0.06 0.019948 0.002582 0.003781 0.019577

0.14 0.019101 0.044947 0.009337 0.026017

True Manning’s n 0.03

Initial guess Final estimate Relative error RMSE (SEIK) RMSE (baseline)

0.01 0.029357 0.021433 0.003870 0.024598

0.02 0.030051 0.001706 0.001264 0.008966

0.06 0.030207 0.006885 0.002745 0.011161

0.14 0.030323 0.010777 0.009045 0.018304

True Manning’s n 0.06

Initial guess Final estimate Relative error RMSE (SEIK) RMSE (baseline)

0.01 0.046579 0.223687 0.007690 0.033514

0.02 0.055006 0.083235 0.004010 0.019577

0.03 0.053981 0.100310 0.002911 0.011161

0.14 0.069481 0.158009 0.003770 0.007729

True Manning’s n 0.14

Initial guess Final estimate Relative error RMSE (SEIK) RMSE (baseline)

0.01 0.051604 0.631399 0.013696 0.038380

0.02 0.072078 0.485155 0.008704 0.026017

0.03 0.087391 0.375777 0.006591 0.018304

0.06 0.112387 0.197233 0.002847 0.007729

This is reasonable since, although this side of the domain contains more nodes (1,073), it also has deeper bathymetry, reducing the value of the drag law coefficient, cf in (1). 5.2. Galveston Bay We proceed with the goal of successfully estimating a field of Manning’s n coefficients on a more realistic domain, Galveston Bay. We use a grid containing 2,113 nodes and 3,397 elements (Fig. 7). There is an open ocean boundary on the southeast side of the domain, and land boundaries surrounding the remainder. There are also island boundaries surrounding 17 islands located within the bay. The bathymetry of the grid ranges from 0.354 m inside the bay to 17.244 m in the Gulf of Mexico. We run ADCIRC on this grid using a 2 s time step, forced by the M2 tidal constituent with an amplitude of 0.1 m measured relative to the geoid. 5.2.1. Estimating a 1-D parameterized Manning’s n field The mean amplitude of the tides generated by running the ADCIRC model with a Manning’s n coefficient of 0.005 in Galveston Bay is 0.099569 m at station 1. We use this value to determine the five classes of Manning’s n coefficients (Table 6). 0.06

We choose a Manning’s n coefficient from the approximate center of each of the five classes to generate synthetic water elevation data at 21 observation stations throughout Galveston Bay (Table 7), denoted by the red circles in Fig. 7. In five separate experiments, we let the true Manning’s n coefficient equal 0.01, 0.02, 0.03, 0.06, and 0.14, respectively. We then use the SEIK filter to estimate the true Manning’s n coefficient from initial guesses equal to each of the four remaining values. We again use an ensemble size of 10 and a value of r ¼ 0:01. We aim to obtain an estimate that lies within the same class as that of the true value (Table 6). As we did when estimating the Manning’s n coefficient for the idealized inlet case, we run baseline cases, where no data is assimilated, to quantify the impact of the SEIK filter. We again find that we are able to accurately estimate nearly all of the Manning’s n coefficients from any initial guess (Table 8). The only exception occurs in the estimation of 0.14. In this case, final estimates belonging to the correct class are obtained from all but the smallest of the initial guesses, 0.01 and 0.02; accurate estimates are recovered from the two larger initial guesses, which come from the classes that are closest to the true coefficient (classes C and D). In addition to estimating the model parameters, we again aim to improve forecasts of the model state. We compute the global 0.045

0.02 0.03 0.06 0.14

0.05

0.035 0.03 RMSE

0.04 RMSE

0.01 0.02 0.03 0.06

0.04

0.03

0.025 0.02

0.02

0.015 0.01

0.01 0.005 0

0 0

20

40

60 assimilation #

80

100

0

20

40

60

80

100

assimilation #

Fig. 9. The time evolution of the RMSE of the forecasted water elevations computed globally over all vertices of the Galveston domain for the estimation of 0.01 (left) and 0.14 (right).

T. Mayo et al. / Ocean Modelling 76 (2014) 43–58 0.2

0.2

0.15

0.15

0.1

0.1

0.05

0.05

elevation (m)

elevation (m)

54

0 −0.05

0 −0.05

−0.1

−0.1

−0.15

−0.15

−0.2

−0.2 0

10

20

30

40

50 60 time (h)

70

80

90

100

0

10

20

30

40

50 60 time (h)

70

80

90

100

Fig. 10. Water elevations that result at station 1 when one of the two parameters defining the Manning’s n field for the Galveston Bay domain is fixed while the other is varied (fixed a and varying b (left), fixed b varying a (right)).

5.2.2. Estimating a 2-D parameterized Manning’s n field To estimate two parameters over Galveston Bay (Fig. 7), we again define a field of Manning’s n coefficients using the values a and b. We set the Manning’s n coefficient to a in the Gulf of Mexico (green), and b inside the bay (blue). We generate synthetic water elevations data letting a ¼ 0:005 and b ¼ 0:1. We denote this field by n0:005;0:1 . We observe that variation in the value of a does not cause notable variation in the modeled water elevations; the water elevation data is not sensitive to variations in the value of the Manning’s n coefficient in the deep water of the Gulf of Mexico (Fig. 10). Variation in the value of b, on the other hand, significantly impacts the

0.25

0.01 0.05 0.15 0.2 truth

0.2 Manning’s n estimate

RMSE in the water elevations across all nodes at every assimilation time and find that it is significantly improved by the SEIK filtering over the baseline case in every instance. The largest RMSE is 0.015166. This is the error in the case when the smallest Manning’s n coefficient (0.01) is estimated from the largest initial guess (0.14). Consequently this RMSE is inflated from the initial error in the water elevations (see Fig. 9). Nevertheless, the error is reduced by over 60% over the baseline case, and it is small in comparison to the mean tidal amplitudes generated using Manning’s n coefficients in this class. In estimating a constant Manning’s n coefficient in Galveston Bay, we observe the same trends seen when estimating the coefficient in the idealized inlet. Again, we see that the size of the classes of Manning’s n coefficients (Table 6) increases with the values, supporting the idea that the water elevation data is not as sensitive to small variations in the larger Manning’s n coefficients in this domain. We expect that for larger values of the Manning’s n coefficients, the SEIK filter will make smaller updates to the model parameters and thus require more time to converge to final estimates that lie within the desired class. This is observed in Fig. 8, as the variation in the final estimates of each Manning’s n coefficient increases with the true value, noticeably for 0.03, 0.06, and 0.14. The small RMSEs in the computed water elevations also support the idea that the updates made to the parameters are small. In the two cases where the Manning’s n coefficients are not accurately recovered (i.e. when 0.14 is estimated from 0.01 and 0.02), the RMSEs are less than 0.015. Because this error is computed globally over space and time, the elevation errors at the later assimilation times are smaller than those at earlier times. This results in negligible updates to the estimates of the Manning’s n coefficients after the initial assimilation cycles. Relative to the idealized inlet case, here the parameters require a longer time to converge to the true value. For the estimates of 0.06 and 0.14, there is still a notable amount of variation in the estimates of the coefficients even through the fifth day of assimilation.

0.15

0.1

0.05

0 0

100

200

300 400 assimilation #

500

600

700

Fig. 11. Estimates of b, a parameter that partially defines a true Manning’s n field defined on Galveston Bay.

modeled water elevations. For this reason, we choose to begin the data assimilation with the assumption that a ¼ 0:005 is known, and use the SEIK filter to estimate the parameter b only. We assimilate data at the 20 observation stations located in Galveston Bay, discarding data collected from station 17, the station located in the Gulf of Mexico. Including a 12 h ramp up period, we perform the data assimilation over 30 days, beginning with the initial guesses n0:005;0:01 ; n0:005;0:05 ; n0:005;0:15 , and n0:005;0:2 (Fig. 11). We are able to accurately estimate b from each of the initial guesses with the exception of the smallest initial guess, b ¼ 0:01 (Table 9). However, we find that the relative error of the final estimate of b is significantly reduced from that of the initial guess in each of the four experiments. In the worst case, the error is still reduced by more than 54%. We do not see a direct relationship between the accuracy of the initial guesses and that of the final estimates. Good final estimates can be obtained even when the initial guesses are as much as 100% inaccurate. The parameter is underestimated for two of the cases and overestimated for the remaining two. We see a direct relationship between this and the initial guesses, i.e. the underestimates result from the initial guesses that are less than 0.1 and the overestimates result from the initial guesses that are greater than 0.1. We see the impact of the SEIK filter when estimating the model state. The RMSEs of the estimated water elevations are reduced in every case respective to the baseline cases. Overall, the RMSEs of the estimated elevations are small, with the largest RMSE in the four experiments being less than 0.01. We again see the interdependence between ADCIRC model input and output, as the accuracy of the estimated Manning’s n coefficients closely relates to

55

T. Mayo et al. / Ocean Modelling 76 (2014) 43–58

Table 9 Final estimate errors of various true Manning’s n fields parameterized by two values. The RMSEs are of forecasted water elevations for the parameter estimation (SEIK filter) and baseline cases. Initial guess

n0:005;0:01

n0:005;0:05

n0:005;0:15

n0:005;0:2

Manning’s n estimate of b ¼ 0:1

0.059391

0.093012

0.104188

0.103496

Relative error (b) RMSE (SEIK) RMSE (baseline)

0.406090 0.008030 0.036703

0.069876 0.001151 0.007887

0.041880 0.000881 0.003102

0.034966 0.001054 0.004869

Mean tidal amplitude (sta. 1) Relative error (Sta. 9) Relative error (Sta. 21) Relative error

0.030060 0.916570 0.023167 0.339812 0.100270 0.006732

0.017511 0.116444 0.018212 0.053279 0.100877 0.000711

0.014246 0.091687 0.016467 0.047695 0.100999 0.000493

0.014294 0.088623 0.016466 0.047717 0.100995 0.000445

the accuracy of the estimated elevations. We obtain the mean tidal amplitude at stations 1, 9, and 21 (Fig. 7). For the true field of Manning’s n coefficients, the amplitudes are equal to 0.015684, 0.017291, and 0.100950, respectively. The mean tidal amplitudes that result when estimating the Manning’s n coefficient in the bay are very close to these values for most of the cases. The relative error is less than 0.12 in every experiment except when the initial guess is b ¼ 0:01. In this case the tidal amplitude is overestimated at stations 1 and 9. At these locations, the mean tidal amplitude is greatly impacted by the initial overestimate of this value that results from the underestimation of the bottom stress. In these experiments, significantly more assimilation cycles are required to accurately estimate the Manning’s n coefficient in the bay. With respect to the classes determined in the first set of OSSEs (Table 6), acceptable estimates are nearly attained after 108 updates, but more time is required for convergence to the true value as conjectured above. Class E, the class containing the true value of b ¼ 0:1, is the largest; Manning’s n coefficients ranging from 0.08 to 0.2 all produce similar elevation output. This implies that the incremental updates made to the parameter by the data assimilation method are small, and the longer convergence times are expected as a result. By the thirtieth day, final estimates within the correct class result from each of the initial guesses except for b ¼ 0:01, the initial guess that comes from class A, the class furthest from the true class. This idea of small incremental updates is also supported by the baseline water elevation data (Table 9). The RMSEs in the elevations for the baseline cases are small even though the Manning’s n coefficients in the bay are as much as 100% inaccurate. Finally, with this experiment we again find that the water elevation data is less sensitive to a, the value of the Manning’s n coefficient in deep water, as was the implication from the idealized inlet case.

6. Conclusions Through the implementation of the SEIK filter for parameter estimation to the ADCIRC model, we are able to accurately estimate Manning’s n coefficients in coastal lagoons and estuaries. We perform the parameter estimation using two spatial domains, an idealized inlet with an ebb shoal, as well as a coarse discretization of Galveston Bay. For both cases, we are able to accurately estimate one parameter, a Manning’s n coefficient held constant across the domain. We are also able to estimate two parameters used to define a field of Manning’s n coefficients in the idealized inlet. When estimating two parameters in Galveston Bay, we find that the parameter corresponding to the Manning’s n coefficient in shallow water areas of the domain can be accurately estimated. However, the water elevation data proves insufficient for estimating the parameter corresponding to the coefficient in the deeper water.

This is due to lack of sensitivity of the data to perturbations in this parameter. We expand on these findings below and suggest several areas for future research. In general, we are able to accurately recover smaller Manning’s n coefficients. Larger Manning’s n coefficients, however, require either a relatively accurate initial guess or many data assimilation cycles (updates by data) to converge to the true value. This is likely because our model output, water elevation data, is not as sensitive to variations in larger values of the Manning’s n coefficient; the number of coefficients that produce similar water elevation data increases with the value of the coefficient. As a result, the residuals between the measured and modeled data used in the Kalman filter equations are small and the updates made to the model forecasts are small as well. Larger values are generally underestimated, particularly when the initial guess is much less than the true value. In all of our experiments, we are able to improve the estimation of the model state. The global space–time RMSE of the water elevations is improved over the baseline case in every instance. In estimating two model parameters, where the Manning’s n coefficient is set to a small value in the deep water of the domain and a larger value in the shallow water, we find that the bathymetry is a significant factor in the sensitivity of the modeled water elevations to the bottom stress parameters. For areas with large bathymetric depths, the value of the Manning’s n coefficient has less of an effect on the corresponding water elevation. Parameter estimation efforts should thus be focused on areas of shallow water if water elevation data is to be used for the assimilation. The results of this study show that statistical data assimilation is a promising approach to parameter estimation. By using statistical data assimilation methods, we are able to efficiently estimate model parameters using only 10 ensemble members, i.e. 10 model runs. Furthermore, these methods eliminate the need for an adjoint model, and provide estimates of the uncertainty in the parameters. There are several areas for future research. These include exploring methods of determining the optimal number of observation stations, as well as overcoming the difficulties encountered when the observations are distorted by measurement error. Future research areas will also include more complex configurations of bottom stress parameters. Additionally, we will seek to estimate parameters for models where the water elevation data is more sensitive to the Manning’s n coefficient, e.g. sediment transport models and hydrodynamic models that allow for variable land boundaries (i.e. wetting and drying). We are also exploring the use of other types of observation data entirely. Finally, the method is well-suited for parallel implementation as it operates by projecting ensembles of parameters into the state space. This translates to running multiple instances of the forward model, for which parallel implementation would be ideal, particularly in the case when real-time state and parameter estimation is desired. For example, we are greatly interested in estimating parameters in

56

T. Mayo et al. / Ocean Modelling 76 (2014) 43–58

hydrodynamic models forced by hurricane wind data, which would allow for more accurate real-time forecasting of hurricane storm surges. These and many other research questions can be informed by statistical data assimilation methods for parameter estimation.

Appendix C. The SEIK filter framework Here we briefly describe details of notation and implementation of the three stages specific to the low-rank SEIK filter used in this paper.

Acknowledgments C.1. Sampling step This work was supported by the King Abdullah University of Science and Technology and the Gulf of Mexico Research Initiative Center for Advanced Research on Transport of Hydrocarbons in the Environment. This support is gratefully acknowledged. Appendix A. The shallow water equations

Pak1 ¼ Lk1 Uk1 LTk1 ;

Depth integration of the incompressible Navier–Stokes equations, assuming hydrostatic pressure, produce the shallow water equations (SWEs). The SWEs appropriately model fluid flow of processes with large horizontal length scales relative to the vertical length scales. The result of the depth integration is a first-order hyperbolic continuity equation for water elevation coupled to the momentum equations for horizontal depth-averaged velocities given by

@H @ @ þ ðQ Þ þ ðQ y Þ ¼ 0; @t @x x @y

ðA:1Þ

@Q x @UQ x @VQ x @½f þ P s =g q0  ag ssx þ þ  fQ y ¼ gH þ @x @t @x @y q0 

sbx þ Mx  Dx  Bx ; q0

In the sampling stage, an ensemble of size r þ 1 is drawn from the initial n-dimensional analyzed state, xak1 . Here, r  n is the rank of the associated error covariance matrix, Pak1 . Pak1 is factored into the matrices Lk1 and Uk1 ,

where Lk1 is n  r and Uk1 is r  r. A Cholesky decomposition is used to factor Uk1 , T U1 k1 ¼ Ck1 Ck1 ;

which allows Pak1 to be written as 1 T T Pak1 ¼ Lk1 CT k1 Xk1 Xk1 Ck1 L k1 :

Here, Xk1 is a random ðr þ 1Þ  r matrix with orthonormal columns and zero column sums. This expression allows an ensemble of analyzed states to be drawn in such a way that its mean and sample covariance are xak1 and Pak1 , respectively:

xai;k1 ¼ xak1 þ ðA:2Þ

@Q y @UQ y @VQ y @½f þ Ps =g q0  ag ssy sby þ þ  fQ x ¼ gH þ  @y @t @x @y q0 q0 þ M y  Dy  By :

ðC:1Þ

pffiffiffiffiffiffiffiffiffiffiffi T r þ 1Lk1 ðXk1;i C1 k1 Þ :

ðC:2Þ

Xk1;i is the ith row of Xk1 . Its inclusion in (C.2) introduces stochastic variation to the sampled states. C.2. Forecast step In the forecast stage, each of the analyzed states xai;k1 in (C.2) is forecasted using the numerical model until time tk , the time at which observed data is available:

Here, f is the free surface departure from the geoid, h is the bottom surface departure from the geoid (bathymetric depth), H ¼ f þ h is the total water column height, U i is the depth-averaged velocity in the xi direction, Q xi ¼ U xi H is the flux per unit width in the xi direction, f is the Coriolis parameter, Ps is the atmospheric pressure at the free surface, q0 is the reference density of water, a is the effective earth elasticity factor, g is the Newtonian equilibrium tide potential, ssxi are the imposed surface stresses, sbxi are the bottom stress components, M xi is the vertically-integrated lateral stress gradient, Dxi is the momentum dispersion, and Bxi is the vertically-integrated baroclinic pressure gradient.

Here, T is an ðr þ 1Þ  r full rank matrix with zero column sums. Its inclusion allows the equality in (C.3). Pf can then be written as

Appendix B. Derivation of the GWCE for ADCIRC

Pfk ¼ Lk Uk1 LTk ;

Here, we describe the derivation of the second-order, hyperbolic generalized wave continuity equation (GWCE), which replaces the continuity equation (A.1) in the SWEs. First, a multiple, s0 P 0, of the continuity equation is added to the time derivative of the continuity equation. Then, it is assumed @H ¼ @f , i.e. the bathymetric depth is constant. Substituting the @t @t momentum equations into this equation completes the derivation:

 @2f @f @ @UQ x @VQ x @½f þ Ps =g q0  ag þ s þ   þ fQ y  gH 0 @t @x @x @x @y @t 2   @UQ y @VQ y ssx sbx @  þ M x  Dx  Bx þ s0 Q x þ   þ @y q0 q0 @x @y  @½f þ Ps =g q0  ag ssy sby þ  þ M y  Dy  By þ s0 Q y þfQ x  gH @y q0 q0 @ s0 @ s0  UH  VH ¼ 0: ðB:1Þ @x @y

xfi;k ¼ Mk xai;k1 : The forecasted state, xfk , is taken as the mean of this ensemble. To forecast the error covariance, Pak1 ; Lk is computed using the ensemble of forecasted states, xfi;k :

h i h i Lk ¼ xf1;k  xfk . . . xfrþ1;k  xfk T ¼ xf1;k . . . xfrþ1;k T:

ðC:3Þ

where 1

Uk1 ¼ ½ðr þ 1ÞTT T : C.3. Analysis step The forecasted state, xfk , is updated using the observed data, yok , to form an analyzed state, xak . First Uk is updated: T 1 1 U1 k ¼ Uk1 þ ðHLÞk R ðHLÞk ;

where

h i ðHLÞk ¼ Hk xf1;k . . . Hk xfrþ1;k T: It is then used to compute the Kalman gain matrix:

Gk ¼ Lk Uk ðHLÞTk R1 :

57

T. Mayo et al. / Ocean Modelling 76 (2014) 43–58

The analyzed state is

Appendix E. The effect of an increased number of observation stations

  xak ¼ xfk þ Gk yok  Hk xfk and the updated error covariance is

Pak ¼ Lk Uk LTk :

Appendix D. The effect of observation error In the experiments presented in Section 5, the data assimilated into the forward model contains no observation error. The purpose of this paper is to present and investigate a new method of parameter estimation, and thus the initial test cases were simplified so that initial questions concerning the effectiveness of the method could be answered. However, in practice, it can be expected that data will always be subject to measurement error. Here we discuss the impacts of adding random noise to the observations. The results for the estimation of the smallest and largest classes are shown in Table D.1. When estimating 0.015, the observation error seems to have little effect on the ability of the filter to recover the accurate parameter; the correct class is attained from all of the initial guesses as is the case when the observations contain no error. The relative errors are slightly larger for all cases other than that when the initial guess 0.17, however the relative errors are still small (less than 0.033). In the estimation of the model state, the maximum RMSE is 0.026, and is reduced by over 90% compared to the baseline in three cases. For the fourth case, the reduction is 67%. None of the RMSEs are significantly different than those for the cases with no observation errors (Table 4) as they are within approximately 1 cm. When estimating a constant Manning’s n coefficient from the largest class, 0.17, we are unable to accurately recover the parameter. The correct class is only obtained from one of the initial estimates, 0.035, and the relative error is as much as 0.765. The RMSEs of the estimated states are less than 9 cm, which is an improvement over the baseline cases for all initial guesses other than 0.105. The RMSEs are comparable to the cases with noise-free observations again with the exception of the initial guess 0.105. We find that while we are able to recover the smaller class of Manning’s n coefficients when the observations contain random noise, it is very difficult to get good estimates of the parameter from the larger class. The initial guess does not seem to impact the final estimate in this case. By introducing random noise to the observations, the experiments contain a large amount of variability, i.e. multiple simulations can produce vastly different results. Additionally, it is likely that the specification of the observation error covariance plays a more significant role in these experiments. As data will often contain some measurement error, this is an area for future research.

In practice, observation stations are chosen based on the locations where data is truly available. The locations chosen for the results presented in Section 5, particularly those for the Galveston Bay case, are representative of what could be expected in practice, and we are able to obtain good parameter estimates with these limited number of observation stations. However, more observations stations will become available with technological advances, e.g. as sensor networks are deployed. Here we discuss the impact of increasing the number of observation stations in the idealized inlet domain. We estimated Manning’s n coefficients from the smallest and largest classes, increasing the number of observation stations from 15 to 24 (Fig. E.1). The results of the experiments are shown in Table E.2. The final estimates of the constant Manning’s n coefficient 0.015 are very accurate. A value from the correct class is recovered from every initial guess. Compared to the experiments using fewer observation stations, the parameters are not underestimated but slightly overestimated. We see improvement in the relative error of the final estimate for every initial guess as expected. The RMSEs are slightly worse except for the case when the initial guess is 0.17, but the difference is not significant and all of the RMSEs are within 1 cm of each other. For the estimates of 0.17, the correct class is only obtained from the initial guess 0.105. The relative error is improved over the cases with fewer observation stations for all initial guesses except 0.035. The RMSEs of the estimated states are also improved with the exception of this case. By increasing the number of observation stations in the domain, we see improvements in most of the parameter estimates. However, we expect that there is a limit on the additional information

1500 1000 500

nodes elements stations station 1 open ocean

0 −500 −1000 −1500 −3000

−2000

−1000

0

1000

2000

3000

Fig. E.1. Locations of 24 observation stations in idealized inlet with ebb shoal.

Table D.1 Final estimate errors of various true Manning’s n fields parameterized by a single value for idealized inlet with ebb shoal case. Random noise, not exceeding 10%, has been added to each of the observations. The RMSEs are of forecasted water elevations for the parameter estimation (SEIK filter) and baseline cases. True Manning’s n 0.015

initial guess Final estimate Relative error RMSE (SEIK) RMSE (baseline)

0.035 0.014820 0.011983 0.025556 0.076792

0.06 0.015492 0.032776 0.010444 0.115541

0.105 0.014581 0.027946 0.013910 0.140996

0.17 0.015235 0.015646 0.015104 0.157038

True Manning’s n 0.17

initial guess Final estimate Relative error RMSE (SEIK) RMSE (baseline)

0.015 0.040681 0.760699 0.083359 0.157038

0.035 0.174588 0.026986 0.022672 0.094232

0.06 0.080805 0.524678 0.045081 0.056430

0.105 0.040034 0.764505 0.086874 0.025872

58

T. Mayo et al. / Ocean Modelling 76 (2014) 43–58

Table E.2 Final estimate errors of various true Manning’s n fields parameterized by a single value when 24 observation stations are used for idealized inlet with ebb shoal case. The RMSEs are of forecasted water elevations for the parameter estimation (SEIK filter) and baseline cases. True Manning’s n 0.015

initial guess Final estimate Relative error RMSE (SEIK) RMSE (baseline)

0.035 0.015001 0.000067 0.014962 0.076792

0.06 0.015007 0.000489 0.024880 0.115541

0.105 0.015009 0.000572 0.019377 0.140996

0.17 0.015038 0.002549 0.016722 0.157038

True Manning’s n 0.17

initial guess Final estimate Relative error RMSE (SEIK) RMSE (baseline)

0.015 0.072624 0.572802 0.066872 0.157038

0.035 0.118255 0.304382 0.032531 0.094232

0.06 0.123587 0.273018 0.028701 0.056430

0.105 0.162138 0.046247 0.010227 0.025872

gained as the number of stations increases. At some point, the information added to the system through the filtering becomes redundant. Determining the optimal number of observation stations for a given domain of interest remains an active research area. References Aksoy, A., Zhang, F., Nielsen-Gammon, J., 2006. Ensemble-based simultaneous state and parameter estimation in a two-dimensional sea-breeze model. Mon. Weather Rev. 134 (10), 2951–2970. Altaf, M., Butler, T., Luo, X., Dawson, C., Mayo, T., Hoteit, I., 2013. Improving short range ensemble Kalman storm surge forecasting using robust adaptive inflation. Mon. Weather Rev. 141 (2013), 2705–2720. Anderson, J., 2001. An ensemble adjustment Kalman filter for data assimilation. Mon. Weather Rev. 129 (12), 2884–2903. Bacopoulos, P., Hagen, S.C., Cox, A.T., Dally, W.R., Bratos, S.M., 2012. Observation and simulation of winds and hydrodynamics in St. Johns and Nassau rivers. J. Hydrol. 420, 391–402. Berg, R., 2009. Tropical cyclone report: hurricane Ike (al092008), 1–14 september 2008, National Hurricane Center. Blanton, B., McGee, J., Fleming, J., Kaiser, C., Kaiser, H., Lander, H., Luettich, R., Dresback, K., Kolar, R., 2012. Urgent computing of storm surge for north carolina’s coast. Proc. Comput. Sci. 9, 1677–1686. Bunya, S., Dietrich, J., Westerink, J., Ebersole, B., Smith, J., Atkinson, J., Jensen, R., Resio, D., Luettich, R., Dawson, C., et al., 2010. A high-resolution coupled riverine flow, tide, wind, wind wave, and storm surge model for southern Louisiana and Mississippi. Part I: model development and validation. Mon. Weather Rev. 138 (2), 345–377. Butler, T., Altaf, M.U., Dawson, C., Hoteit, I., Luo, X., Mayo, T., 2012. Data assimilation within the advanced circulation (ADCIRC) modeling framework for hurricane storm surge forecasting. Mon. Weather Rev. 140, 2215–2231. Dawson, C., Proft, J., 2004. Coupled discontinuous and continuous Galerkin finite element methods for the depth-integrated shallow water equations. Comput. Methods Appl. Mech. Eng. 193 (3), 289–318. Dietrich, J., Bunya, S., Westerink, J., Ebersole, B., Smith, J., Atkinson, J., Jensen, R., Resio, D., Luettich, R., Dawson, C., et al., 2010. A high-resolution coupled riverine flow, tide, wind, wind wave, and storm surge model for southern Louisiana and Mississippi. part II: synoptic description and analysis of hurricanes Katrina and Rita. Mon. Weather Rev. 138 (2), 378–404. Dietrich, J., Westerink, J., Kennedy, A., Smith, J., Jensen, R., Zijlema, M., Holthuijsen, L., Dawson, C., Luettich Jr, R., Powell, M., et al., 2011. Hurricane Gustav (2008) waves and storm surge: hindcast, synoptic analysis, and validation in southern Louisiana. Mon. Weather Rev. 139 (8), 2488–2522. Dietrich, J., Trahan, C., Howard, M., Fleming, J., Weaver, R., Tanaka, S., Yu, L., Luettich, R., Dawson, C., Westerink, J., et al., 2012. Surface trajectories of oil transport along the northern coastline of the Gulf of Mexico. Continental Shelf Research 41, 17–47.

Evensen, G., 2003. The ensemble Kalman filter: theoretical formulation and practical implementation. Ocean Dyn. 53 (4), 343–367. Evensen, G., 2009a. Data Assimilation: The Ensemble Kalman Filter. Springer Verlag. Evensen, G., 2009b. The ensemble Kalman filter for combined state and parameter estimation. Control Syst., IEEE 29 (3), 83–104. Ghanem, R., Spanos, P., 2003. Stochastic Finite Elements: A Spectral Approach. Dover Pubns. Hagen, S.C., Bacopoulos, P., 2012. Coastal flooding in Florida’s big bend region with application to sea level rise based on synthetic storms analysis. Terr. Atmos. Oceanic Sci. 23 (5), 481–500. Hoteit, I., Pham, D., Blum, J., 2002. A simplified reduced order Kalman filtering and application to altimetric data assimilation in tropical pacific. J. Mar. syst. 36 (1), 101–127. Kalman, R.E. et al., 1960. A new approach to linear filtering and prediction problems. J. Basic Eng. 82 (1), 35–45. Kennedy, A.B., Gravois, U., Zachry, B.C., Westerink, J.J., Hope, M.E., Dietrich, J.C., Powell, M.D., Cox, A.T., Luettich, R.A., Dean, R.G., 2011. Origin of the hurricane Ike forerunner surge. Geophys. Res. Lett. 38 (8). Kinnmark, I., 1986. The Shallow Water Wave Equations: Formulation, Analysis, and Application, vol. 15. Springer-Verlag, Berlin, Heidelberg. Kubatko, E.J., Westerink, J.J., Dawson, C., 2006. Discontinuous Galerkin methods for advection dominated problems in shallow water flow. Comput. Methods Appl. Mech. Eng. 196 (1), 437–451. Luettich, R., Westerink, J., 2004. Formulation and numerical implementation of the 2D/3D ADCIRC finite element model version 44.XX. Lynch, D.R., Gray, W.G., 1979. A wave equation model for finite element tidal computations. Computers & fluids 7 (3), 207–228. Medeiros, S.C., Hagen, S.C., Weishampel, J.F., 2012. Comparison of floodplain surface roughness parameters derived from land cover data and field measurements. J. Hydrol. 452, 139–149. Nerger, L., Janjic, T., Schröter, J., Hiller, W., 2012. A unification of ensemble square root Kalman filters. Mon. Weather Rev. 140 (7), 2335–2345. Pham, D., 2001. Stochastic methods for sequential data assimilation in strongly nonlinear systems. Mon. Weather Rev. 129 (5), 1194–1207. Tippett, M.K., Anderson, J.L., Bishop, C.H., Hamill, T.M., Whitaker, J.S., 2003. Ensemble square root filters. Mon. Weather Rev. 131 (7), 1485–1490. Triantafyllou, G., Hoteit, I., Petihakis, G., 2003. A singular evolutive interpolated Kalman filter for efficient data assimilation in a 3-d complex physical– biogeochemical model of the cretan sea. J. Mar. Syst. 40, 213–231. Tuan Pham, D., Verron, J., Christine Roubaud, M., 1998. A singular evolutive extended Kalman filter for data assimilation in oceanography. J. Mar. syst. 16 (3), 323–340. Westerink, J.J., Luettich, R.A., Feyen, J.C., Atkinson, J.H., Dawson, C., Roberts, H.J., Powell, M.D., Dunion, J.P., Kubatko, E.J., Pourtaheri, H., 2008. A basin-to channelscale unstructured grid hurricane storm surge model applied to southern Louisiana. Mon. Weather Rev. 136 (3), 833–864. Wiener, N., 1938. The homogeneous chaos. Am. J. Math. 60 (4), 897–936. Yanagi, T., 1999. Coastal Oceanography, vol. 1. Springer.