Inversion Study on Pollutant Discharges in the Bohai Sea withthe Adjoint Method

SHEN Youli, WANG Chunhui, WANG Yonggang, and LV Xianqing), *

?

Inversion Study on Pollutant Discharges in the Bohai Sea withthe Adjoint Method

SHEN Youli1), 2), WANG Chunhui3), WANG Yonggang4), and LV Xianqing1), *

1),,266100,2),536000,3),,266033,4),,266061,

The adjoint methodis presented which could be used to estimate the time-varying contamination concentration (CC) from pollution source (PS). Then the pollutant flux is calculated.In order to guarantee the continuity of pollutant distribution and make the calculated results more accurate,the independent point (IP)scheme is proposed. The contamination concentrations (CCs)atsome time steps are selected as the independent points (IPs), and only CCs at these IPs are optimized while CCsat other points are calculated through linear interpolation of the independent CCs.In twin numerical experiments, all the given distributions are successfully inverted with the adjoint method. The cost functions and the mean absolute errors (MAEs) in concentrations and pollutant fluxes decrease greatly after assimilation,and the cost functions are reduced by about 5 orders of magnitude compared with their initial values. The results indicate that the adjoint method is computationally efficient to recover CCs from PS. It is easier to invert the given distribution which is less complex.The inversion efficiency with IP scheme is raised compared to that without this scheme. The IP scheme is significant for the inversion result, in which appropriate IP number couldyieldbetter inversion results. More work will be done to apply this method to real experiment.

adjoint method; time-varying variable; contamination concentration; pollutant flux; independent point scheme

1 Introduction

The Bohai Sea is China’s only inland sea, which receives large volumes of domestic and industrial discharges and even accidental or intentional spills from moving ships. According to statistics, there are more than 80 rivers flowing into the Bohai Sea, including 40 rivers which have important effects on the Bohai Sea. Due to limited flushing capacity, the multiple pollutant discharges in the Bohai Sea can deteriorate the water quality significantly (Zheng, 2011). Internationally there is a clear consensus that source pollution, especially land source pollution, is having a profound effect on coastal and marine ecosystems(Waterhouse, 2012). Preventing point source pollution is among the most cost-effective means of marine management.

Obtaining precise information on pollutant discharge from point sources such as wastewater effluent and runoff from oil fields and unsewered industrial sites, is important forpollution prediction and effective control. However, the information is hard to obtainbecause of the complex variation of sources, especially for large-scale and natural sources (Zhu and Wang, 2006). One way to solve this problem is to study variety of sources and to construct pollutant dischargeor emission model. Gupta(2004) carried out water quality simulation of Thane creek using a two-dimensional advection-diffusion model to determine the wastewater assimilative capacity of the creek considering 11 point-source discharges. Shulka(2002) considered advection and reaction when studying the steady nonconservative pollutant with time-dependent periodic waste discharge concentration. Milnes and Perrochet (2007) predicted the concentration at the unknown PS when a reversed flow field transport simulation was performed. However, with limited observations, this approachusually causes the uncertainty and even deviation from the reality.

Another approach is to invert the pollutant discharges or emissions using the observation and the pollutant discharge or emission model.This approachmainly includes two categories: gradient based methods and gradient free methods. In the first category, there is variational method based on an adjoint, reduced adjoint, simple finite difference method and stochastic estimation of gradient. For thesecond category, there is ensemble Kalman filter (EnKF), Bayesian methodsand sequential Monte-Carlo methods like particle filtering. Gilliland(2001) used the discrete Kalman filter (DKF) in an adaptive-iterative mode to deduce time-varying air quality emission.Their result shows that if the initial condition assumption is violated, the adaptive-iterative DKF can produce biased emissions to compensatefor the initial modeled and observed concentration difference. Heeminkand Segers (2002) used the Kalman filter to estimate and predict ozone concentrations in a part of North West Europe. Their work focused onthe development of fast filter algorithms to make Kalman filtering feasible for high dimensional space-time models. Zhu and Wang (2006) presented some ensemble based statistical estimation methods for inverting modeling of pollution emissions. The results show that EnKF can estimate the time-variant emission effectively at every time step when the observations are not available. Bousquet(1999) presented a 3-D Bayesian inverse methodology to infer annual net carbon dioxide (CO2) sources and sinks using atmospheric measurements of CO2. Elbern(2000) studied emissions using a regulatory air quality model with four dimension variational (4DVAR) data assimilation. It is shown that the space-time variational method is able to analyze emission rates of Nitric oxide (NO) directly, as for volatile organic compounds (VOC), regularization techniques must be introduced. 4DVAR has the advantage of directly assimilating variousobservations distributed in time and space into numerical models while maintaining dynamicaland physical consistency with the model. In this paperthe adjoint method is used to invert CC from point source and the pollutant flux is also computed. The CC in this paper is taken as a time-varying variable, rather than a constant in traditional studies.

The paper is organized as follows. In the next section the model is presented, leading to the relationship between the target concentrations and the flow field. In Section 3 the method we used in this paper is illustrated. The selectionand the optimization of the independent CCs areshown in Section 4. In Section 5 numericalexperiments arecarried out and the results are analyzed. Finally conclusions are presented in Section 6.

2 Model Introduction

2.1 The Background Flow Field

The background flow field is provided by the 3-D Regional Ocean Model System (ROMS). ROMS is a free-surface, hydrostatic, and primitive equation ocean model that uses stretched, terrain-following coordinates in the vertical and orthogonal curvilinear coordinates in the hori- zontal. It has been expanded to include a variety of new features including high-order advection schemes, accurate pressure gradient algorithms, several subgrid-scale param- eterizations, atmospheric, oceanic, benthic boundary layers, biological modules, and radiation boundary conditions (Shchepetkin and McWilliams, 2005).

In the present study the horizon resolution, the vertical resolutionand the integral time step is similar to Wang(2013). The open boundary is installed on 122.5?E. The open boundary condition is from the measured data.Surface forcing functions are drawn from the COADS climatology of Da Silva. (1994). These forcing fields include surface wind stresses, heat and freshwater fluxes, and heat flux sensitivity to sea surface temperature. The model is run for 10 years and result of the last year serves as the hydrodynamic background field for ecological model. The simulated circulation is Eulerian residual current. Mean surface and bottomcirculation of the Bohai Sea in July are shown in Fig.1. As for the mean surface circulation, the maximum current is about8.0cms?1, while the minimum current is less than 2.0cms?1.A counterclockwise gyre is in the northeast corner of the Liao- dong Bay and the flow direction is from south to north in the Bohai Bay. In the Laizhou Bay the flow becomes weaker. Meanwhile the circulation near the Bohai Strait is inflow in the north and outflow in the southern part of the Bohai Strait, which is similar to Wei(2001) and Wei(2003). While the bottom circulation is much weaker than that of surface except in the Bohai Strait, and the maximum current is less than 1.0cms?1.

Fig.1 Mean surface and bottom circulation of the Bohai Sea in July. (a) is the mean surface circulation, and (b) is the mean bottom (the third layer) circulation.

2.2 Governing Equation and Its Difference Scheme

Based on the hydrodynamic background field, an advection-diffusion model is constructed, and the basic equation system is written in the following form:

, (1)

whereis the concentration of contamination,,,are components of the Cartesian coordinate system,,,are velocities in the Cartesian directions,is time,AandKare the horizontal and vertical turbulent diffusion coefficients, respectively.

The petroleum pollutants are considered as pollutants for example, and the initial concentration is set to be 25μgL?1. Radiation boundary condition is given on the open boundary, which is described as:

wheredenotes the normal direction of the open boundary. While the closed wall boundary condition is in the following form:

. (3)

The difference equation of governing Eq. (1) is written as follows:

, (4)

where,

, (6)

, (7)

The stabilization condition is:

where,,are components of the Cartesian coordinate system, ?x, ?(where ?is constant),?zare grid spacings in the Cartesian directions, respectively.

2.3 The Adjoint Equation and Its Difference Scheme

The cost function which is usually used to measure the misfit between the observations and the simulation results is defined as follows:

where Tdenotes transposition andthe observations.Kis the weight of the observation, which is defined as:

. (11)

By introducing Lagrangianfunction, the adjoint model, which provides a method of calculating the gradient of cost function with respect to CC from PS, can be constructed.is defined as:

, (12)

where*is called adjoint variable of. In order to minimize, let ?/?=0, then the adjoint equation can be obtained:

. (13)

The corresponding open boundary condition is defined as:

While the closed wall boundary condition is in the fol- lowing form:

. (15)

The difference equation of adjoint Eq. (13) is asfollows:

, (16)

where

, (18)

, (19)

2.4 Test of the Adjoint Model

It is essential to verify the adjoint model before proceeding with the minimization runs. As the work of Zhang and Lv (2008) and Chen(2013), let

be a Taylor expansion of the cost function, which is defined in Section 2.3, whereis a general value of CC,is a small scalar andis an arbitrary vector of unit length. Rewriting (21), a function ofcan be defined as

. (22)

As mentioned in the work of Smedstad and O’Brien (1991), Das and Lardner (1991), Lardner and Song (1995) and Navon(1991), ifis chosen close to machine zero, one cannot expect to be able to verify that the correct gradient has been found. For values ofwhich are not too close to the machine zero, one should expect to obtain a value ofΦ() which is close to 1. In Fig.2, the values Φ() are plotted when the CCs are taken as an example. In this test, the basic pointis set to be 25μgL?1, which is equal to the initial values of CCs in the following numerical experiments. It shows that forbetween 10?1and 10?8, Φ() is close to 1 (the dashed line), which means the adjoint model is verified.

Fig.2 Values of the function Φ(α)versus α.

3 The Adjoint Method

The variational method with adjoint equations was first demonstrated by Penenko and Obraztsov (1976) and has been used in manyprevious studies(Lv and Zhang, 2006; Seiler, 1993; Fan and Lv, 2009; Wang, 2013; Cao., 2013; Chen., 2013). It is an advanced data assimilationmethod that involves the adjoint method and has the advantage of directly assimilating variousobser- vations distributed in time and space into numerical mod- els while maintaining dynamicaland physical consistency with the model (Zhang and Lv, 2008). The basic idea of variational adjoint method is to define a cost function in a desirable way that measures the misfit of the model con- sidered (Robertson and Langner, 1998).

In this study, the calculation process of the adjoint method (Fig.3) is as follows:

1) An initial value of the CCs from PS is given empiri- cally;

2) Perform the simulation by running the forward model, and the simulation results are obtained;

3) The adjoint model is driven backward in time with the misfit between the model result and the observation;

4) Calculate the gradient of cost function with respect to the CCs from PS;

5) Adjust the CCs from PS on the basis of the gradient above;

6) Return to step (1) while the initial values are updated with the aboveconcentrations, and repeat the iteration procedure;

7) The process is stopped when the number of iterations is 200 (We tried different numbers of iterationsand found that 200 is a good choice mainly based on the following considerations: (1) The cost functions and the MAEs are almost no longer falling after 200 numbers of iterations; (2) The calculation required is acceptable for the 200 numbers of iterations.) or the target<(whereis a small number, such as 0.0001) is reached.

Fig.3 The calculation process of the adjoint method.

4 Twin Experiment (TE)

4.1 Experimental Design

The calculated regionis the Bohai Sea (37?N–41?N, 117.5?E–122.5?E). The horizontal resolution of the model is 4?×4?, and the open boundary is installed on 122.5?E. In the vertical direction the water is divided into 6 layers, and the thicknesses of each from top to bottom are 10m, 10m, 10m, 20m, 25m and 25m, respectively. The integral time step is 1 hour and the integral period is 7 days.

A bathymetry map of the Bohai Sea is shown in Fig.4, in which the PS and the observations are also located. Observations are marked from 1 to 20, Nos. 1–4 are observed once per day during 8:00–18:00from the second to the fourth day, Nos. 5–12 are observed once per day during 8:00–18:00 on the fourth, the fifth and the seventh day, and Nos. 13–20 are observed once per day during 8:00–18:00on the sixth and seventh day.

Fig.4 The topography of the Bohai Sea and the locations of the observation sites. ‘+’is the location of the PS, and the dots are the observations.

4.2 Correction of CCs at IPs

In order to guarantee the continuity of pollutant distribution and make the calculated results more accurate,the IP scheme is proposed.is the value of CC at the IPs,is the CC from the PS after linear interpola-tion (denotesthe calculated point,denotesthe IP), and the relationship between them is as follows:

wherewis the weight coefficient in the Cressman form (Cressman, 1959) and the specific form is as follows:

, where, (24)

whereis the influence radius,ris the distance from the IP to thecalculated point.

After derivation, the gradient of cost functionwithrespecttois

, (26)

, (28)

, (29)

Since the variable is adjusted in the inverse direction of the gradient above, thecalibration relationship of independent CC from PS is as follows:

, (31)

, (32)

4.3 The Given Variation of CC from PS

Considering the real variation of CC from PS, three types of time-varying CCs are prescribed as follows:

Type 1:; (34)

Type 2:; (35)

Type 3:. (36)

4.4 Design Strategies of IP

Six IP strategies are designed to invert the three types of prescribed CCs (Fig.5). The strategies are presented as follows:

Strategy (A): Three evenly distributed grid points are taken as IPs, and the influence radius is the length of 72 grids.

Strategy (B) to Strategy (F): The numbers of IPs are 5, 7, 9, 13 and 19, respectively, which are distributed evenly. Radiuses are the lengths of 38, 24, 18, 12 and 8 grids corresponding to 5 strategies, respectively.

All strategies are used to invert CCs with the distributions of Type 1 and 2, and Strategies (B)–(F) are used to invert CCs with the distributions of Type 3.

Fig.5 Six IP strategies. Larger points and smaller ones represent IPs and other grid points, respectively.

5 Experimental Results and Discussion

5.1 The Inversion Result of Type 1

The cost functions and MAEsin CCs and pollutant fluxes are two important convergence criteria for data assimilation in this model. The errors after assimilation for Type 1 are shown in Table 1, in which ‘None’ denotes the result without IP scheme. The results show that the cost functions, and the MAEs inCCs and pollutant fluxeswith strategies (A)–(F) decrease significantly, and the cost functions are reduced by about 5 orders of magnitude compared with their initial values, indicating that the given time-varying distribution is successfully inverted with these strategies. The result with Strategy (A) is the most satisfactory.The MAE in CC declines by 97.7% from 14.98μgL?1to 0.41μgL?1. Meanwhile the MAE in pollutant flux decreases by 98.9%from 22.74t to 0.25t. The MAEs with ‘None’ strategy are generally bigger than those with IP strategies except strategy (C), which shows that IP scheme is a more effective technique to invert the prescribed CCs.

Table 1 Errors after assimilation

The inverted distributions of CCs and pollutant fluxes with using Strategy (A) are shown in Fig.6. In this figure, the inverted distribution of CCs shows a good agreement with the prescribed one, and the maximum deviation is 1.54μgL?1. The pollutant flux is also well inverted, and the error at any moment is within 1.73t. It is concluded that the given time-varying distribution is successfully inverted using Strategy (A).

Fig.6 The given distribution (solid lines) and the corresponding inversion results (dotted lines). (a) is the given time-varying concentration and the inversion result, (b) is the given time-varying flux and the inversion result.

Then observations are increased, and Nos. 1–4 are observed from the second to the seventh day, while other observations are the same as the former. The Strategy (A) is taken as an example, and the inversion result is given in Fig.7. The original observations are marked observed (1), and latter ones are marked observed (2). It is clear that the concentration after assimilation using observed (2) in the last few hours is better than that using observed (1), and the MAE inCC is reduced to 0.32μgL?1. Apparently, the inversion efficiency is raised by increasing observations, which is in accordance with the reality.

5.2 The Inversion Result of Type2

The errors after assimilation in Type 2 are shown in Table 2. The MAEsin CCs with strategies (A)–(F) decline from 21.62μgL?1to 1.45μgL?1or smaller, while the MAEsin pollutant fluxes with these strategies decline from 37.63t to 0.81t or smaller, indicating that the given distribution is successfully inverted. The result with strategy(C) is the best. The errors without IP scheme are a little bigger than those with this strategy, which also shows the inversion result can be improved with IP scheme.

Fig.7 The given distribution and the inversion results.

Table 2 Errors after assimilation

The given distribution and the inversion result with strategy (B) are shown in Fig.8, which shows the given concentration is well inverted with deviation remaining within 2% except those in the last few hours. From Fig.8 it is clear that the adjoint method also makes it possible to invert the given pollutants flux with errors less than 4t.

Fig.8 The given distributions (solid lines) and the corresponding inversion results (dotted lines). (a) is the Given time-varying concentration and the inversion result, (b) is the given time-varying flux and the inversion result.

5.3 The Inversion Result of Type 3

The errors after assimilation in Type 3 are shown in Table 3 from which it is seen that the MAEsinCCs and pollutant fluxes using Strategies (B)–(F) can reach 2.52μgL?1and 1.90tor smaller, respectively.The cost function can even reach the magnitude of 10?5of its initial value, indicating that the given distribution is successfully inverted with all these strategies. However, the MAE in CC and pollutant flux without IP scheme are 3.32μgL?1and 1.79t, respectively.

The inversion results of the given time-varying con-centration and errors between the given time-varying concentration and the inversion results using Strategies (B)–(F) are shown in Fig.9, in which deviation is smaller than 5% except in the last few hours, indicating that the given time-varying concentration is successfully inverted with all these strategies. But the result is not as good as the above inversion result, because the variation is more complex and it is hard to invert. From Fig.9 it also can be seen that the IP strategy has significant influence on inversion results, for inappropriate IP number can cause the result to deviate from the given distribution greatly.

Table 3 Errors after assimilation

Fig.9 The inversion results of the given time-varying concentration (a) and the errors between the given time-varying concentration and the inversion results (b).

The inversion results of the given time-varying pollutant flux and errors between the given time-varying pollutant flux and the inversion results using Strategies (B)–(F) are shown in Fig.10, where deviation is smaller than 3% except in the last few hours, indicating that the given time-varying pollutant flux is also successfully inverted. Asfor the time-varying CC, the inversion result is also a little worse than the above two results.

The MAEsinCCsand pollutant fluxes versus assimilation step using Strategies (B)-(F) are shown in Fig.11. From this figure, we can find the MAEs in CCs decrease rapidly in the beginning, and can decrease to 2.52or smaller after 200 interaction steps. The MAEs in pollutant fluxes are similar to those in CCs. It is evident that the given distributions are successfully inverted using all these strategies.

Fig.10 The inversion results of the given time-varying pollutant flux (a) and the errors between the given time-varying pollutant flux and the inversion results at any moment (b).

Fig.11 TheMAEs in CCs (a) and pollutant fluxes (b) versus assimilation step.

6 Conclusions

The PS information is hard to get due to various reasons, so one of the objectives of this study is to invert the time-varying CCs and then calculate the pollutant fluxes from PS by using the adjoint method. Based on the hydrodynamic background field, anadjoint assimilation modelof pollutantin the Bohai Sea is presented. In order to guarantee the continuity of pollutant distribution over time and make the calculated results more reasonable,theIP scheme is proposed.The CCs atsome time steps are selected as the IPs, and only CCs at these IPs are optimized while thoseat other points are calculated through linear interpolation of the independent CCs. In this way, the number of control variables is reduced.

It is shown that all the given distributions could be estimated by the adjoint method. The cost functions and the MAEs decrease greatly after assimilation, and the cost functions are reduced by about 5 orders of magnitude compared with their initial values. The results indicate that the adjoint method is an effective tool to recover CCs from PS. It is also shown that the IP scheme is significant for the inversionresult, in which appropriate IP number couldyieldbetter inversion results.

The adjoint method can invert the time-varying point source information, which is the first application in oceanography. In the present study only ideal experiments are carried out in which observations are obtained by the model itself, and data in our hand are available to support this study. In the further study we would estimate the source information by assimilating the routine monitoring data in the Bohai Sea, which gives a reference to the scientific countermeasures to control pollution discharge.

Acknowledgements

Partial support for this research was provided by the National Natural Science Foundation of China (Grant Nos. 41072176 and 41371496), the State Ministry of Science and Technology of China (Grant Nos. 2013AA121203 and 2013BAK05B04), and the Fundamental Research Funds for the Central Universities (201262007).

Bousquet, P., Ciais, P., Peylin, P., Ramonet, M., and Monfray, P., 1999. Inverse modeling of annual atmospheric CO2sources and sinks. 1. Method and control inversion., 104 (21): 26161-26178.

Cao, A. Z., Chen, H. B., Zhang, J. C., and Lv, X. Q., 2013. Optimization of open boundary conditions in a 3D internal tidal model with the adjoint method around Hawaii., DOI: 10.1155/2013/950926.

Chen, H. B., Cao, A. Z., Zhang, J. C., Miao, C. B., and Lv, X. Q., 2013. Estimation of spatially varying open boundary conditions for a numerical internal tidal model with adjoint method., DOI: 10.1016/j.matcom.2013.08.005.

Cressman, G. P., 1959. An operational objective analysis system., 87: 367-374.

Da Silva, A. M., Young, C. C., and Levitus, S., 1994.. NOAA Atlas NESDIS 6, US Department of Commerce, NOAA, NESDIS, USA.

Das, S. K., and Lardner, R. W., 1991. On the estimation of parameters of hydraulic models by assimilation of periodic tidal data., 96: 15187-15196.

Elbern, H., Schmidt, H., Talagrand, O., and Ebel, A.,2000. 4D-variational data assimilation with an adjoint air quality model for emission analysis.,15: 539-548.

Fan, W., and Lv, X. Q., 2009.Data assimilation in a simple marine ecosystem model based on spatial biological parameterizations.g, 220(17): 1997-2008.

Gilliland, A., 2001. A sensitivity study of the discrete Kalman filter (DKF) to initial condition discrepancies., 106 (16): 17939-17952.

Gupta, I., Dhage, S., Chandorkar, A.A., and Srivastav, A., 2004. Numerical modeling for Thane creek.,19: 571-579.

Heemink, A. W., and Segers, A. J., 2002. Modeling and prediction of environmental data in space and time using Kalman filtering., 16: 225-240.

Lardner, R. W., and Song, Y., 1995. Optimal estimation of eddy viscosity and friction coefficients for a quasi-three-dimensional numerical tidal model.,33 (3): 581-611.

Lv, X. Q., and Zhang, J. C., 2006. Numerical study on spatially varying bottom friction coefficient of a 2D tidal model with adjoint method., 26: 1905-1923.

Milnes, E.,and Perrochet, P., 2007. Simultaneous identification of a single pollution point-source location and contamination time under known flow field conditions., 30: 2439-2446.

Navon, I. M., Zou, X., Derber, J., and Sela, J., 1991. Variational data assimilation with an adiabatic version of the NMC spectral model., 120: 1433-1446.

Penenko, V. V., and Obraztsov, N. N., 1976. A variational initialization method for the fields of meteorological elements (English translation)., 1: 1-11.

Robertson, L., and Langner, J., 1998. Source function estimate by means of variational data assimilation applied to the ETEX-I tracer experiment., 32 (24): 4219-4225.

Seiler, U., 1993. Estimation of open boundary conditions with the adjoint method., 98: 22855-22870.

Shchepetkin, A. F., and McWilliams, J. C., 2005. The regional oceanic modeling system (ROMS): A split-explicit, free-surface, topography-following-coordinate oceanic model.,9: 347-404.

Shulka, P., 2002. Analytical solutions for steady transport dispersion of nonconservative pollutant with time-dependent periodic waste discharge concentration.,129 (9): 866-869.

Smedstad, O. M., and O’Brien, J. J., 1991. Variational data assimilation and parameter estimation in an equatorial Pacific Ocean model., 26: 179-241.

Wang, C. H., Li, X. Y., and Lv, X. Q., 2013. Numerical study on initial field of pollution in the Bohai Sea with an adjoint method., DOI: 10.1155/2013/104591.

Waterhouse, J., Brodie, J., Lewis, S., and Mitchell, A., 2012. Quantifying the sources of pollutants in the Great Barrier Reef catchmentsand the relative risk to reef ecosystems., 65: 394-406.

Wei, H., Wu, J. P., and Pohlmann, T., 2001. A simulation on the seasonal variation of the circulation and transport in the Bohai Sea., 19 (2): 1-9.

Wei, Z. X., Li, C.Y., Fang, G. H., and Wang, X. Y., 2003. Numerical diagnosticstudyof thesummertimecirculationinthe Bohai Seaandthe water transport inthe Bohai Strait., 21 (4): 454-464 (in Chinese with English abstract).

Wolfe, P., 1969. Convergence conditions for ascent methods., 11: 226-235.

Wolfe, P., 1971. Convergence conditions for ascent methods. II. Some corrections., 13: 185-188.

Zhang, J. C., and Lv, X. Q., 2008. Parameter estimation for a three-dimensional numerical barotropic tidal model with adjoint method., 57 (1): 47-92.

Zheng, B. H., Zhao, X. R., Liu, L. S., Li, Z. C., Lei, K., Zhang, L., Qin, Y. W., Gan, Z., Gao, S. Z., and Jiao, L. X.,2011. Effects of hydrodynamics on the distribution of trace persistent organic pollutants and macrobenthic communities in Bohai Bay., 84: 336-341.

Zhu, J.,and Wang, P., 2006. Ensemble Kalman smoother and ensemble Kalman filter approaches to the joint air quality state and emission estimation problem., 30: 5871-882 (in Chinense).

(Edited by Xie Jun)

DOI 10.1007/s11802-015-2501-8

ISSN 1672-5182, 2015 14 (6): 941-950

? Ocean University of China, Science Press and Springer-Verlag Berlin Heidelberg 2015

(September 30, 2013; revised January 14, 2014; accepted September 2, 2015)

* Corresponding author. Tel: 0086-532-66782971 E-mail: xqinglv@ouc.edu.cn