A three-dimensional cohesive sediment transport ...

Estuarine, Coastal and Shelf Science xxx (2016) 1e14

Contents lists available at ScienceDirect

Estuarine, Coastal and Shelf Science journal homepage: www.elsevier.com/locate/ecss

A three-dimensional cohesive sediment transport model with data assimilation: Model development, sensitivity analysis and parameter estimation Daosheng Wang a, Anzhou Cao b, Jicai Zhang b, c, *, Daidu Fan d, Yongzhi Liu a, Yue Zhang d a

Key Laboratory of Physical Oceanography (Ocean University of China), Ministry of Education, Qingdao 266100, China Institute of Physical Oceanography, Ocean College, Zhejiang University, Zhoushan 316021, China State Key Laboratory of Satellite Ocean Environment Dynamics, Second Institute of Oceanography, Hangzhou 310012, China d State Key Laboratory of Marine Geology, Tongji University, Shanghai 200092, China b c

a r t i c l e i n f o

a b s t r a c t

Article history: Received 23 March 2016 Received in revised form 21 August 2016 Accepted 25 August 2016 Available online xxx

Based on the theory of inverse problems, a three-dimensional sigma-coordinate cohesive sediment transport model with the adjoint data assimilation is developed. In this model, the physical processes of cohesive sediment transport, including deposition, erosion and advection-diffusion, are parameterized by corresponding model parameters. These parameters are usually poorly known and have traditionally been assigned empirically. By assimilating observations into the model, the model parameters can be estimated using the adjoint method; meanwhile, the data misfit between model results and observations can be decreased. The model developed in this work contains numerous parameters; therefore, it is necessary to investigate the parameter sensitivity of the model, which is assessed by calculating a relative sensitivity function and the gradient of the cost function with respect to each parameter. The results of parameter sensitivity analysis indicate that the model is sensitive to the initial conditions, inflow open boundary conditions, suspended sediment settling velocity and resuspension rate, while the model is insensitive to horizontal and vertical diffusivity coefficients. A detailed explanation of the pattern of sensitivity analysis is also given. In ideal twin experiments, constant parameters are estimated by assimilating ‘pseudo’ observations. The results show that the sensitive parameters are estimated more easily than the insensitive parameters. The conclusions of this work can provide guidance for the practical applications of this model to simulate sediment transport in the study area. © 2016 Elsevier Ltd. All rights reserved.

Keywords: Cohesive sediment Transport model Data assimilation Adjoint method Sensitivity analysis Parameter estimation

1. Introduction Coastal waters are often characterized by the high concentration of suspended sediments derived from discharge of rivers or seabed resuspension (Miller and McKee, 2004). High suspended sediment concentrations (SSCs) have significant impact on phytoplankton productivity (Cloern, 1987; May et al., 2003), coral growth (McLaughlin et al., 2003), productivity of submerged aquatic vegetation (Dennison et al., 1993), the transport of pollutants (Martin and Windom, 1991), nutrient dynamics (Mayer et al., 1998) and geomorphic evolution (Jia et al., 2006; Fan et al., 2014). Consequently, a better understanding of suspended sediment

* Corresponding author. Current address: Ocean College, Zhejiang University, 1 Zheda Road, Zhoushan, Zhejiang Province, China. E-mail address: [email protected] (J. Zhang).

transport process is essential. The dynamics of suspended sediment transport is very complex in estuaries and coastal areas, and SSC varies as a function of tides, wind, circulation, etc. (Fettweis et al., 2007). In addition, the distribution of suspended sediments in coastal environments has large spatial and temporal variability (Miller and McKee, 2004; He et al., 2013). In-situ measurement is the most direct method to investigate the mechanism of suspended sediment transport; however, it only provides a local description (Amoudry and Souza, 2011). Therefore, the traditional field sampling method is inadequate to synoptically map suspended sediments. Over past decades, polar-orbiting satellite ocean color data (e.g., Myint and Walker, 2002; Warrick et al., 2004; Petus et al., 2010; Zhang et al., 2010) and geostationary meteorological satellite data (e.g., Neukermans et al., 2009, 2012; Salama and Shen, 2010) have been used to map suspended sediments in various coastal regions. However, only the

http://dx.doi.org/10.1016/j.ecss.2016.08.027 0272-7714/© 2016 Elsevier Ltd. All rights reserved.

Please cite this article in press as: Wang, D., et al., A three-dimensional cohesive sediment transport model with data assimilation: Model development, sensitivity analysis and parameter estimation, Estuarine, Coastal and Shelf Science (2016), http://dx.doi.org/10.1016/ j.ecss.2016.08.027

2

D. Wang et al. / Estuarine, Coastal and Shelf Science xxx (2016) 1e14

surface SSC can be obtained from the satellite data, while the vertical structure of SSC is still challenging our modelling and observing capabilities. Numerical suspended sediment transport models have gradually become a powerful tool to study suspended sediments (e.g., Hayter and Mehta, 1986; Wang and Pinardi, 2002; Guan et al., 2005; Leupi et al., 2008; Hu et al., 2009; Amoudry, 2014). Bagnold (1935, 1937) firstly attempted to correlate sediment transport with fluid dynamics. The mathematical modelling of sediment transport has become a major subject in coastal sciences since the work carried out by Einstein (1950) and higher development of the computing capacity. Broadly speaking, based on spatial scale, sediment transport models can be divided into onedimensional models (e.g., Thomas and Prasuhn, 1977; Krishnappan, 1981; Holly, 1990; Papanicolaou et al., 2004), twodimensional models (e.g., Onishi and Wise, 1982; Van Rijn and Tan, 1985; Jia and Wang, 1999) and three-dimensional models (e.g., Blumberg and Mellor, 1987; Hamrick, 1992; Song and Haidvogel, 1994). The details about the development stages of current representative models and their main applications, strengths and limitations are available in the review article by Papanicolaou et al. (2008). Amoudry and Souza (2011) reviewed state-of-the-art Eulerian implementations of bottom-up sediment transport and morphological change in coastal ocean hydrodynamic models. They noted that a very important characteristic of present sediment transport models is the high degree of empiricism, and another issue is the mismatch between model results and the experimental data. There are a large number of parameters in sediment transport models, including the initial conditions, inflow open boundary conditions, settling velocity, resuspension rate, and horizontal and vertical diffusivity coefficients. Unfortunately, it is very difficult to assign reasonable values to these poorly constrained parameters, which leads to large discrepancies between modelling results and observational data. This is the one of the largest limitations for the application of sediment transport models. As we know, in oceanic and atmospheric studies, data assimilation methods have been widely used to reduce the data misfit between models and observations (Zhang and Wang, 2014). For sediment transport models, data assimilation methods also provide possibilities to achieve this goal, which will be studied in this work. The four-dimensional variational (4DVAR) data assimilation is one of the most powerful and effective methods among all of the data assimilation approaches. The adjoint method has been successfully implemented in parameter estimation for different models (e.g., Navon et al., 1992; Zou et al., 1995; Alekseev et al., 2009; Zhang and Lu, 2010; Zhang et al., 2011; Cao et al., 2013). Navon (1998) presented a significant overview on the state of the art of parameter estimation in meteorology and oceanography with respect to applications of 4DVAR data assimilation techniques to inverse parameter estimation problems. In this paper, based on the theory of inverse problems, a threedimensional sigma-coordinate cohesive sediment transport model with the adjoint data assimilation is developed. In Section 2, the cohesive sediment transport model and the adjoint model are constructed, and the model settings are described. The parameter sensitivity analysis and the estimation of several constant parameters are performed in Sections 3 and 4, respectively. Finally, a summary is provided in Section 5 of this paper. 2. Models 2.1. The cohesive sediment transport model The governing equation of the cohesive sediment transport

model is a three-dimensional advection diffusion equation with vertical particle settling, which is written as:

vC vC vC vC v vC v vC þu þv þw KH þ KH ¼ vt vx vy Hvs vx vx vy vy v vC vC þ KV þ ws Hvs Hvs Hvs

(1)

where C represents the SSC; t is the time; x and y are the horizontal coordinates; s is the vertical coordinate; H is the total water depth, including the undisturbed water depth and sea surface elevation; u, v and w are the flow velocity components in the x, y and s directions; KH and KV are the horizontal and vertical diffusivity coefficients; and ws denotes the cohesive sediment settling velocity. At the bottom of the water column, the cohesive sediment flux is equal to the combination of deposition and erosion, which is described as:

ws $C

KV vC ¼ E D; H vs

s¼0

(2)

where E and D are the erosion rate and deposition rate, respectively. Erosion occurs when the bottom shear stress exceeds the threshold of erosion (Yang and Hamrick, 2003, hereinafter referred to as YH03):

E¼

M0 $ðtb =tce 1Þ; 0 ; tb tce

tb > tce

(3)

Conversely, deposition occurs when bottom shear stress is smaller than a critical value (YH03):

D¼

ws $C1 $ð1 tb =tcd Þ; 0; tb tcd

tb < tcd

(4)

where M0 is the resuspension rate; tb is the bottom shear stress; tce and tcd are the critical shear stress for erosion and deposition, respectively; and C1 is the SSC near the bottom. Generally, the settling velocity ws is a function of salinity and sediment concentration (Owen, 1970; Hayter, 1983; Dyer, 1986; Wang et al., 2013). For simplicity, constant settling velocity is applied in this study, without considering the effects of salinity and sediment concentration. According to Einstein (1950), erosion and deposition cannot occur simultaneously; therefore, the same critical shear stress for erosion and deposition are employed in this study. At the sea surface, the no flux boundary condition is applied, which is written as:

ws $C

KV vC ¼ 0; H vs

s¼1

(5)

In addition, the model is subject to the first boundary conditions at the inflow boundary U1, no gradient boundary conditions at the outflow boundary U2 and the solid boundary U3, which are given by:

CjU1 ¼ Cobc

(6)

! vC=v n jU2 ¼ 0

(7)

! vC=v n jU3 ¼ 0

(8)

! where Cobc denotes the SSC at the inflow boundary and n is the vector normal to the boundary. Furthermore, the model has initial conditions C0, which are described by:



Cjt¼0 ¼ C0

(9)

The equations are discretized with Arakawa A grid to solve the model numerically. The finite difference schemes of Eq. (1), Eq. (2) and Eq. (5) are displayed in Appendix A.1. 2.2. The adjoint model Both continuous and discrete adjoint models have been used in previous works. Continuous adjoint models are derived from continuous forward equations, which are then discretized, whereas discrete adjoint models are directly derived from the discretization of the forward model (Sirkes and Tziperman, 1997). Sandu et al. (2005) noted that discrete adjoint models are preferred for data assimilation problems because they provide the exact derivative of the discrete function being minimized. Therefore, the discrete adjoint model is adopted for this study. To construct the discrete adjoint model, the cost function is first defined as:

J¼

2 W X n ~n Ci;j;k C i;j;k 2

3

The first-order derivative of the Lagrangian function with respect to the variables and parameters should be zero to minimize the cost function J:

vL

¼0

(12)

vL n ¼0 vCi;j;k

(13)

vL ¼0 vpni;j;k

(14)

vlni;j;k

where p denotes the parameters in the cohesive sediment transport model, including KH, KV, ws, M0, Cobc and C0. In fact, Eq. (12) is the same as Eq. (A.1), which is the discrete form of the governing equation. From Eq. (13), the discrete adjoint model can be derived, which is displayed in Appendix A.2. From Eq. (14), gradients of the cost function with respect to the parameters can be obtained, which are shown in Appendix A.3.

(10)

i;j;k;n

n and C ~ n are the simulated and observed SSC at the (i, j, where Ci;j;k i;j;k k) grid point at the n time step, and W is the weighting matrix, which should theoretically be the inverse of the observational error covariance matrix. W can be simplified assuming that the errors of the data are uncorrelated and equally weighted (Yu and O'Brien, 1992). In the present study, W is 1 where observations are available, and 0 otherwise. Based on the Lagrangian multiplier method, the Lagrangian function is defined as:

2.3. Model settings 2.3.1. Study area and background flow field The numerical experiments in this work are performed to test the feasibility and ability of the developed model. The horizontal resolution of the study area of Hangzhou Bay (Fig. 1) is 0.01 0.01, and there are five uniform s-layers in the vertical direction. The hydrodynamic background field in the model domain is calculated by the Finite Volume Coastal Ocean Model (FVCOM,

2 nþ1 3 n nþ1 nþ1 n n n n Ci;j;k Ci;j;k wnþ1 Ci;j;kþ1 Ci;j;k1 Ciþ1;j;k Ci1;j;k Ci;jþ1;k Ci;j1;k 9 8 i;j;k n n þ ui;j;k $ þ vi;j;k $ þ nþ1 $ 6 7> > > 6 7> D t 2Dx 2Dy 2Ds H > > > 6 7> i;j > > > 6 7> > > > 6 7> > > ! > > 6 7 n n n n n n n n > > > 6 1 ðKH Þiþ1;j;k þ ðKH Þi;j;k Ciþ1;j;k Ci;j;k ðKH Þi;j;k þ ðKH Þi1;j;k Ci;j;k Ci1;j;k 7> > > > 6 7> > $ $ > > 6 Dx 7> > > 2 D x 2 D x > > 6 7 > > > > 6 7 > > > > 6 7 > > ! > > 6 7 > > n n n n n n n n > > 6 7 C C C C ðK Þ þ ðK Þ ðK Þ þ ðK Þ > > 1 H H H H P i;jþ1;k i;j;k i;jþ1;k i;j;k i;j;k i;j1;k i;j;k i;j1;k > > 6 7 n > > $ $ > > l $ 6 7 > > i;j;k > > D y 2 D y 2 D y 6 7 > > k > > 6 7 > > > > 6 7 > > > > 6 7 > > ! = nþ1 nþ1 nþ1 nþ1 nþ1 nþ1 nþ1 nþ1 6 7 X< ðKV Þi;j;kþ1 þ ðKV Þi;j;k Ci;j;kþ1 Ci;j;k ðKV Þi;j;k þ ðKV Þi;j;k1 Ci;j;k Ci;j;k1 7 1 6 L¼Jþ 6 7 $ $ > nþ1 nþ1 6 H nþ1 Ds 7> 2 2 i;j;n > Hi;j Ds Hi;j Ds > > 6 7> > > i;j > 6 7> > > > 6 7> > > > 6 7> > > > > nþ1 nþ1 nþ1 6 7 > > C C ðw Þ > > s i;j;k i;j;kþ1 i;j;k1 4 5 > > > > nþ1 $ > > > > 2 Ds > > H > > i;j > > > > > > > > " # > > nþ1 nþ1 nþ1 nþ1 nþ1 n n > > > > C þ C C C ðK Þ l þ l V i;j;1 > > i;j;0 i;j;0 i;j;1 i;j;0 i;j;1 nþ1 i;j;1 nþ1 nþ1 > > > > þ ðws Þi;j;1 $ þ $ D þ E > > i;j i;j nþ1 > > 2 2 Ds > > Hi;j > > > > > > > > > > n " # > > n nþ1 nþ1 nþ1 nþ1 nþ1 > > C þ C C C ðK Þ ; : li;j;K þ li;j;Kþ1 V i;j;K i;j;Kþ1 i;j;Kþ1 i;j;K nþ1 i;j;K ðws Þi;j;K $ þ $ þ nþ1 2 2 Ds Hi;j

where l denotes the adjoint variable of C; K is the total number of vertical layers, and there is one false layer at the bottom, numbered 0, while Kþ1 represents that at the surface.

(11)

Chen et al., 2003). The computational domain is configured with non-overlapped triangular mesh, with horizontal resolution varying from 0.4 to 0.6 km. Spherical and sigma coordinates are


4


30.8

11

30.7

10

30.6

9

30.5

8

30.4

7

30.3

6

30.2

5

30.1

4

30

3

29.9

120.6

120.8

121

121.2 Longitude(E)

121.4

121.6

Depth (m)

Latitude(N)

12

2

Fig. 1. Model domain and bathymetry of the regional Hangzhou Bay.

selected, in which the water column is divided into five equal vertical sigma layers. The open boundary is set at 121.8 E, where tidal boundary conditions including M2, S2, K1 and O1 are implemented. Temperature and salinity open boundary conditions are set as Blumberg and Kantha implicit radiation conditions. The Qiantang River is divided into three contiguous grids from 120.58 E, 30.36 N to 120.58 E, 30.40 N as freshwater sources. The river boundary forcing is specified by using the multi-year monthly averaged river runoff in July 1200 m3/s (Du et al., 2010). The model totally runs for 50 days. The simulated results in the initial 40 days are discarded to make sure that the simulation is stable, and the results of the final 10 days are used to interpolate the hydrodynamic field, which is then used to force the cohesive sediment transport model. 2.3.2. Initial and open boundary conditions The initial and open boundary conditions are critical for the present model. Cai et al. (2015) analyzed the influence of both seabed topography and tidal currents on SSC distribution using in-

situ SSC, remote sensing reflectance of sea water, water depth and simulated currents. They found that SSC has a negative correlation with bathymetry, which is especially obvious during mid to late period of flood tide. Surface initial conditions are shown in Fig. 2 according to the approximated relationship between SSC and water depth. The initial conditions at other layers are given by multiplying the surface initial conditions with coefficients (3.5e0.5k), where k is the number of the layer and the values of k are equal to 1 and 5 for the bottom layer and surface layer, respectively. Additionally, the inflow open boundary is set at 120.58 E and the outflow open boundary is set at 121.8 E. The sediment concentrations at the west inflow boundary are set to 1.5 kg/m3 according to Du et al. (2010). 2.3.3. Parameter settings For this study, the integral time step is 600 s, and the total simulation time is 120 h, which is enough to test the ability of the developed model. As stated by YH03, the settling velocity and resuspension rate

0.65 30.8

0.6 sediment concentration (kg/m3)

Latitude(N)

30.7 30.6

0.55

30.5

0.5

30.4

0.45

30.3 30.2

0.4

30.1

0.35

30 29.9

0.3 120.6

120.8

121

121.2 Longitude(E)

121.4

121.6

Fig. 2. The initial sediment concentration at the surface layer. Black dotted lines show locations of the three vertical cross sections used for parameter sensitivity analysis.



are two critical parameters controlling the sediment exchange process between the water column and the sediment bed, and their values are important in numerical modelling of the cohesive sediment transport. The settling velocity and resuspension rate of cohesive sediment vary over the range of 107 to 103 m/s and 103 to 102 g/m2/s, respectively, according to the literature (e.g., Krone, 1962; Gibbs, 1985; Mehta et al., 1989; Sheng et al., 1992; Sanford and Halka, 1993). These two parameters, especially the settling velocity, vary over a wide range, which significantly increases the difficulties of adjustment if data assimilation is not used. Agrawal and Pottsmith (2000) described two sensors for the measurement of particle size-distribution and settling-velocity distribution, and found that the settling velocity data from a field experiment off the New Jersey coast fit the model:

ws ¼ 0:45 103 a1:2 n

(15)

where an is the radius in microns and the settling velocity ws is in cm/s. The median grain size of bottom material in Yangtze Estuary and Hangzhou Bay is about 0.008 mm (Hu et al., 2009). The settling velocity is calculated to be about 5.0 105 m/s by Eq. (15), which agrees with the value in YH03. According to Hu et al. (2009), the resuspension rate is set to 1.0 106 kg/m2/s. As previously stated, the critical erosion stress and critical deposition stress are equal in this model. The critical stress is set to 0.1 N/m2, which is consistent with the range of critical stress in Hu et al. (2009). According to the dimensional analysis, the horizontal and vertical diffusivity coefficients are 80 and 1 103 m2/s, respectively. 3. Parameter sensitivity analysis Amoudry and Souza (2011) recommended that modelling studies should include a sensitivity analysis. There are a large number of parameters in the cohesive sediment transport model; therefore, it is necessary to investigate which parameters are sensitive before performing parameter estimation. A relative sensitivity function (Zou et al., 1993a; Cacuci et al., 2005) is used to quantify the effects of each parameter on the model results, which is written as follows:

SV;p ¼

Vpþ% Vp Vp ðpþ% pÞ=p

(16)

where SV,p is the relative sensitivity, p is the value of a model parameter, and pþ% indicates the value of the parameter increased by a percentage of p. V is the total SSC at the final time step in a certain area, and Vpþ% is the analogous value obtained with pþ%. It is obvious that the dimensionless quantity SV,p can serve as a guide to ranking the importance of the different model parameters. To evaluate the results of the sensitivity analysis comprehensively, the total SSC at three vertical cross sections at 120.8 E, 121.2 E and 121.6 E (as shown in Fig. 2), five s layers and the whole region are selected to explore how the parameters affect the simulation of SSC. To calculate sensitivity values, each parameter is increased by 25%. These experiments are classified as Group 1 and the sensitivity results (i.e., SV,p) are listed in Table 1. From Table 1, the resulting variations of SSC at all selected regions and the variations of C0 or M0 have a positive correlation. Because C0 is the initial condition, SSCs in all areas will increase when C0 becomes larger. M0 is the resuspension rate and the sediment flux from the bottom will be increased with larger M0, which will lead to higher SSCs in all areas. Additionally, the currents at approximately 120.8 E are stronger than other areas; therefore, the increased SSC at section 120.8 E is much larger. When Cobc is

5

Table 1 Values of sensitivity analysis to parameter variations (25% perturbation) in Group 1. Areas

120.8 E 121.2 E 121.6 E s¼1 s ¼ 0.8 s ¼ 0.6 s ¼ 0.4 s ¼ 0.2 Whole

Parameters C0

Cobc

KH

KV

M0

ws

6.31E-02 8.36E-01 8.41E-01 7.47E-01 7.49E-01 7.51E-01 7.52E-01 7.53E-01 7.51E-01

7.69E-01 7.66E-07 7.32E-07 1.01E-01 9.88E-02 9.68E-02 9.47E-02 9.26E-02 9.65E-02

1.82E-02 1.08E-02 2.17E-02 1.07E-02 1.08E-02 1.08E-02 1.08E-02 1.11E-02 1.08E-02

2.77E-02 3.27E-02 3.56E-02 1.18E-01 8.01E-02 4.27E-02 5.33E-03 3.10E-02 3.95E-02

1.68E-01 1.64E-01 1.59E-01 1.52E-01 1.52E-01 1.53E-01 1.53E-01 1.54E-01 1.53E-01

1.02Eþ00 8.14E-01 8.69E-01 8.78E-01 8.41E-01 8.03E-01 7.65E-01 7.27E-01 7.99E-01

increased, the SSCs in all regions are also increased, but the increases at the vertical cross sections at 121.2 E and 121.6 E are much smaller than others. It is likely that the suspended sediments from the inflow open boundary have limited influence on remote areas due to the transport capacity of the background flow. In contrast with C0, M0 and Cobc, the SSCs at all selected areas are negatively correlated with ws. A larger value of ws will enhance sedimentation, which then decreases the SSCs. The SSCs at all selected areas are negatively correlated with KH, except the vertical cross section at 121.6 E. When KH becomes larger, the effect of horizontal diffusion is enhanced, which discharges the suspended sediment at the outflow boundary condition. At the vertical cross section at 121.6 E, the increases are maybe affected by the local gradient of SSCs. When KV becomes larger, the SSCs exhibit different patterns depending on water depth, namely decreased values at the bottom and increased values in the upper part of the water column. However, the SSCs become larger at three vertical cross sections and in the whole region. These patterns occur because KV is the vertical diffusivity coefficient, which transports the suspended sediment at deep depths to shallower depths, and maintains the balance of SSC vertically with ws. On the whole, parameters can be divided into two groups. In the first group, including C0, Cobc, KV and M0, variations of parameters are positively correlated with SSC. On the contrary, KH and ws belong to the second group for which this correlation is negative. The values of relative sensitivity indicate that KH and KV are the least sensitive parameters of the two groups. For the adjoint method, the cost function described by Eq. (10) is the implicit function of the model parameters above. The gradients of cost function with respect to these parameters can also be used to do sensitivity analysis. According to Li and Wunsch (2004), if the impact of model parameters on the cost function were larger, that parameter was more confined to the observations, meaning that particular parameters can be estimated more accurately. Fan and Lv (2009) further suggested that the magnitudes of gradients of the cost function with respect to the parameters could be used to evaluate the sensitivity of different types of parameters. Based on this method, we made another group of experiments, classified as Group 2. All parameters are perturbed by 25%, and the average values at each grid point of the gradients of the cost function with respect to each parameter after one assimilation step are then calculated. The magnitudes of the results are displayed in Table 2. From Table 2, it is stated that the gradients of the cost function with

Table 2 The magnitudes of the mean value of the gradients of the cost function with respect to each parameter at every grid point in Group 2. Parameters

C0

Cobc

KH

KV

M0

ws

Magnitudes

100

101

107

101

101

101


6


respect to M0 and ws at every grid point are at least one magnitude higher than those of other parameters, indicating that the model results are sensitive to them. Conversely, the magnitude of the gradient of the cost function with respect to KH is at a minimum, at least 6 orders of magnitude lower than other parameters. Clearly, the results of Group 2 exhibit the similar tends as Group 1. Therefore, combining the results of Groups 1 and 2 reveals that C0, M0, Cobc and ws are sensitive parameters for the present model, while KH and KV are insensitive parameters. 4. Parameter estimation 4.1. Description of twin experiments To test the abilities of the adjoint model developed in this paper, the constant parameters ws, M0, Cobc, KH and KV are estimated in ideal twin experiments (TEs). In these experiments, hourly SSCs at the surface are assumed to be ‘pseudo’ observations which are obtained by running the forward model with the prescribed value of a specific parameter. To test the validity of data assimilation, onetenth of the total observations are randomly selected as independent observations, and the corresponding grid points are termed ‘checking points’ (CPs). At each CP, the observations are used to verify model results, rather than being assimilated. The remaining nine-tenths of the surface grid points are ‘assimilating points’ (APs), at which observations are assimilated into the model. Specifically, the difference in SSC between simulated values and ‘pseudo’ observations serves as the external force of the adjoint model. By assimilating the ‘pseudo’ observations at APs, the parameter is adjusted using the steepest descent method from an initial guess value, which is written as follows:

pkþ1 ¼ pk a$VJðpÞ

(17)

Where, p is the parameter which is adjusted, VJðpÞ is the gradient of cost function with respect to this parameter, and a is step length of the steepest descent method. By assimilating the ‘pseudo’ observations, the parameters will be optimized continuously, and the difference between simulated values and ‘pseudo’ observations can be reduced. In this methodology, the gradients of cost function with respect to a certain parameter calculated by the equations in Appendix A.3 are important because the direction of greatest decrease of the cost function is determined by the direction opposite to the gradient vector. The number of iteration steps of assimilation is an important variable for parameter estimation. Too many iteration steps will introduce waste of computing resources. The number of iteration steps depends on many factors, such as the characteristic of parameters, the magnitude of smoothness, and the initial guess value. The iteration of optimization will terminate once the following criterion is achieved:

kGk < eps maxð1; kXkÞ

(18)

where kGk is the L2 norm of the gradients of cost function with respect to control variables (i.e., model parameters), eps is a positive variable that determines the accuracy with which the solution is to be found, and kXk is the L2 norm of control variables. Both the values of kGk and kXk vary with the iteration steps. For a correct adjoint model, kGk will gradually decrease as the number of iteration steps is increased, and the inverted values of model parameters must gradually approach the prescribed ‘true values’. In this work, the assimilation process is terminated at the same iteration step in all numerical experiments to compare the results of

different parameters. Through trial and error, 20 iteration steps are selected, which satisfy Eq. (18) for all experiments. More details on the description of twin experiments can be found in Alekseev et al. (2009), and Zhang and Lu (2010). 4.2. Benchmark experiments For benchmark experiments (TE1), the prescribed constant values are 5 105 m/s, 1.0 106 kg/m2/s, 1.5 kg/m3, 80 m2/s and 1 103 m2/s for ws, M0, Cobc, KH and KV, respectively. In TE1, a certain parameter is estimated solely and the initial guess value is set to 0. The step lengths of the steepest descent method are 1.8 1015, 5 1016, 1 104, 1.275 1023 and 5.28 1025 for ws, M0, Cobc, KH and KV, respectively. The cost function that is normalized by its initial value (i.e., the value at the first iteration step) and the values of the estimated parameter versus the iteration steps are shown in Fig. 3. Form Fig. 3A, it is found that the normalized cost function tends to be stable after 14 iterations, and also reach its minimum. Correspondingly, the inverted ws nearly perfectly matches the prescribed value after 14 iterations. Fig. 3B indicates that the normalized cost function becomes stable after 7 iterations, while the inverted M0 converges to the prescribed value. The inverted Cobc is approximately equal to the prescribed value after 9 iterations in Fig. 3C. After two iteration steps, the normalized cost function is reduced by an order of 1020, while the inverted KH approximately converges to the given value (Fig. 3D). Similar results have been obtained for KV, as shown in Fig. 3E. In addition, the normalized cost function at the last iteration step, the mean absolute errors (MAEs) between simulated results and ‘observations’ at the APs and CPs, and the MAEs between the estimated values and the prescribed values are calculated, which are shown in Table 3. The resulting decrease of MAEs between simulated results and ‘observation’ at the CPs is almost equal to that of the APs in all experiments. In addition, Table 3 indicates that all parameters are estimated successfully, which is identical to the results obtained from Fig. 3. Thus, for simplicity, error statistics will not be shown for the latter experiments, and the inversion results will be shown by the figures. 4.3. Influence of step length As shown by the results of sensitivity analysis, KH and KV are insensitive parameters, which generally are difficult to be estimate. To identify this conjecture, two new groups of twin experiments (TE2 and TE3) are carried out, where the magnitude of smoothness (i.e., the step length of the steepest decent method) is changed. Compared with TE1, the step length in TE2 is decreased by 50% in the estimation of each parameter, and in TE3 the step length is increased by 50%. The results of TE2 and TE3 are exhibited in Figs. 4 and 5, respectively. In Fig. 4A, it is shown that the normalized cost function reduces at all of the iteration steps, and reaches the order of 108 in the last step. Correspondingly, the inverted ws gradually converges to the prescribed value. Fig. 4B and C indicate that the inverted M0 and Cobc at the last step are also approximately equal to the prescribed values. These results demonstrate that ws, M0 and Cobs which are sensitive parameters can be inverted successfully with different choices of step length. However, this conclusion is different for insensitive parameters. As shown by Fig. 4D, it is obvious that the normalized cost function is reduced in the first two iteration steps, and the inverted value of KH is far from the prescribed one, which is the same for the inversion of KV in Fig. 4E. Although the estimated ws is greater than the prescribed value in the second iteration step, the estimated value is adjusted to be close to the prescribed value in the latter steps (Fig. 5A). From



10 Iteration Steps

Cost function

(B)

Normalized cost function Inverted M0

−9

10

−18

10

0

5

10 Iteration Steps

Cost function

(C)

Normalized cost function Inverted Cobc

−9

10

5

10 Iteration Steps

−3

0

10

(E)

Normalized cost function Inverted KV

−9

10

10

KH (m2/s)

0 20

15

x 10 2

1

0

5

10 Iteration Steps

0 20

15

3

1.5

−18

10

0

160

80

−18

0 20

15

0

10

−11

10

x 10 2

1

Normalized cost function Inverted KH

10

−6

0

10

(D)

−22

0 20

15

Cost function

5

M0 (kg/m2/s)

0

0

10 Cost function

0.5

−20

10

ws (m/s)

−10

10

x 10 1

Cobc (kg/m3)

Cost function

(A)

Normalized cost function Inverted ws

KV (m2/s)

−4

0

10

7

0

5

10 Iteration Steps

0 20

15

Fig. 3. Results of TE1. (A) Read line marked by triangles shows the variations of the normalized cost function, and blue line with asterisks denotes variations of inverted settling velocity ws. (B) Same as (A) but for the resuspension rate M0. (C) Same as (A) but for the inflow open boundary condition Cobc. (D) Same as (A) but for the horizontal diffusivity coefficient KH. (E) Same as (A) but for the vertical diffusivity coefficient KV. Note that, in all panels, the black dotted horizontal lines delineate the prescribed values. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Table 3 Error statistics of TE1 before and after assimilation. Parameter

K1a

K2a

K3a

K4a

Before

After

Before

After

Before

After

ws M0 Cobc KH KV

6.81E-16 1.38E-13 1.22E-14 1.47E-20 9.75E-17

5.61E-01 2.91E-02 2.73E-02 1.03Eþ03 7.48Eþ02

1.28E-08 1.28E-08 1.28E-08 2.94E-06 1.17E-04

5.62E-01 2.90E-02 2.17E-02 8.86Eþ02 7.82Eþ02

1.27E-08 1.27E-08 1.27E-08 2.77E-06 1.17E-04

5.00E-05 1.00E-06 1.50Eþ00 8.00Eþ01 1.00E-03

1.10E-15 1.51E-16 2.16E-10 9.59E-03 1.63E-06

a K1 is the normalized cost function at the last iteration step; K2 and K3 are the MAEs between simulated results and ‘observations’ at the APs and the CPs, respectively, whose units are in kg/m3; and K4 is the MAEs between prescribed and inverted parameter values, whose units are the same as the units of the corresponding parameters.

Fig. 5B and C, it is shown that the normalized cost functions are stable after several iteration steps, meanwhile the inverted M0 and Cobc converge to their prescribed values. However, the estimated KH and KV cannot be adjusted to be close to their prescribed values (Fig. 5D and E). Overall, the results of TE2 and TE3 reveal that the sensitive parameters, including ws, M0 and Cobc, can be inverted successfully in any reasonable step length. However, the estimation of insensitive parameters, including KH and KV, will be easily affected by step length and can be inverted successfully only with a specific step length.

KH and KV are equal to the initial guesses during all iteration steps (Fig. 6D and E). Most probably, the reason is that KH and KV are insensitive which leads to a small difference between the ‘pseudo’ observations generated using the prescribed values and the simulated results using the initial guesses. The results of TE5, shown by Fig. 7, indicate that ws, M0 and Cobc can be inverted successfully, while the estimations of KH and KV failed again. Therefore, the results of TE4 and TE5 prove that the sensitive parameters can be estimated successfully with a wide range of initial guesses, while the choice of initial guess of insensitive parameters is strict. 4.5. Comparison with limited-memory BFGS algorithm1

4.4. Influence of initial guess The value of initial guess can also affect the estimation of parameters, which is discussed in TE4 and TE5. The initial guess in each experiment of TE4 and TE5 is set to 50% and 150% of the prescribed value, respectively, while the other settings remain the same as TE1. The results of TE4 are displayed in Fig. 6. From Fig. 6A, it is seen that the normalized cost function tends to be stable and reduces to the order of 1015 after 14 steps, while the inverted ws approaches the prescribed value. Fig. 6B and C demonstrate that M0 and Cobc can also be estimated successfully. However, the inverted

The Broyden-Fletcher-Goldfarb-Shanno (BFGS) method is an efficient quasi-Newton conjugate-gradient algorithm, which has been widely used in the unconstrained minimization of the cost function in the inverse problem (Zou et al., 1993b; Alekseev et al., 2009). The limited-memory BFGS (L-BFGS) algorithm is an adaptation of the BFGS method to large minimization problem (Zhang and Lu, 2010). Zou et al. (1993b) carried out a study with several

1

Suggested by an anonymous reviewer.



0.5

−10

0

5

10 Iteration Steps

15

Cost function

(B)

−6


−6

10

1

−12

10

0

5

10 Iteration Steps

15

5

10 Iteration Steps

15

0 20 −3

0

(E)


−6

10

x 10 2

1

−12

0

5

10 Iteration Steps

15

0 20

3


−6

10

1.5

−12

10

80

0

10

160

Cobc (kg/m3)

Cost function

(C)

−6

10

0 20

0

10


10

10

x 10 2

0

10

(D)

−12

0 20

Cost function

10

Cost function

ws (m/s)

−5

10

0

10

M0 (kg/m2/s)

Cost function

(A)


KH (m2/s)

−4

x 10 1

0

10

KV (m2/s)

8

0

5

10 Iteration Steps

15

0 20

Fig. 4. Same as Fig. 3, but for TE2 in which the step length is decreased by 50%.

−18

5

10 Iteration Steps

15

Cost function

(B)


−8

10

0

5

10 Iteration Steps

15

Cost function

(C)


−9

10

80 Normalized cost function Inverted KH

0

5

10 Iteration Steps

15

−3

(E)

−10

10

10

0

0 20 x 10 2

0

10

1 Normalized cost function Inverted KV

−20

0 20

0

10

5

10 Iteration Steps

15

0 20

3

1.5

−18

10

x 10 2

1

−16

10

10 −6

0

10

−11

10

−22

0 20

Cost function

0

M0 (kg/m2/s)

10

160

(D)

KH (m2/s)

0.5

0

10 Cost function

−9

10

x 10 1

ws (m/s)


Cobc (kg/m3)

Cost function

(A)

KV (m2/s)

−4

0

10

0

5

10 Iteration Steps

15

0 20

Fig. 5. Same as Fig. 3, but for TE3 in which the step length is increased by 50%.

limited-memory quasi-Newton methods, and concluded that the LBFGS had the best overall performance for the examined problems. In TE6, several experiments using the L-BFGS version of Liu and Nocedal (1989) are carried out to repeat the experiments in TE1. The number of corrections used in the BFGS update is set to 7. It is found that if the estimated parameters are not scaled, the minimization process cannot converge and the prescribed parameters cannot be estimated, which is similar to that in Yang and Hamrick (2005). The convergence rate is related to the Hessian matrix which is the matrix of second derivatives of the cost function with respect to the estimated parameter (Yang and Hamrick, 2005). If the Hessian matrix is well-conditioned, the L-BFGS will converge rapidly

and the parameters will be well estimated. Therefore, an efficient preconditioning method is necessary to speed up the convergence rate. As Navon et al. (1992) pointed out, the control variables in the minimization algorithm should be scaled to similar magnitudes on the order of unity because convergence tolerances and other criteria are based on an implicit definition of small and large in the optimization algorithm. Therefore, we scaled these parameters by a factor of the reciprocal of the prescribed values, respectively. After two iteration steps, the estimated sensitive parameters, including ws, M0 and Cobc, are close to the prescribed values. Furthermore, all the normalized MAEs between the prescribed and estimated



10 Iteration Steps

15

Cost function

(B)


−8

10

−16

10

0

5

10 Iteration Steps

15

Cost function

(C)


−9

10

5

10 Iteration Steps

15

KH (m2/s)

0 20

(E)

−3


0

10

x 10 2

1

−1

0

5

10 Iteration Steps

15

0 20

3

1.5

−18

10

0

1

10

160

80

10

0 20

0

10

0

10

x 10 2

1


10

−6

0

10

(D)

−1

0 20

Cost function

5

M0 (kg/m2/s)

0

1

10 Cost function

0.5

−18

10

ws (m/s)

−9

10

x 10 1

Cobc (kg/m3)

Cost function

(A)


KV (m2/s)

−4

0

10

9

0

5

10 Iteration Steps

15

0 20

Fig. 6. Same as Fig. 3, but for TE4 in which the initial guess is set to 50% of the prescribed value.

0.5

−18

5

10 Iteration Steps

15

Cost function

(B)


−8

10

−16

10

0

5

10 Iteration Steps

15

Cost function

(C)


−9

10

0

5

10 Iteration Steps

15

0 20 −3

x 10 2

1

10

(E)

0

10

1 Normalized cost function Inverted KV

−1

0 20

0

10

10

0

5

10 Iteration Steps

15

0 20

3

1.5

−18

10

x 10 2

1

80 Normalized cost function Inverted KH

10 −6

0

10

0

10

−1

0 20

Cost function

0

M0 (kg/m2/s)

10

160

(D)

KH (m2/s)

−8

10

1

10 Cost function


ws (m/s)

(A)

Cobc (kg/m3)

Cost function

10

x 10 1

KV (m2/s)

−4

0

0

5

10 Iteration Steps

15

0 20

Fig. 7. Same as Fig. 3, but for TE5 in which the initial guess is set to 150% of the prescribed value.

parameter values, which is normalized by its initial value, are less than 107. On the contrary, the insensitive parameters, including KH and KV, are still not estimated successfully. The results indicate that the L-BFGS is more effective than the steepest descent method. In addition, the same conclusion, i.e., the sensitive parameters are estimated more easily than the insensitive parameters, is obtained using the L-BFGS. 5. Discussion The dynamics of sediment transport are of great importance for coastal waters. The parameters in sediment transport models are

characterized by complexity and large quantities, especially when they are assumed to be spatially or temporally varying. For most parameters, a better understanding is required to accurately map SSCs. Data assimilation methods provide possibilities to achieve this purpose by assimilating observations from multiple sources, including ship-based measurements and remote sensing data. However, few sediment transport models have been developed with data assimilation technology. The model developed in YH03 is one of them, but it only estimates constant settling velocity and resuspension rate. In fact, other parameters are also very important for coastal waters, such as the SSC along the river mouth and open boundaries. Therefore, the model developed in the present work


10


tried to include and to estimate all kinds of parameters by assimilating not only conventional SSC observations but also data from new technologies. Additionally, the continuous adjoint model was set up in YH03, while the discrete adjoint model is developed in the present study. The latter one is in principle preferred for data assimilation problems because it provides the exact derivatives of the discrete function being minimized (Sandu et al., 2005). Besides, further application of the model in YH03 was not found. In this work, through the parameter sensitivity analysis and the estimation of constant parameters, we find several sensitive parameters which are more important for the present sediment transport model. Based on the conclusions of this work, we will try to estimate the space-dependent and time-dependent parameters by assimilating the real observations in the future, which is very important for accurate simulation of SSC. In addition, the parameters in the sediment transport model have upper and lower limits, a constrained minimization algorithm, such as interior penalty and barrier methods, will improve the efficiency of the minimization procedure and produce better results for the estimated parameters (Zhu and Navon, 1999). Therefore, we will also try to improve the minimization algorithm in the future work. 6. Summary To improve the accuracy of SSC simulation, 4DVAR data assimilation is implemented in numerical modelling. A threedimensional sigma-coordinate cohesive sediment transport model with the adjoint method is set up in this paper. The main physical processes of cohesive sediment transport are parameterized, including advection-diffusion, deposition and erosion. Tradi-

nþ1 n Ci;j;k Ci;j;k

In the twin experiments, the constant parameters, including inflow open boundary conditions, suspended sediment settling velocity, resuspension rate, and the horizontal and vertical diffusivity coefficients, are estimated with different step lengths or different initial guesses using the steepest descent method. The results show that all sensitive parameters can be estimated successfully in all of the experiments. Therefore, the reasonability and feasibility of the model developed in this paper are verified. Conversely, the horizontal and vertical diffusivity coefficients, which are insensitive parameters, can be inverted successfully only in specific experiment. In addition, these parameters are estimated using the L-BFGS method and the similar results are obtained. Therefore, we conclude that the sensitive parameters are estimated more easily than the insensitive parameters. This conclusion will benefit the practical application of this model in future work. Acknowledgements We would like to deeply thank the reviewers for the constructive suggestions to greatly improve the manuscript. Thanks are extended to Professor Jorge Nocedal at Northwestern University for sharing the source code of L-BFGS. This work is supported by the National Natural Science Foundation of China through grant 41206001, the Natural Science Foundation of Zhejiang Province through Grant LY15D060001, and the State Key Laboratory of Marine Geology, Tongji University. Appendix A A.1. Numerical scheme of Eq. (1), Eq. (2) and Eq. (5).

nþ1 nþ1 n n wnþ1 Ci;j;kþ1 Ci;j;k1 Ci;jþ1;k Ci;j1;k i;j;k þ vni;j;k $ þ nþ1 $ Dt 2Dx 2Dy 2Ds Hi;j ! n n n n n n n n 1 ðKH Þiþ1;j;k þ ðKH Þi;j;k Ciþ1;j;k Ci;j;k ðKH Þi;j;k þ ðKH Þi1;j;k Ci;j;k Ci1;j;k $ $ Dx 2 Dx 2 Dx ! n n n n n n n n 1 ðKH Þi;jþ1;k þ ðKH Þi;j;k Ci;jþ1;k Ci;j;k ðKH Þi;j;k þ ðKH Þi;j1;k Ci;j;k Ci;j1;k $ $ Dy 2 Dy 2 Dy

þ uni;j;k $

n n Ciþ1;j;k Ci1;j;k

nþ1 nþ1 nþ1 nþ1 nþ1 Ci;j;kþ1 Ci;j;k C nþ1 Ci;j;k1 ðKV Þnþ1 ðKV Þnþ1 i;j;kþ1 þ ðKV Þi;j;k i;j;k þ ðKV Þi;j;k1 i;j;k $ $ nþ1 nþ1 nþ1 2 2 Hi;j Ds Hi;j Ds Hi;j Ds

1

(A.1) !

nþ1 nþ1 Ci;j;kþ1 Ci;j;k1 ðws Þnþ1 i;j;k ¼0 nþ1 $ 2Ds Hi;j

tionally, these parameters are poorly known and usually are empirically assigned. However, for this model, these parameters are estimated using the adjoint method by assimilating SSC observations. Two methods of parameter sensitivity analysis, including the relative sensitivity function and the method using the gradients of cost function with respect to the parameters, are applied to evaluate the importance of each parameter. The results of sensitivity analysis indicate that the initial conditions, inflow open boundary conditions, suspended sediment settling velocity and resuspension rate are sensitive for the present model, while the horizontal and vertical diffusivity coefficients are insensitive.

ðws Þnþ1 i;j;K $

ðws Þnþ1 i;j;1 $

nþ1 nþ1 Ci;j;K þ Ci;j;Kþ1

2 nþ1 nþ1 Ci;j;1 þ Ci;j;0

2

nþ1 nþ1 Ci;j;Kþ1 Ci;j;K ðKV Þnþ1 i;j;K $ ¼0 nþ1 Ds Hi;j

(A.2)

nþ1 nþ1 Ci;j;1 Ci;j;0 ðKV Þnþ1 i;j;1 nþ1 $ ¼ Ei;j Dnþ1 i;j nþ1 Ds Hi;j

(A.3) A.2. The discrete adjoint model.



11

nþ1 nþ1 nþ1 nþ1 nþ1 nþ1 lnþ1 lnþ1 lni;j;k1 wnþ1 lni;j;kþ1 wnþ1 lni;j;k lnþ1 i1;j;k ui1;j;k liþ1;j;k uiþ1;j;k i;j1;k vi;j1;k li;jþ1;k vi;jþ1;k i;j;k i;j;k1 i;j;kþ1 þ þ þ nþ1 Dt 2Dx 2Dy 2Hi;j Ds nþ1 nþ1 nþ1 nþ1 nþ1 nþ1 nþ1 ðKH Þnþ1 1 ðKH Þiþ1;j;k þ ðKH Þi;j;k liþ1;j;k li;j;k i;j;k þ ðKH Þi1;j;k li;j;k li1;j;k Dx 2 Dx 2 Dx

!

nþ1 nþ1 nþ1 nþ1 nþ1 nþ1 nþ1 ðKH Þnþ1 1 ðKH Þi;jþ1;k þ ðKH Þi;j;k li;jþ1;k li;j;k i;j;k þ ðKH Þi;j1;k li;j;k li;j1;k Dy 2 Dy 2 Dy

!

nþ1 n n ðKV Þnþ1 i;j;kþ1 þ ðKV Þi;j;k li;j;kþ1 li;j;k

1

nþ1 Hi;j Ds

nþ1 Hi;j Ds

2

n nþ1 lni;j;k1 ðws Þnþ1 i;j;k1 li;j;kþ1 ðws Þi;j;kþ1

nþ1 2Hi;j Ds

(A.4)

nþ1 n n ðKV Þnþ1 i;j;k þ ðKV Þi;j;k1 li;j;k li;j;k1

!

nþ1 Hi;j Ds

2

nþ1 ~ nþ1 ¼ 0 þ W Ci;j;k C i;j;k

The corresponding surface and bottom boundary conditions in the adjoint model are:

li1;j;k n li;j;k n vJ n n Ci;j;k Ci1;j;k þ C ¼ 2Ci;j;k n 2 2Dx 2Dx2 iþ1;j;k vðKH Þi;j;k n

lni;j;K ðws Þnþ1 $ i;j;K

þ

lni;j;Kþ1

2

lni;j;1 þ lni;j;0

ðws Þnþ1 i;j;1 $

2

n n ðKV Þnþ1 i;j;K li;j;Kþ1 li;j;K $ nþ1 Ds Hi;j

n n ðKV Þnþ1 i;j;1 li;j;1 li;j;0 $ ¼0 nþ1 Ds Hi;j

¼0

(A.5)

(A.6)

Furthermore, other boundary conditions and initial conditions in the adjoint model are equal to zero. A.3. The gradients of the cost function with respect to parameters.

vJ vðws Þnþ1 i;j;k

8 ! nþ1 nþ1 n nþ1 > lni;j;1 Ci;j;2 Ci;j;0 > l þ lni;j;1 nþ1 t > i;j;0 b > Ci;j;1 1 þ > > nþ1 > 2 tcd Ds 2Hi;j > > > > > > > ln þ ln C nþ1 þ C nþ1 > > i;j;0 i;j;0 i;j;1 i;j;1 > > ; k ¼ 1&tnþ1 < tcd > b > 2 2 > > > > > n nþ1 nþ1 > nþ1 nþ1 > þ Ci;j;0 lni;j;0 þ lni;j;1 Ci;j;1 < li;j;1 Ci;j;2 Ci;j;0 ; k ¼ 1&tnþ1 tcd ¼ b nþ1 2 2 Ds 2H > i;j > > > > > n > nþ1 nþ1 > C l C > i;j;k i;j;kþ1 i;j;k1 > > > ; k ¼ 2; 3; :::; K 1 > nþ1 > > Ds 2Hi;j > > > > > > nþ1 nþ1 > nþ1 nþ1 > lni;j;K Ci;j;Kþ1 Ci;j;K1 þ Ci;j;Kþ1 lni;j;K þ lni;j;Kþ1 Ci;j;K > > > ; k¼K > : nþ1 2 2 Ds 2Hi;j

n

ln ln iþ1;j;k i;j1;k n n n n Ciþ1;j;k þ Ci;j;k þ Ci1;j;k Ci;j;k 2 2 2Dx 2Dy lni;j;k n n n n þ Ci;jþ1;k 2Ci;j;k Ci;j1;k þ Ci;j1;k 2 2Dy lni;jþ1;k n n Ci;jþ1;k Ci;j;k 2Dy2 (A.8)

(A.7)


12


vJ vðKV Þnþ1 i;j;k

¼

8 > lni;j;0 lni;j;1 > nþ1 nþ1 nþ1 nþ1 nþ1 > > C þ C C 2C þ C > i;j;1 i;j;0 i;j;2 i;j;1 i;j;0 2 2 > > nþ1 nþ1 > 2 Hi;j Ds 2 Hi;j Ds > > > > > > > ln þ ln > ln > i;j;0 i;j;1 nþ1 nþ1 nþ1 nþ1 > > i;j;2 Ci;j;2 C k¼1 C C > i;j;1 i;j;1 i;j;0 ; 2 nþ1 > > 2Hi;j Ds nþ1 > > 2 Hi;j Ds > > > > > > > ln lni;j;k > nþ1 nþ1 nþ1 nþ1 nþ1 > > i;j;k1 Ci;j;k þ C 2 Ci;j;kþ1 2Ci;j;k þ Ci;j;k1 > i;j;k1 2 > > nþ1 > Ds 2 Hnþ1 Ds > lni;j;kþ1 > nþ1 nþ1 > > 2 Ci;j;kþ1 Ci;j;k ; k ¼ 2; 3; :::; K 1 > > > nþ1 > > > 2 Hi;j Ds > > > > > > lni;j;K1 nþ1 lni;j;K > nþ1 nþ1 nþ1 nþ1 > > Ci;j;K Ci;j;K1 þ Ci;j;Kþ1 2Ci;j;K þ Ci;j;K1 > 2 2 > > nþ1 nþ1 > 2 Hi;j Ds 2 Hi;j Ds > > > > > > > ln þ ln > lni;j;Kþ1 nþ1 > i;j;K i;j;Kþ1 nþ1 nþ1 nþ1 > > C C ; k¼K C C > i;j;Kþ1 i;j;K i;j;Kþ1 i;j;K 2 nþ1 > 2Hi;j Ds > : 2 Hnþ1 Ds i;j

vJ vðM0 Þnþ1 i;j

vJ vðCobc Þni;j;k

! 8 n l þ lni;j;1 tnþ1 > > b < i;j;0 1 ; 2 tce ¼ > > : 0; tnþ1 tce b

tnþ1 > tce b

(A.10)

8 n n n n n ðK l Þ þ ðK Þ > l u H H iþ1;j;k i;j;k iþ1;j;k > iþ1;j;k iþ1;j;k > > þ ; for west > > 2Dx > 2Dx2 > > > > n > n n > > lni1;j;k uni1;j;k li1;j;k ðKH Þi;j;k þ ðKH Þi1;j;k > > > þ ; for east > < 2D x 2Dx2 ¼ n n n > n > li;jþ1;k vni;jþ1;k li;jþ1;k ðKH Þi;j;k þ ðKH Þi;jþ1;k > > > þ ; for south > > > 2Dy 2Dy2 > > > > > > n n n > > lni;j1;k vni;j1;k li;j1;k ðKH Þi;j;k þ ðKH Þi;j1;k > > : ; for north 2Dy 2Dy2

(A.11)

0 0 0 0 l0i;j;k li1;j;k u0i1;j;k liþ1;j;k u0iþ1;j;k li;j1;k u0i;j1;k li;jþ1;k u0i;jþ1;k 1 ðKH Þ0iþ1;j;k þ ðKH Þ0i;j;k l0iþ1;j;k l0i;j;k vJ þ ¼ vðC0 Þi;j;k Dt 2Dx 2Dy Dx 2 Dx ! ! 0 0 0 0 0 0 0 0 ðKH Þ0i;j;k þ ðKH Þ0i1;j;k l0i;j;k l0i1;j;k 1 ðKH Þi;jþ1;k þ ðKH Þi;j;k li;jþ1;k li;j;k ðKH Þi;j;k þ ðKH Þi;j1;k li;j;k li;j1;k þ 2 Dx Dy 2 Dy 2 Dy

(A.12)

References Agrawal, Y.C., Pottsmith, H.C., 2000. Instruments for particle size and settling velocity observations in sediment transport. Mar. Geol. 168, 89e114. Alekseev, A.K., Navon, I.M., Steward, J.L., 2009. Comparison of advanced large-scale

minimization algorithms for the solution of inverse ill-posed problems. Optim. Methods & Softw. 24, 63e87. Amoudry, L.O., 2014. Extension of k-u turbulence closure to two-phase sediment transport modelling: application to oscillatory sheet flows. Adv. Water Resour. 72, 110e121.


D. Wang et al. / Estuarine, Coastal and Shelf Science xxx (2016) 1e14 Amoudry, L.O., Souza, A.J., 2011. Deterministic coastal morphological and sediment transport modeling: a review and discussion. Rev. Geophys. 49, RG2002. Bagnold, R.A., 1935. The movement of desert sand. Geogr. J. 85, 342e365. Bagnold, R.A., 1937. The size-grading of sand by wind. Proc. R. Soc. Lond. Ser. A, Math. Phys. Sci. 250e264. Blumberg, A.F., Mellor, G.L., 1987. A description of a three-dimensional coastal ocean circulation model. In: Three-dimensional Coastal Ocean Models, pp. 1e16. Cacuci, D.G., Ionescu-Bujor, M., Navon, I.M., 2005. Sensitivity and Uncertainty Analysis, Volume II: Applications to Large-scale Systems. CRC press, Boca Raton, Florida, USA. Cai, L., Tang, D., Li, X., Zheng, H., Shao, W., 2015. Remote sensing of spatial-temporal distribution of suspended sediment and analysis of related environmental factors in Hangzhou Bay, China. Remote Sens. Lett. 6, 597e603. Cao, A., Gao, Y., Zhang, J., Lv, X., 2013. Trajectory estimation of aircraft in a doublesatellite passive positioning system with the adjoint method. Math. Problems Eng. 2013. Article ID 502610. Chen, C., Liu, H., Beardsley, R.C., 2003. An unstructured grid, finite-volume, threedimensional, primitive equations ocean model: application to coastal ocean and estuaries. J. Atmos. Ocean. Technol. 20, 159e186. Cloern, J.E., 1987. Turbidity as a control on phytoplankton biomass and productivity in estuaries. Cont. Shelf Res. 7, 1367e1381. Dennison, W.C., Orth, R.J., Moore, K.A., Stevenson, J.C., Carter, V., Kollar, S., Bergstrom, P.W., Batiuk, R.A., 1993. Assessing water quality with submersed aquatic vegetation. BioScience 43, 86e94. Du, P., Ding, P., Hu, K., 2010. Simulation of three-dimensional cohesive sediment transport in Hangzhou Bay, China. Acta Oceanol. Sin. 29, 98e106. Dyer, K., 1986. Coastal and Estuarine Sediment Dynamics. Wiley, Chichester. Einstein, H.A., 1950. The Bed-Load Function for Sediment Transportation in Open Channel Flows. US Department of Agriculture. Fan, D., Tu, J., Shang, S., Cai, G., 2014. Characteristics of tidal-bore deposits and facies associations in the Qiantang Estuary, China. Mar. Geol. 348, 1e14. Fan, W., Lv, X., 2009. Data assimilation in a simple marine ecosystem model based on spatial biological parameterizations. Ecol. Model. 220, 1997e2008. Fettweis, M., Nechad, B., Van den Eynde, D., 2007. An estimate of the suspended particulate matter (SPM) transport in the southern North Sea using SeaWiFS images, in situ measurements and numerical model results. Cont. Shelf Res. 27, 1568e1583. Gibbs, R.J., 1985. Estuarine flocs: their size, settling velocity and density. J. Geophys. Res. Oceans 90 (C2), 3249e3251. Guan, W.B., Kot, S.C., Wolanski, E., 2005. 3-D fluid-mud dynamics in the Jiaojiang Estuary, China. Estuarine. Coast. Shelf Sci. 65, 747e762. Hamrick, J.M., 1992. A Three-dimensional Environmental Fluid Dynamics Computer Code: Theoretical and Computational Aspects. Virginia Institute of Marine Science, College of William and Mary. Hayter, E.J., 1983. Prediction of Cohesive Sediment Transport in Estuarial Waters. Unpublished PhD Dissertation. University of Florida, Gainesville, Florida. Hayter, E.J., Mehta, A.J., 1986. Modelling cohesive sediment transport in estuarial waters. Appl. Math. Model. 10, 294e303. He, X., Bai, Y., Pan, D., Huang, N., Dong, X., Chen, J., Chen, C.-T.A., Cui, Q., 2013. Using geostationary satellite ocean color data to map the diurnal dynamics of suspended particulate matter in coastal waters. Remote Sens. Environ. 133, 225e239. Holly, F.M., 1990. Charima: Numerical Simulation of Unsteady Water and Sediment Movement in Multiply Connected Networks of Mobile-bed Channels. Iowa Institute of Hydraulic Research, the University of Iowa. Hu, K., Ding, P., Wang, Z., Yang, S., 2009. A 2D/3D hydrodynamic and sediment transport model for the Yangtze Estuary, China. J. Mar. Syst. 77, 114e136. Jia, J., Wang, Y., Gao, S., Wang, A., Li, Z., 2006. Interpreting grain-size trends associated with bedload transport on the intertidal flats at Dafeng, central Jiangsu coast. Chin. Sci. Bull. 51, 341e351. Jia, Y., Wang, S.S., 1999. Numerical model for channel flow and morphological change studies. J. Hydraul. Eng. 125, 924e933. Krishnappan, B.G., 1981. Users Manual: Unsteady, Non-uniform, Mobile Boundary Flow Model-MOBED. Hydraulic Division, National Water Research Institute, CCIW, Burlington, Ontario, Canada. Krone, R.B., 1962. Flume Studies of the Transport of Sediment in Estuarial Shoaling Processes. U.C. Berkrley, Berkeley. Leupi, C., Altinakar, M.S., Deville, M., 2008. Numerical modeling of cohesive sediments dynamics in estuaries: Part IdDescription of the model and simulations in the Po River Estuary. Int. J. Numer. Methods Fluids 57, 237e263. Li, X., Wunsch, C., 2004. An adjoint sensitivity study of chlorofluorocarbons in the North Atlantic. J. Geophys. Res. Oceans 109, C01007. Liu, D.C., Nocedal, J., 1989. On the limited memory method for large scale optimization. Math. Program. 45, 503e528. Martin, J.M., Windom, H.L., 1991. Present and future roles of ocean margins in regulating marine biogeochemical cycles of trace elements. Ocean Margin Process. Glob. Change 1991, 45e67. May, C.L., Koseff, J.R., Lucas, L.V., Cloern, J.E., Schoellhamer, D.H., 2003. Effects of spatial and temporal variability of turbidity on phytoplankton blooms. Mar. Ecol. Prog. Ser. 254, 111e128. Mayer, L.M., Keil, R.G., Macko, S.A., Joye, S.B., Ruttenberg, K.C., Aller, R.C., 1998. Importance of suspended participates in riverine delivery of bioavailable nitrogen to coastal zones. Glob. Biogeochem. Cycles 12, 573. McLaughlin, C.J., Smith, C.A., Buddemeier, R.W., Bartley, J.D., Maxwell, B.A., 2003. Rivers, runoff, and reefs. Glob. Planet. Change 39, 191e199.

13

Mehta, A.J., Hayter, E.J., Parker, W.R., Krone, R.B., Teeter, A.M., 1989. Cohesive sediment transport. I: process description. J. Hydraul. Eng. 115, 1076e1093. Miller, R.L., McKee, B.A., 2004. Using MODIS Terra 250 m imagery to map concentrations of total suspended matter in coastal waters. Remote Sens. Environ. 93, 259e266. Myint, S.W., Walker, N.D., 2002. Quantification of surface suspended sediments along a river dominated coast with NOAA AVHRR and SeaWiFS measurements: Louisiana, USA. Int. J. Remote Sens. 23, 3229e3249. Navon, I.M., 1998. Practical and theoretical aspects of adjoint parameter estimation and identifiability in meteorology and oceanography. Dyn. Atmos. Oceans (Special Issue honor Richard Pfeffer) 27, 55e79. Navon, I.M., Zou, X., Derber, J., Sela, J., 1992. Variational data assimilation with an adiabatic version of the NMC spectral model. Mon. Weather Rev. 120, 1433e1446. Neukermans, G., Ruddick, K., Bernard, E., Ramon, D., Nechad, B., Deschamps, P.-Y., 2009. Mapping total suspended matter from geostationary satellites: a feasibility study with SEVIRI in the Southern North Sea. Opt. Express 17, 14029e14052. Neukermans, G., Ruddick, K.G., Greenwood, N., 2012. Diurnal variability of turbidity and light attenuation in the southern North Sea from the SEVIRI geostationary sensor. Remote Sens. Environ. 124, 564e580. Onishi, Y., Wise, S.E., 1982. SERATRA: User's Manual for the Instream Sedimentcontaminant Transport Model. Technical Rep. No. PB-83e122739 (Pacific Northwest Laboratory, Battelle-Northwest, Richland, Wash). Owen, M.W., 1970. A Detailed Study of the Settling Velocities of an Estuary Mud. Hydraulics Research Station. Papanicolaou, A.N., Bdour, A., Wicklein, E., 2004. A numerical model for the study of sediment transport in steep mountain streams. J. Hydraul. Res. 42, 357e366. Papanicolaou, A.N., Elhakeem, M., Krallis, G., Prakash, S., Edinger, J., 2008. Sediment transport modeling reviewdcurrent and future developments. J. Hydraul. Eng. 134, 1e14. Petus, C., Chust, G., Gohin, F., Doxaran, D., Froidefond, J.M., Sagarminaga, Y., 2010. Estimating turbidity and total suspended matter in the Adour River plume (South Bay of Biscay) using MODIS 250-m imagery. Cont. Shelf Res. 30, 379e392. Salama, M.S., Shen, F., 2010. Simultaneous atmospheric correction and quantification of suspended particulate matters from orbital and geostationary earth observation sensors. Estuarine. Coast. Shelf Sci. 86, 499e511. Sandu, A., Daescu, D.N., Carmichael, G.R., Chai, T., 2005. Adjoint sensitivity analysis of regional air quality models. J. Comput. Phys. 204, 222e252. Sanford, L.P., Halka, J.P., 1993. Assessing the paradigm of mutually exclusive erosion and deposition of mud, with examples from upper Chesapeake Bay. Mar. Geol. 114, 37e57. Sheng, Y.P., Eliason, D.E., Chen, X.J., 1992. . Modeling three-dimensional circulation and sediment transport in lakes and estuaries. Estuar. Coast. Model. 1991, 105e115. ASCE. Sirkes, Z., Tziperman, E., 1997. Finite difference of adjoint or adjoint of finite difference? Mon. Weather Rev. 125, 3373e3378. Song, Y., Haidvogel, D., 1994. A semi-implicit ocean circulation model using a generalized topography-following coordinate system. J. Comput. Phys. 115, 228e244. Thomas, W.A., Prasuhn, A.L., 1977. Mathematical modeling of scour and deposition. J. Hydraul. Div. 103, 851e863. Van Rijn, L.C., Tan, G.L., 1985. Sutrench Model: Two-dimensional Vertical Mathematical Model for Sedimentation in Dredged Channels and Trenches by Currents and Waves. Rijkswaterstaat. Wang, X.H., Pinardi, N., 2002. Modeling the dynamics of sediment transport and resuspension in the northern Adriatic Sea. J. Geophys. Res. Oceans 107 (C12), 3225. Wang, Y.P., Voulgaris, G., Li, Y., Yang, Y., Gao, J., Chen, J., Gao, S., 2013. Sediment resuspension, flocculation and settling in a macrotidal estuary. J. Geophys. Res. Oceans 118, 5591e5608. Warrick, J.A., Mertes, L.A.K., Siegel, D.A., Mackenzie, C., 2004. Estimating suspended sediment concentrations in turbid coastal waters of the Santa Barbara Channel with SeaWiFS. Int. J. Remote Sens. 25, 1995e2002. Yang, Z., Hamrick, J.M., 2003. Variational inverse parameter estimation in a cohesive sediment transport model: an adjoint approach. J. Geophys. Res. Oceans 108 (C2), 3055. Yang, Z., Hamrick, J.M., 2005. Optimal control of salinity boundary condition in a tidal model using a variational inverse method. Estuarine. Coast. Shelf Sci. 62, 13e24. Yu, L., O'Brien, J.J., 1992. On the initial condition in parameter estimation. J. Phys. Oceanogr. 22, 1361e1364. Zhang, J., Lu, X., 2010. Inversion of three-dimensional tidal currents in marginal seas by assimilating satellite altimetry. Comput. Methods Appl. Mech. Eng. 199, 3125e3136. Zhang, J., Lu, X., Wang, P., Wang, Y.P., 2011. Study on linear and nonlinear bottom friction parameterizations for regional tidal models using data assimilation. Cont. Shelf Res. 31, 555e573. Zhang, J., Wang, Y.P., 2014. A method for inversion of periodic open boundary conditions in two-dimensional tidal models. Comput. Methods Appl. Mech. Eng. 275, 20e38. Zhang, M., Tang, J., Dong, Q., Song, Q., Ding, J., 2010. Retrieval of total suspended matter concentration in the Yellow and East China Seas from MODIS imagery. Remote Sens. Environ. 114, 392e403.


14


Zhu, Y., Navon, I.M., 1999. Impact of parameter estimation on the performance of the FSU Global Spectral Model using its full-physics adjoint. Mon. Weather Rev. 127, 1497e1517. Zou, J., Hsieh, W.W., Navon, I.M., 1995. Sequential open-boundary control by data assimilation in a limited-area model. Mon. Weather Rev. 123, 2899e2909. Zou, X., Barcilon, A., Navon, I.M., Whitaker, J., Cacuci, D.G., 1993a. An adjoint

sensitivity study of blocking in a two-layer isentropic model. Mon. Weather Rev. 121, 2833e2857. Zou, X., Navon, I.M., Berger, M., Phua, K.H., Schlick, T., Le Dimet, F.X., 1993b. Numerical experience with limited-memory quasi-Newton and truncated Newton methods. SIAM J. Optim. 3, 582e608.
