Abstract
Based on the simulation of a marine ecosystem dynamical model in the Bohai Sea, routine monitoring data are assimilated to study the initial field of pollution by using the adjoint method. In order to reduce variables that need to be optimized and make the simulation results more reasonable, an independent grid is selected every four grids both in longitude and latitude, and only the pollutant concentrations of these independent grids needed to be optimized while the other grids were calculated by interpolation method. Based on this method, the stability and reliability of this model were proved by a set of twin experiments. Therefore, this model could be applied in real experiment to simulate the initial field of the total nitrogen (totalN) in May, 2009. Moreover, the distribution of totalN in any time step could be calculated by this model, and the monthly mean distribution in May in the Bohai Sea could be obtained.
1. Introduction
The Bohai Sea is the only inland sea in China with a maritime area being km2. Its mean depth is 18 m, and the deepest point is located in the west of the Lao Tie Shan channel. Because the runoff of the Yellow River, the Haihe River and the Liao River are all discharged into the Bohai Sea; the organic pollutants are tremendous. Unfortunately, water exchange of the Bohai Sea is very weak and its physical self-cleaning capacity is poor due to its special geographical position. The Bohai Sea is surrounded by land on three sides, and it is only connected with the Yellow Sea through the Bohai strait on the east side. Therefore, it is hard to recover if the Bohai Sea is polluted. A mount of industrial effluent, living sewage, and aquaculture waste water are released into the Bohai Sea with rapid development of the economy along the Bohai Sea, which can cause accumulation of nutrient substances, such as nitrogen and phosphorus. The accumulation of nutrient substance brings a series of ecoenvironmental degradation, including red tide, decrease of seawater oxygen content, decline of biodiversity, and decrease in fish catch. In order to protect and recover the ecoenvironment of the Bohai Sea and to coordinate and improve coastal economy, an accurate simulation of the time-varying pollutant distribution is needed. Only in this way can we achieve rational utilization of marine resources and sustainable development of economy.
Recently, more and more numerical models (e.g., Chen et al. [1], Duan and Nanda [2], Lee and Seo [3], Liu et al. [4], Gupta et al. [5], Periáñez [6, 7], Rajar et al. [8], Rajar and Cetina [9]) for simulating pollutant dispersion have been actually developed since they can be used for decision making after releases of contaminants into the marine environment (Periáñez [10]). Gupta et al. [5] applied a two-dimensional numerical model considering organized wastewater discharges to determine the wastewater assimilative capacity of Thane creek. They found that on the basis of monitoring and simulation, the water quality had been deteriorated significantly due to limited flushing capacity. The volumetric load in the creek needed to be restricted because projected wastewater flows and loads for 2015 were above the assimilative capacity of the creek. Huang et al. [11] investigated the distribution characteristics of heavier or lighter pollutants released at different cross-sectional positions of a wide river with a well-tested three-dimensional numerical model, and their findings assisted in cost-effective countermeasures to be taken for accidental or planned pollutant releases into a wide river. All the three major Siberian rivers, including Ob, Yenisei, and Lena, flow northward into Arctic, and they are supposed to be important sources for various contaminants, so Harms et al. [12] applied a three-dimensional coupled ice-ocean-models of different horizontal resolution to simulate the dispersion of water from these rivers. The model study confirmed that contaminant transport through sediment-laden sea ice offers a short and effective pathway for pollutant transport from Siberian River to the Barents and Nordic Seas. Different methodologies for coupling hydrodynamic submodels with mass transport submodels into integrated water quality models are described in a companion paper (Rajar et al. [8]). The conclusion is that the choice of the methodology depends on the space and time scales, on the prevalent forcing factors and on the nature of the contaminant.
The initial condition has dramatic influence on the simulated results when we use model both in meteorology and oceanography. But in most cases, we do not know the initial field in advance. Therefore, it is important to simulate the initial field accurately. For example, the pollutant distribution characteristics are very different if either the pollutant density or the release location is changed when pollutants are released into a river (Huang et al. [11]). A suite of experimental results of Pe ng and Xie [13] show that although forecast errors due to deficiencies in model physics or numerics cannot always be effectively corrected through improving initial conditions alone, the four-dimensional variational data assimilation algorithm based on Princeton Ocean Model (POM) leads to effective convergence between the forecasts and the ‘‘observations’’ by finding an “optimal” initial condition for the storm surge forecasting. Allen et al. [14] found the combination of source location, source strength, and surface wind direction that best matched the dispersion model output to the receptor data by using genetic algorithm (GA) and demonstrated that the GA was capable of computing the correct solution as long as the magnitude of the noise did not exceed that of the receptor data.
In order to obtain the initial condition or the average distribution of pollutant concentration in a certain period, the traditional way is to use all the observations within this period by interpolation method, such as Cressman and Kriging. However, this method only takes the spatial information of the observations into consideration but ignores the time information, which may lead to a big difference between the simulation results and the actual situation. Adjoint method minimizes a predetermined cost function which defines differences between model-derived quantities and measured quantities by using the existing observations to the maximum extent. Adjoint method is an effective variational assimilation technique based on mathematics strictly. It takes ocean model equations, initial conditions, and boundary conditions as the constraint conditions and combines variational principle and the theory of optimal control. With the accessible oceanic elements, some inaccessible oceanic elements can be obtained by optimizing initial field and/or parameters (Sasaki [15]; Lu and Zhang [16]).
Adjoint method was widely used in both meteorology and oceanography. Zhang et al. [17] applied this method to study the similarities and the differences between the Ekman (linear) and the Quadratic (nonlinear) bottom friction parameterizations for a two-dimensional tidal model. The simulation results indicated that the nonlinear Quadratic parameterization is more accurate than the linear Ekman parameterization if the traditional constant boundary friction coefficient is used. However, when the spatially varying boundary friction coefficients were used, the differences between the Ekman and the Quadratic approaches diminished. In the study of Fan and Lv [18], SeaWiFS chlorophyll-a data were assimilated into a simple NPZD model by the adjoint method in a climatological physical environment provided by FOAM. The results showed that the values of the selected sensitive parameters were spatially variable and the application of spatial parameterizations could improve the assimilation results significantly. Many researches (Yu and O’Brien [19], Lawson et al. [20], Zhao et al. [21], Zhao and Lu [22], Qi et al. [23]) have proved the validity and rationality of the adjoint method. Therefore, in this paper, we apply this method to simulate the distribution of totalN in any time step, and the monthly mean distribution in May in the Bohai Sea can be obtained.
The contents of this paper are organized as follows. Section 2 describes the ecosystem model and the database. Section 3 illustrates the adjoint method and independent grids briefly. Section 4 describes the twin experiment to validate the model’s capability of inversing pollutant initial field and finds the optimal strategy of setting independent grids. Based on twin experiments, practical experiment is performed in Section 5 in order to obtain the initial field of totalN in the Bohai Sea in May. The conclusions of our work are presented in Section 6.
2. Model and Data
2.1. Model Equations
Based on hydrodynamic model, the transporting diffusion process of pollution can be written as follows: where and represent horizontal and vertical eddy diffusivities, respectively, is the concentration of pollution, and is the degradation coefficient of pollution. When the pollution is conservative substance, ; otherwise, . In this paper, we treat the pollution as conservative substance, so ; the finite-difference form can be seen in Appendix A.
The three-dimensional Regional Ocean Model System (ROMS) is used to calculate the ambient physical velocities and temperature in the Bohai Sea (37°N~41°N, 122.5°E~127.5°E, Figure 1), and the horizontal resolution is 4 second in both latitude and longitude. The thickness of each layer from top to bottom is 10 m, 10 m, 10 m, 20 m, 25 m, and 25 m, respectively, and the integral time step is 6 hours.
Monthly mean horizontal currents in May in the Bohai Sea obtained by ROMS in 5 m depth, 15 m depth, and 25 m depth are shown in Figures 2(a), 2(b), and 2(c), respectively. On the whole, mean flow velocity decreases gradually from surface to bottom in the Bohai Sea. There is a clockwise vortex in the central Bohai Sea, and an anticlockwise vortex can be seen in the east part of Liaodong Bay. The circulation in Laizhou Bay is very weak, and there is a clockwise vortex in the mouth. Meanwhile, water flows out of Bohai Sea through north of Bohai Strait and into Bohai Sea through south of Bohai Strait.
(a)
(b)
(c)
2.2. Observations and Model-Generated Observations
The distribution of the conventional monitoring stations is depicted in Figure 3. We can see that most of the stations are located in the Bohai strait and the coastal areas while only a few of them are located in the central Bohai Sea. The observations are needed in practical experiments and some model-generated observations are needed in the twin experiments. We can choose some model-generated observations through the following methods: first, an initial pollutant distribution (initial field) in the Bohai Sea is assigned. Then, the forward model is run for 30 days, so the pollutant distribution at every time step can be obtained. As the sampling locations of the conventional monitoring stations have been known from Figure 3, we can pick up model-generated observations according to the following method: the sampling locations of the model-generated observations are the same as the conventional monitoring stations. The total number of the conventional monitoring stations is 121. Since the total calculating step is also 121, we prescribe that the number of model-generated stations are in one-to-one correspondence with the sampling time for the sake of simplicity. The observation numbers of each station depend on the depth of water. If the water depth is within three layers, the pollutant concentration of every layer is chosen as model-generated observations; otherwise, only the pollutant concentration of upper three layers is chosen as model-generated observations.
If a guess initial field of pollution is given, the initial field of the pollution can be optimized by using the observations (practical experiment) or the model-generated observations (twin experiment) through the adjoint method.
3. Method
3.1. Adjoint Assimilation Method
Adjoint assimilation method treats all the practical problems as minimum problems. It takes model equations, initial conditions, and boundary conditions as constraint conditions and minimizes a predetermined cost function which defines differences between model-derived quantities and measured quantities. The flowchart in Figure 4 summarizes the steps that make up the adjoint method.
Step 1. an initial distribution of pollution is given empirically (guess distribution);
Step 2. perform the simulation by running the forward model, and the simulation results are obtained;
Step 3. Calculate the cost function which defines the misfit between simulation results and observations. The cost function is defined by where represents the simulation result, is the observation, the index triplet is a pointer to certain grid cell, and is the weight of the observation, which is defined as follows:
Step 4. the adjoint of the model (Appendix B) is run backward in time;
Step 5. calculate the gradient of the cost function with respect to initial field;
Step 6. update the initial pollutant distribution closer to the minimum of the cost function;
Step 7. return to Step 1; repeat the iteration with the updated pollutant distribution;
Step 8. end this procedure after a specific number of iterations or the cost function is small enough to meet the criteria ( is a small real number, such as 0.01), and the optimized initial field of pollution is obtained. In this paper, we choose the former for easiness to compare the simulation results in twin experiment 1.
3.2. Independent Grids
If the pollutant concentration of each grid is optimized independently, then there are too many variables to constraint, and the pollutant distribution is not continuous, which is not reasonable. So, we can reduce variables that need to be optimized and guarantee that the simulation result coincides with the law of physics by using independent grids. Several grids are selected as independent grids and only pollutant concentrations of these independent grids need to be optimized while those of other grids are calculated by Cressman method [18]. Since cost function declines in the inverse direction of its gradient, the gradient is used to determine the direction to optimize the pollutant concentration. In this paper, the distributions of independent grids in longitude are the same as in latitude.
4. Numerical Experiments and Result Analysis
4.1. Twin Experiment 1: The Strategy of Independent Grids
The simulated results are affected by the number of independent grids, so we will discuss this factor in twin experiment 1. The experiment is designed as follows. Assume that the initial distribution of pollutant concentration shows a parabolic surface with upward convex, which means it is high in the centre and low in the surroundings. The independent grids are selected every 2 to 9 common grids. The influence radius is 1.2 times of the distance between adjacent independent grids, and the number of iterations is set to 100 mainly based on the following considerations: (1) the misfit between “observations” and the simulated results is very small and approximately constant when it declines to a certain value after 100 iterations; (2) the cost functions are almost no longer falling after 100 iterations; and (3) the calculation amount is acceptable for the 100 iterations. The results are given in Figure 5.
In Figure 5, X-axis represents that there is one independent grid every 2 to 9 common grids. The solid line shows the ratio of cost function to its initial value while the dotted line indicates mean absolute error (MAE) of the model-generated observations. It shows that when independent grids get fewer, both the cost function and the MAE of the model-generated observations get smaller at first but larger after 1.2°. This means that choosing an independent grid every four grids is the best choice. Therefore, in the following experiments, we select an independent grid every four grids both in longitude and latitude.
4.2. Twin Experiment 2: The Pollutant Concentration Shows a Parabolic Surface
4.2.1. Parabolic Surface with Upward Convex
Assume that the initial distribution of pollutant concentration shows a parabolic surface with upward convex, which means it is high in the centre and low in the surroundings, and the pollutant concentration at any grid can be calculated by where and indicate the longitude and latitude of grid , respectively. Equation (4) shows that the pollutant concentration varies from 0.05 to 2.10.
We guess that the initial pollutant concentration at any grid is 0.05 mg/L, and the number of iterations is the same as twin experiment 1.
As we can see from Figures 6 and 7, the cost function can reduce to 3.7 percent of its initial value. Table 1 shows that the MAE of the model-generated observations declines from 0.53 mg/L to 0.06 mg/L, which decreases by 88.7 percent, and Table 2 shows that the mean relative error (MAE) also declined obviously after adjoint assimilation. The results indicate that through adjoint assimilation, the model-generated observations have been effectively used, and the given distribution can be inverted successfully. The given distribution and inversion results can be seen in Figures 8 and 9. The inversion results of the central Bohai Sea and the Laizhou Bay are very satisfactory while those of the transition zones between the central Bohai Sea and the three bays are a little worse. The inversion results are basically the same as the prescribed distribution.
4.2.2. Parabolic Surface with Downward Convex
Assume that the initial distribution of pollutant concentration shows a parabolic surface with downward convex, which means it is low in the centre and high in the surroundings, and the pollutant concentration at any grid can be calculated by
Repeat the process of Section 4.2.1. The cost function can reduce to 8.3 percent of its initial value, and the MAE of the model-generated observations declines from 1.5 mg/L to 0.2 mg/L, indicating that the model is stable and valid. No matter which convex is given, the initial distribution can be inverted successfully.
4.3. Twin Experiment 3: The Pollutant Concentration Shows a Conical Surface
Assume that the initial pollutant concentration shows a conical surface with upward convex, and the pollutant concentration at any grid can be calculated by
Assume that the initial pollutant concentration shows a conical surface with downward convex and the pollutant concentration at any grid can be calculated by
Twin experiments can also be evaluated by the absolute and the relative difference between the prescribed and the inversion results.
In these circumstances, the initial pollutant concentration also varies from 0.05 to 2.10. Repeat the process of twin experiment 2, and the results are similar to those of twin experiment 2. The detailed information is presented in Table 1. These results once again verify the stability and validation of this model, indicating that the coupled model can be used in practical experiment. In other words, the initial field of the Bohai Sea can be inverted through adjoint method by using the regular monitoring observations.
5. Practical Experiment and Result Analysis
The distribution of totalN monitoring stations in May, 2009 is shown in Figure 10. The larger the dot is, the greater the observation it represents. According to the traditional method, the monthly mean distribution of totalN in surface layer is always obtained by interpolation method, such as Cressman and Kriging. The monthly mean distribution of totalN in surface layer calculated by Cressman method is depicted in Figure 11. When we make use of the coupled model mentioned in this paper, not only the initial totalN distribution in May, but also the distribution in any time step can be obtained. And then, we can get the monthly mean distribution of totalN in surface layer, which is depicted in Figure 12.
In the Bohai strait and the central Bohai Sea, the distribution of totalN calculated by the Cressman method is similar to the result obtained by the adjoint method. However, the latter is higher in the Laizhou Bay and lower in the Liaodong Bay than the former. The reason may be that the sampling time of each routine monitor station is not the same, and the observation cannot always reflect the mean value in a month. Moreover, the Cressman method only takes the spatial information of the stations into consideration and ignores the time information.
The average concentration of totalN in surface layer in the Bohai Sea is 0.49 mg/L by using the Cressman method while it is only 0.41 mg/L by using the adjoint method. The former is almost twenty percent greater than the latter. The reason is probably that most of routine monitoring stations are located in the Bohai strait and the coastal areas while only a few of those located in the central Bohai Sea. In the meantime, the observations in the Bohai strait and the coastal areas are of relativly high values while the observations in the central Bohai Sea are of relative low values. When we calculate the pollutant concentration between high-value observations and low-value observations by using the Cressman method, there may be more high-value observations involved in the calculation. Therefore, it will lead to a larger pollutant concentration than the real value in most areas. Moreover, although the interpolation results can guarantee the continuity of the totalN distribution, the totalN concentration at any grid cannot be larger or smaller than the observations we already have.
The adjoint method can make use of the observations to a maximum extent. By using this method, not only the values and locations of the observations, but also the sampling time is taken into consideration. The distribution characteristics of totalN in the Bohai Sea in May are depicted in Figure 12. The concentration of totalN is lowest in the central Bohai Sea and Bohai strait. In fact, the concentration is almost zero in these areas, and the highest concentration is no more than 0.5 mg/L. The totalN concentrations of three bays are higher than the central Bohai Sea, among them the concentration of Bohai Bay is the lowest, while Liaodong Bay takes the second place, and the concentration of most areas is less than 1.0 mg/L. The totalN concentration of the Laizhou is the highest. It is more than 1.0 mg in half of areas, and it even reaches 1.4 mg/L in some areas.
6. Conclusion
Recently, more and more numerical models for simulating pollutant dispersion have been actually developed since they can be used for decision-making purposes after releases of contaminants into the marine environment (Periáñez [10]). Regardless of which numerical model we employ, the initial field of pollution has dramatic influence on the simulated results. Therefore, in order to simulate the dispersion process of the pollution accurately, we need not only a good pollutant model and reasonable model parameters, but also an accurate initial distribution of pollution. Although some researchers have tried to inverse the pollutant location and strength, this is the first time to inverse the initial distribution of pollution as far as we know.
Based on the simulation of a marine ecosystem dynamical model in the Bohai Sea, routine monitoring data were assimilated to study the initial field of pollution in this paper. Firstly, in order to reduce variables that need to be optimized and make the simulation results more reasonable, an independent grid is selected every four grids both in longitude and latitude, and only the pollutant concentrations of these independent grids needed to be optimized while the other grids were calculated by interpolation method. Based on this method, the stability and reliability of this model were proved by a set of twin experiments. Therefore, this model could be applied in real experiment to simulate the initial field of the totalN in May, 2009. Furthermore, the distribution of totalN in any time step could be calculated by this model, and the monthly mean distribution in May in the Bohai Sea could be obtained. Compared with the Cressman method, the adjoint method can make use of the existing observations to the maximum extent. Not only the values and locations of the observations, but also the sampling time are taken into consideration by using this method. Therefore, it can reduce the misfit between inversion results and observations significantly and make the simulated results closer to reality. This method can also be used to simulate the distribution of other pollutions, such as total phosphorus, chemical oxygen demand, and petroleum hydrocarbon. So, it could most probably be used in solving environmental problems in the future.
Appendices
A. Finite-Difference Form of Equation (1)
Consider the following: where and are horizontal velocities and is vertical velocity.
B. Adjoint Equation
Consider the following: where , represent observation and corresponding simulating result, respectively. is the adjoint operator of . The finite-difference form of adjoint equation is as follows:
Acknowledgments
The authors acknowledge the support of the Major State Basic Research Development Program of China through Grant 2013CB956500, the National Natural Science Foundation of China through Grants 41076006 and 41206001, the Natural Science Foundation of Jiangsu Province of China through Grant BK2012315, and the Fundamental Research Funds for the Central Universities Grants 201261006 and 201262007.