PARAMETER IDENTIFICATION AND NUMERICAL SIMULATION FOR THE EXCHANGE COEFFICIENT OF DISSOLVED OXYGEN CONCENTRATION UNDER ICE IN A BOREAL LAKE

. Dissolved oxygen (DO) is one of the main parameters to assess the quality of lake water. This study is intended to construct a parabolic distributed parameter system to describe the variation of DO under the ice, and identify the vertical exchange coeﬃcient K of DO with the ﬁeld data. Based on the existence and uniqueness of the weak solution of this system, the ﬁxed solution problem of the parabolic equation is transformed into a parameter identiﬁcation model, which takes K as the identiﬁcation parameter, and the deviation of the simulated and measured DO as the performance index. We prove the existence of the optimal parameter of the identiﬁcation model, and derive the ﬁrst order optimality conditions. Finally, we construct the optimization algorithm, and have carried out numerical simulation. According to the measured DO data in Lake Valkea-Kotinen (Finland), it can be found that the orders of magnitude of the coeﬃcient K varying from 10 − 6 to 10 − 1 m 2 s − 1 , the calculated and measured DO values are in good agreement. Within this range of K values, the overall trends are very similar. In order to get better ﬁtting, the formula needs to be adjusted based on microbial and chemical consumption rates of DO.

1. Introduction. Dissolved oxygen (DO) is a useful ecological index of lake water, especially for the frozen lakes. All multicellular organisms in a lake, from fish to insects and to microscopic zooplankton, need oxygen for respiration [1]. While ice free, photosynthesis and atmospheric exchange, largely via wind aeration, are the primary inputs of oxygen into lakes and losses are due to organismal respiration plus chemical oxidation of reduced materials [15]. When lake surface is frozen, the water body is isolated from the atmosphere and oxygen supply from there is blocked. Oxygen-consuming processes in the lake rely predominantly on the amount of oxygen dissolved in water before freezing. Once the oxygen supplies produced by water aeration and photosynthesis in the water body are limited, oxygen depletion may occur under the ice. As a result anoxic conditions may be generated with several ecological consequences, such as fish kills, loss of benthic organisms, and dominance of anaerobic processes with the accumulation of hydrogen sulfide and other harmful compounds [4,6,13,19]. Fish kill can occur in shallow lakes near the end of the ice cover period when DO concentrations drop below 2 mg L −1 at all depths over several days [8,18].
Many researchers have been devoted to the DO study in lakes and different models have been developed and used to analyze the change of DO, mainly including deterministic and stochastic/statistical models. And computational-intelligence techniques such as neural networks, fuzzy logic and combined neuro-fuzzy systems have become very effective tools to identify and model the DO system [7,12,17]. Zhang et al (2015) [23] examined the relationships between oxycline depth and thermocline depth, daily average air temperature and Secchi disk depth (SDD), and further assessed the response of DO stratification to climate warming and lake eutrophication. Patterson et al (1985) [16] extended a vertical mixing model of Lake Erie with an oxygen budget model to simulate vertical distributions of temperature and dissolved oxygen over relatively short time periods with data gathered in the summers of 1979-80. Stefan and Fang (1994) [20] formulated a one-dimensional, unsteady DO model to simulate summer DO conditions in a wide range of lakes of the north central United States and to study potential impacts of global climate change. Mathias and Barica (1980) [14] analyzed the winter oxygen depletion rates from four sets of Canadian lakes and found that there was an inverse relationship between oxygen depletion rate and mean depth. Babin and Prepas (1985) [4] determined winter oxygen depletion rates for thirteen ice-covered lakes in central Alberta in Canada , and predicted DO depletion rates by a combination of mean summer total phosphorus in the euphotic zone and mean depth. Foley et al (2012) [9] analyzed the consequences of eutrophication and climate warming on hypolimnetic oxygen depletion based on data of years 1968-2008 in a small lake in Blelham Tarn, UK. Fang and Stefan (1997) [8] simulated climate change effects on DO characteristics in ice-covered lakes in the temperate zone (in Minnesota, USA) with a year-round DO model. Antonopoulos and Gianniou (2003) [2] calibrated and verified a DO model by using data of Lake Vegoritis in 1981 and 1993, which was based on the unsteady diffusion equation with oxygen fluxes through the ice-free surface and bottom of the lake as boundary conditions, photosynthetic oxygen production, biochemical oxygen demand and plant respiration as internal sources and sinks. Golosov et al (2007) [10] prepared a one-dimensional nonstationary equation for the vertical transfer of the non-conservative substance and the expression for the vertical DO flux was written as Where C is the DO concentration (unit: mg L −1 ); z is the depth (unit: m); t is the time (unit: s); Q is the vertical flux of DO (unit: mg L −1 ·m s −1 ); T (z, t) is the vertical distribution of water temperature (unit: Kelvin); γ[T (z, t)] is the total rate of DO consumption (unit: s −1 ) ; and K is the coefficient for the vertical mass exchange. As a rule, there is no preliminary reliable information on the variability in the value of K, especially concerning ice-covered lakes [10]. Hence, in this paper, we aim at assessing the change of K with the field data measured in an ice-covered lake in Finland utilizing the parabolic distributed parameter system. In this paper , we present a parabolic distributed parameter system to describe the variation of DO under ice, and prove the existence and uniqueness of the weak solution of the evolution equation on basis of some important properties of a bilinear functional in Section 2. In addition, we transform the fixed solution problem into a parameter identification model, and prove the existence of the optimal parameter of the identification model, and derive the first order optimality conditions in Section 3. A numerical optimization algorithm is constructed and applied to numerical simulation in Section 4, and numerical results are finally presented in Section 5.
2. DO model under ice with variable boundary.

2.1.
The form of DO model system. The ice-water interface is taken as the coordinate origin O, and the direction of the vertical depth along the origin is positive z axis, as shown in Fig. 1; Let z u , z l be upper and lower boundaries of water depth respectively (unit: m ) and z u < z l , which z l = z(t) is a depth function varying with the time t (unit: s). C(z, t) represents the DO concentration at the depth z and the time t (unit: mg L −1 ), which C 0 (z) is the DO concentration at the initial time, and C u (t) and C l (z(t), t) are the DO concentrations of the upper and lower boundary at the depth z u and z l = z(t). T (z, t) (unit: Kelvin) is water temperature at the depth z and the time t. Let the space domain of the system Ω = (z u , z l ), and the time domain I = (0, t f ). And we denote that Ω = [z u , z l ] , I = [0, t f ] . According to the above model (1)-(2), the DO model with the initial value condition and variable boundary can be described by the following parabolic distributed parameter system where K is the vertical exchange coefficient, and γ[T (z, t)] is DO consumption rate. As the measuring interval is the water layer under ice, the expression of the DO consumption rate is as follows [10]: In the formula (4), T min and T max are the minimal and maximal possible values of the water temperature during the ice period. Usually the vertical distribution of the consumption rate is characterized by minimal values (γ min ) in the upper layers of the water column and maximal values (γ max ) in the upper layer of sediments. These parameters should be determined from observational data.
In the analysis of the model, according to the actual situation, we give the following assumptions: Assumption 1. In practical problems, the coefficient for the vertical mass exchange K is bounded in the discussion area, and there exist K l > 0 and K u > 0 satisfying K l ≤ K u , which makes K l ≤ K ≤ K u .
then the system (3) can be described as the following homogeneous boundary value problem where 2.2. Properties of the DO system. Let U ad (I) = K(t)| K(·) : I → Γ is measuable, K l ≤ K ≤ K u } is the admissible parameter set, K ∈ U ad (I), and Γ is a known nonempty, bounded and closed convex set on R. And let Ω K (t) be the spatial domain of the homogeneous boundary value system (5) at the time t. For any t ∈ I and K ∈ U ad (I), let H K (t) = L 2 (Ω K (t); R) be a Hilbert space, and · V K (t) respectively. Let H K (t) and V K (t) are the dual spaces of H K (t) and V K (t) respectively, and ·, · V K (t), V K (t) the dual paring between V K (t) and V K (t), The inner product and the norm on H K (t) are expressed as And the inner product and the norm on V K (t) are expressed as Then the dual paring ·, We multiply both sides of the first equation of system (5) by any ξ in V K (t), and integrate them over Ω K (t) then By using the partial integral formula and the homogeneous boundary conditions, we can get then the former (13) can be described as For any t ∈ I, the functional δ(t, K; φ, ϕ) on V K (t) × V K (t) can be defined as Three properties of the functional δ(t, K; φ, ϕ) are presented in the following.
is a continuous function on the measurable set I , then it is a measurable function on I. ( Proof. The conclusion can follow directly from the definition of δ(t, K; φ, φ) and Gårding inequality [22].
From Properties 1 and 2 and Lax-Milgram theorem [22], there exists a bounded continuous linear operator A(t, K) ∈ L(V K (t), V K (t)) such that and According to (16) and Properties 1 and 2, we can get the following property.
Property 3. If Assumptions 1-3 hold, then for all t ∈ I and for all φ, ϕ ∈ V K (t) and K ∈ U ad (I), there exist M 1 > 0 and M 2 > 0 such that ∀t ∈ I and K ∈ U ad (I), let P (t, K; φ) = Ω K (t) φp(z, t)dz , which φ ∈ V K (t) and p(z, t) is defined as the equation (6), then P (t, K; φ) is a continuous linear functional on space V K (t). According to Riesz representation theorem [22], we know that there exists a unique element q(t; K) ∈ V K (t) such that From the above analysis, we can know that the system (5) can be described as the following parabolic evolution equation The above reasoning process is a reversible process, so the system (5) is equivalent to the development equation (18).
Next, we will discuss the problem of weak solutions of the development equation (18). Firstly, the definition of weak solutions of (18) is given below.
is a weak solution of the evolution equation (18), if that for ξ ∈ V K (t) and K ∈ U ad (·).
According to the definition of V K (t) , for any given t ∈ I, every ϕ ∈ V K (t) satisfies the Dirichlet condition ϕ(z) = 0, z ∈ ∂Ω K (t)\∂Ω 0 K (t). Then we can get that the evolution equation (18) has a unique weak solution v(z, t; K), which depends continuously on q and v 0 .
3. Parameter identification of the DO system.

3.1.
Definition of identification model. The performance index of the identification model is defined as here v(z, t; which C(z, t; K) is DO concentration function of the system (3), and C(z, t) is measured DO concentration. Let the equation (19) can be expressed as In order to make DO concentrations obtained from the system (3) fit the measured data as far as possible, we can minimize the error of DO concentration between the measured and calculated by system (3) and the following parameter identification model is established where S U ad (Ω K (t)×I) ⊂ L 2 (Ω K (t), I; H K (·))∩C Ω K (t), I; H K (·) is the weak solution set of the evolution equation (18), i.e. S U ad (Ω K (t)×I) = {v(z, t; K) |v(z, t; K) is the solution of system (18) corresponding to K ∈ U ad (I)} and J(K) is defined by the expression (22). By the expressions (19)- (21), identification model (23) can be converted to the following questions where S U ad (Ω K (t) × I) = {C(z, t; K) |C(z, t; K) is the solution of system (3) corresponding to K ∈ U ad (I)}.

3.2.
Existence of optimal parameter and necessary optimality condition.
Lemma 3.1. If Assumptions 1-3 hold, then for all K ∈ U ad (I) , the mapping K → v(z, t; K) is strongly continuous.
Proof. For given K 0 ∈ U ad (I) , let feasible parameter sequence {K n , n = 1, 2, · · ·} in U ad (I) such that K n − K 0 → 0 as n → ∞ , and v n and v 0 are the weak solutions of the evolution equation (18) corresponding to K n and K 0 respectively. Set w n = v n − v 0 , p n = p(z, t; K n ), p 0 = p(z, t; K 0 ), we can get ∂w n ∂t − K ∂ 2 w n ∂z 2 = p n − p 0 , w 0 (z, 0) = 0, w n | z=zu = 0, w n z=z(t) = 0.
and for the second term of the equation (26) on the left side by Property 2, we have t 0 δ(s, K n ; w n , w n )ds + µ here µ ≥ 0 and η > 0 are the same as Property 2.
According to the Cauchy inequality ab ≤ ε 2 a 2 + 1 2ε b 2 , (a > 0, b > 0, ε > 0), the right term of (26) is expressed as Substituting (27)-(29) into (26), then we have By the third and fourth equations (boundary conditions) of the system (5) as well as by Gronwall's inequality, we can get Combined with K n − K 0 → 0, we can get |p n − p 0 | → 0 and so the mapping K → v(z, t; K) is strongly continuous.
Theorem 3.2. If Assumptions 1-3 hold, then for all K ∈ U ad (I), there exists at least one optimal parameter K opt ∈ U ad (I) satisfies that J(K opt ) ≤ J(K), i.e. K opt is the optimal parameter of the identification problem (24).
According to Lemma 3.1, the mapping K → v(z, t; K) is strong continuous on U ad (I) , then K → J(K) is a continuous mapping on U ad (I) .
And Since U ad (I) is a nonempty, bounded and closed set, for all K ∈ U ad (I), there is K opt ∈ U ad (I) such that J(K opt ) ≤ J(K).
The necessary optimality condition of the identification system will be given in the following. We know that U ad (I) is a compact and convex set, and K → J(K) is a continuous mapping on the set U ad (I) . If K opt ∈ U ad (I) is the optimal parameter of the identification problem (24), we can J(K(·)) is Gateaux differentiable at K opt and its Gateaux derivative DJ(K opt (·)) exists. Theorem 3.3. If Assumptions 1-3 hold, K opt ∈ U ad (I) is the optimal parameter of the identification problem (24), then for all K ∈ U ad (I), that K opt (·) satisfies the following inequality 4. Optimization and numerical algorithm.

4.1.
Parameter identification model. In this part, we will construct an optimization algorithm to solve the parameter identification problem (24). The solution domain of system (3) is Ω×I = {(z, t) |z ∈ [z u , z l ] , t ∈ [0, t f ]}, which z l = z(t) is a depth function varying with the measurement time t . There are N z observation depths, i.e. z k ∈ Ω, k ∈ {k 1 , k 2 , k 3 , · · · , k Nz } in the spatial region Ω and N t observation times, i.e. t j ∈ I, j ∈ {j 1 , j 2 , j 3 , · · · , j Nt } in the time region. For a given parameter K ∈ U ad , the C (z k , t j ; K) is DO concentration calculated by the system (3) at the depth z k and the time t j and C(z k , t j ) is measured DO concentration. Then the identification problem (24) can be expressed as where is the solution of system (3) corresponding to K ∈ U ad (I)} In conclusion, according to the theory of partial differential equations, we can have the following theorem.
Theorem 4.1. If Assumptions 1-3 hold, then the system (3) has a unique weak solution C(z, t; K) ∈ C (Ω K (t), I, U ad ; R) which depends continuously on K ∈ U ad (I).

Numerical optimization algorithm.
It is difficult to obtain the analytical solution about the parameter system (3). The semi-implicit difference method is an unconditionally stable finite difference method, and it does not need to solve algebraic equations. In this section, we apply the semi-implicit difference scheme to solve the system (3) and construct the optimization algorithm to solve identification problems (32).
In the plane, we take z as the spatial direction and t as the time direction. And ∆z = h is the spatial step and ∆t = τ is the time step, n z = (z(t) − z u )/h , n t = t f /τ , where · denotes the largest integer less than or equal to the element.
Set C j k = C(z k , t j ), z k−1 ≤ z k ≤ z k+1 , t j+ 1 2 = t j + τ 2 , the first type of system (3) is discretized at the grid node (k, j + 1 2 ), and we approximate the time derivative with central difference ∂C ∂t The spatial partial derivatives are approximated by forward and backward difference schemes. The forward difference scheme is and the backward difference scheme is The formulas (33) and (34) are brought into the first equation of the system (3), and the truncation error The formula of the forward difference scheme is τ h 2 , j = 0, 2, 4, 6, · · · ; k = 1, 2, 3, · · · , n z − 1.
Similarly, we can get the formula of the backward difference scheme is In the process of solving equation (32), it is different in the order of calculation and the iterative equation. When j is even (j = 0, 2, 4, 6, · · · ), the system is calculated by the forward difference method (from top to bottom in Fig. 1), i.e. which is calculated by the formula (36). But when j is odd (j = 1, 3, 5, · · · ), the system is calculated by the backward difference method (from bottom to top in Fig. 1), i.e. which is calculated by the formula (37). The C j+1 k−1 and C j+1 k+1 in formulas (36) and (37) are given through the boundary value or the previous iterative calculation results, which avoids solving algebraic equations, so it has advantage of the explicit difference scheme.
The specific steps of the algorithm can be described as follows.
Step 2. If j is even, then go to Step 3; else, go to Step 4; Step 3. When the time node j = 1, 3, 5, · · · is odd, DO concentrations C j+1 k in each space node is calculated by the formula (36), and DO concentrations in the corresponding position with the measured depths are kept, and set j = j + 1. If j = n t , then go to Step 5; else, go to Step 4.
Step 4. When the time node j = 2, 4, 6, · · · is even, DO concentrations C j+1 k in each space node is calculated by the formula (37), and DO concentrations in the corresponding position with the measured depths are kept, and set j = j + 1. If j = n t , then go to Step 5; else, go to Step 3.
Step 5. Calculate the objective function of the model (32). If i > 0, go to Step 6 , else set i = i + 1, go to Step 2.
Step 6. If J d (K i ) < J d (K i+1 ), set K opt = K i , else set K opt = K i+1 . If i = n, stop, else, go to Step 2.
5. Numerical results. DO concentration data came from Lake Valkea-Kotinen in southern Finland [3,11,21] Fig. 2. The zero reference level of the depth is the ice-water interface. The lake was frozen on November 8, 2010 and the ice broke up on April 26, 2011 [5].
In the process of simulation identification, we compare the measured and calculated DO concentrations at different depths when the vertical exchange coefficient K varies from 10 −8 to 10 m 2 s −1 with a total of ten orders of magnitude. The The relative error (Relerr) is expressed as and the correlation coefficient (Corcoef) is denoted as  Table 1, and correlation coefficients are shown in Table 2. Table 1. Relative errors (%) of the measured and the calculated DO concentrations at different depths when K varying from 10 −8 to 10 m 2 s −1 with ten orders of magnitude   Table 2. Correlation coefficients of the measured and the calculated DO concentrations at different depths when K varying from 10 −8 to 10 m 2 s −1 with ten orders of magnitude  In the following, results of the depth of 0.70 m were taken as an example to discuss the calculated and the measured values with the variation of K, and the contrast images are shown in Fig. 3. For K changes in the two magnitudes, 10 −8 and 10 −7 m 2 s −1 , trend gap is relatively large between the calculated and the measured values judging from the image. And K changes in the other two orders of magnitude, 10 0 and 10 1 m 2 s −1 , observing from the image, the volatility of the calculated values is larger than the measured values. When the K value changes from 10 −6 to 10 −1 m 2 s −1 , from the image, the simulated and the measured values are in good agreement with the trend change, and as the K value becomes larger gradually, the fluctuation of the trend becomes larger. The smaller the K value, the smaller the fluctuation of the simulated data. At the same time, it can be seen from Table 1, the relative errors between the simulated and measured values at the depth of 0.70 m change in 0.16%-0.18%, so the gap is very small. In Table 2, the correlation coefficients vary in the range of 0.96-0.97. In the process of the whole simulation, we observe that the simulation of the others besides 1.95 m is relatively good. It is the worst in simulation result at the depth of 1.95 m, and relative errors even reach 30%, but relative errors have no obvious change in the process of the change of K value; and the correlation coefficients change between 0.96 and 0.97, which shows that the trend of simulation is consistent with the measured values. The above analysis indicated that the value of K varied between 10 −6 to 10 −1 m 2 s −1 . The simulated DO concentrations followed well with the measured ones and both had almost equal trend.
And for the final fitting effect, when K equals 10 −4 m 2 s −1 order of magnitude, the comparison curves between measured and calculated DO values at different depths are given in Fig. 4. 6. Conclusions. According to the variation of DO in ice-covered lakes, we constructed a parabolic distributed parameter system to describe the change of DO under ice, and proved the existence and uniqueness of the weak solution of this system. The fixed solution problem of the parabolic equation was transformed into a parameter identification model. The model took the vertical exchange coefficient K of the DO as the identification parameter, and the deviation between simulated and measured DO as the performance index. We proved the existence of the optimal parameter of the identification model, and derived the first order optimality conditions. Finally, we constructed the optimization algorithm, and carried out numerical simulation to discuss the range of the magnitude of the coefficient K based on the measured DO data and water temperature in Lake Valkea-Kotinen.
There is no reliable information available on the K about DO under ice in lakes. According to the data of Lake Valkea-Kotinen, the magnitude of Kvaried from 10 −6 to 10 −1 m 2 s −1 . Within this range of K values, the overall trends were very similar, and the stability of the data was slightly different.
We can conclude that if we want to get better fitting effect in Lake Valkea-Kotinen, any modification of the K value cannot essentially improve the model fitting. We suggest the reason is that the DO consumption rates at different depths are strongly dependent on the DO consumption caused by microbes and chemical processes in the sediment and water column. And as a result in a small lake under limited exchange strength beneath the ice, and the vertical exchange coefficient K does not have a major role in determining DO concentrations at different depths.