Simulation and Domain Identification of Sea Ice Thermodynamic System

Copyright q 2012 Bing Tan et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. Based on the measured data and characteristics of sea ice temperature distribution in space and time, this study is intended to consider a parabolic partial differential equation of the thermodynamic field of sea ice coupled by snow, ice, and sea water layers with a time-dependent domain and its parameter identification problem. An optimal model with state constraints is presented with the thicknesses of snow and sea ice as parametric variables and the deviation between the calculated and measured sea ice temperatures as the performance criterion. The unique existence of the weak solution of the thermodynamic system is proved. The properties of the identification problem and the existence of the optimal parameter are discussed, and the one-order necessary condition is derived. Finally, based on the nonoverlapping domain decomposition method and semi-implicit difference scheme, an optimization algorithm is proposed for the numerical simulation. Results show that the simulated temperature of sea ice fit well with the measured data, and the better fit is corresponding to the deeper sea ice.


Introduction
Arctic sea ice cover affects the exchange of heat, energy, mass, and momentum between the atmosphere and the Arctic ocean, and it is a major component of the Arctic environment.Parkinson and Cavalieri 1 suggested that sea ice in the Arctic is a key indicator of climate change.Sea ice has been the principal threat to the development of the Arctic cruise which is cheap, reliable, and capable of year-around operation.Extremely formidable sea ice cover can restrict the operation of surface ships to a few months each summer, and even then there are areas where icebreaker assistance is still required.Meanwhile, sea ice is also a significant threat for the safe design of offshore structures and harbor facilities.Therefore, the research on sea ice has aroused interests of the scientists in many fields.Unfortunately, the physical parameters of sea ice have been obtained mainly by field measurement until now, which are spare and unsatisfactory due to the difficulties in situ, especially during the polar winter.The parameter identification method can reinforce the field data for improved understanding of the physical mechanism and variabilities of parameters of sea ice because sea ice temperature can be measured continuously and automatically.Moreover, it can help to forecast the variabilities of sea ice and climate 2-4 .
Many researchers have been devoted to the establishment and improvement of the thermodynamic model of sea ice system.Maykut and Untersteinter 5 first formally proposed a comprehensive one-dimensional sea ice thermodynamic model named MU model, but the corresponding fluxes of atmosphere and ocean were not involved.Semtner 6 simplified MU model, his model has been widely used in climate simulations, while lacks universal and effectiveness.Parkinson and Washington 7 made an improvement of Semtner's model 6 , extended the thermodynamic model to a three-dimensional model on a large scale, and considered especially the impact of the leads, which are stretches of open water within fields of sea ice.Hibler 8 established a thermal-dynamic model and carried out first the numerical simulation using finite difference method.Numerical simulation and parameter identification of sea ice thermodynamic system have been paid increasing attention over the past several years, and the consideration of the thermodynamic process of sea ice has been more detailed and comprehensive 9-12 .The thickness of sea ice is one of the most important parameters in thermodynamic and dynamic models of sea ice, and it describes the vertical scale of sea ice.Unfortunately, owing to the randomness and variability of the distribution of sea ice, the continuous and accurate measurement of the thickness is very difficult.Meanwhile, the continuous changes of sea ice thickness in the numerical simulation of large-scale and long process may lead to comparatively large deviation between the simulated and measured sea ice temperatures.The introduction of the change of sea ice thickness in the thermodynamic system thus cannot only describe the thermodynamic behavior of sea ice better but also make the numerical simulation more accurately, while there is few consideration of it at present 13 .
The thermodynamic process of sea ice is a heat conductivity process in nature which can be described by a parabolic partial differential equation generally.Therefore, we describe the variability of sea ice temperature in space and time by a parabolic distributed parameter system, with sea ice temperature as the state variable.A one-dimensional threelayer thermodynamic system coupled by snow, ice, and sea water is considered, and the thicknesses of snow and sea ice are identified.A distributed parameter system in a timedependent domain is established and transferred to an abstract parabolic evolution system by some mathematical methods, and the existence and uniqueness of the weak solution of the evolution equation are proved.By taking the thicknesses of snow and sea ice as identified variables and the deviation between the calculated and measured temperature of sea ice as the performance criterion, an optimal model with state constraints is presented.Then, the existence of the optimal parameter is discussed, and the one-order necessary condition of the optimality is given.Finally, based on the nonoverlapping domain decomposition method and semi-implicit difference schemes, an optimization algorithm is constructed and applied to study the thickness and temperature of the Arctic sea ice.Through numerical simulation, we obtained the change characteristics of snow and ice thicknesses and temperature distribution in sea ice.This paper is organized as follows.In Section 2, we present a thermodynamic system coupled by snow, ice, and sea water.Section 3 derives some important properties of a bilinear functional and proves the existence and uniqueness of the weak solution of the evolution equation.Section 4 establishes an optimal model with state constraints and proves the strong continuity of the weak solution with respect to the parameters and the existence of the optimal parameter.Moreover, the one-order necessary condition of the optimality is derived.An numerical optimization algorithm is constructed and applied to numerical simulation in Section 5, and numerical results are presented in Section 6.Some conclusions are presented in the final section.

The Coupled Thermodynamic System of Sea Ice
Consider the thermodynamic system of sea ice coupled by snow, ice, and sea water layers Figure 1 , denoted by snow-ice-water system.Since the gradient variation of sea ice temperature in the vertical direction is far greater than that in the horizontal direction, we only consider the heat flux in the vertical direction.Let t denote the time unit: s , T 0 < T < ∞ the final time, and 0, T the observation period.Set I T : 0, T .
Let x axis be the depth direction of the snow-ice-water system unit: m , and the point on snow surface at the initial measured time the coordinate origin Figure 1 .Let h 1 t , h 2 t , and h 3 t be the thicknesses of snow, ice, and sea water layers at the time t ∈ I T , respectively, then the total thickness at the time t can be expressed as Obviously, the thickness of sea water at the time t can be expressed as According to the physical properties of sea ice, we have h 1 t , h 2 t ∈ C 1 I T ; R .Let h 1 0 h 10 and h 2 0 h 20 be the initial thicknesses of snow and sea ice, respectively, which are known constants, then the initial thickness of sea water layer is Let the initial space domains of snow, ice, and sea water layers be Ω 10 0, h 10 , Ω 20 h 10 , h 10 h 20 and Ω 30 h 10 h 20 , L 0 , respectively, then the initial space domain of the system can be described as obviously, Ω 0 is a nonempty, bounded, and connected set.Additionally, let , and Ω 3 t h 1 t h 2 t , L t be the space domains of snow, ice, and sea water layers at the time t ∈ I T , Snow layer

Ice layer
Sea water layer respectively, then the space domain of the snow-ice-water system at the time t ∈ I T can be expressed as Apparently, the measured space-time domain denoted by Q T Ω t × I T is opened and bounded.
By the energy conversation and the Fourier law, the coupled thermodynamic sea ice model is thus described by the following parabolic partial differential dynamic system where ρ ρ i is the density, c c i the specific heat and K K i the thermal conductivity, T x, t the temperature of the snow-ice-water system at the depth x and the time t unit: Kelvin , T 0 x the temperature at the depth x and the time t 0, and g x, t g i x, t the heat source term.x, t ∈ Ω i × I T , i 1, 2, 3, denotes the snow, ice, and sea water layer, respectively.Q 1 t and T 3 t are both known functions on the time t.
The system 2.6 is piecewise smooth owing to the difference of the densities, specific heats, and thermal conductivities in the three layers.By the physical properties of sea ice, we give the following assumptions.
A1 The thermodynamic parameters ρ, c, and K remain unchanged during the measured period.
A2 The upper and lower boundaries remain unmoved, namely, the total thickness of the system is a constant L t L .
A3 T x, t is continuous on Q T .
A4 h 1 t , h 2 t ∈ C 2 I T ; R are bounded functions with positive values, and there exist h jl > 0 and h ju > 0 satisfying h jl ≤ h ju , such that h jl ≤ h j ≤ h ju , j 1, 2.
Additionally, functions in the system 2.6 should be differential on x or t in some sense; thus, we assume that According to the assumptions A1 and A2 , we set α ρc/K, and y x, t T x, t − x − L Q 1 t − T 3 t , then the system 2.6 can be described as the following homogeneous boundary value problem denoted by HBP : where

Properties of the Coupled Thermodynamic System
For convenience, we let , ∀ϕ ∈ H u t .

3.1
The inner product and the induced norm on V u t are defined by where ϕ x and ψ x are the derivatives of ϕ and ψ on the depth x, respectively.The dual paring between V u t and V u t is expressed as Let ξ ∈ V u t , we multiply both sides of 2.7 by ξ and integrate them over Ω u t : Using the method of integration by parts and the homogeneous boundary conditions, 3.4 can be described as For any t ∈ I T , we define a functional on V u t × V u t as a t, u; ϕ, ψ The properties of the functional a t, u; ϕ, ψ are presented in the following.
Property 1. Suppose that assumptions A1 -A5 hold, then 1 a t, u; ϕ, ψ defined by 3.6 is a bilinear functional on V u t × V u t .
2 For all ϕ, ψ ∈ V u t , a •, u; ϕ, ψ is a measurable function on I T , and there exists c > 0, such that |a t, u; ϕ, ψ

3.7
So a t, u; ϕ, ψ is a bilinear functional on V u t × V u t .
2 Obviously, for all ϕ, ψ ∈ V u t , a t, u; ϕ, ψ is a measurable function on I T according to 3.6 , and a t, u; ϕ, ψ
Proof.The conclusion follows directly from 3.6 and Garding inequality 14 .
From the Properties 1 and 2 and Lax-Milgram theorem 15 , there exists a bounded continuous linear operator A t, u ∈ L V u t , V u t , such that

3.10
According to 3.9 , the Properties 1 and 2, we can get the following property.
Property 3. Suppose that assumptions A1 -A5 hold, then for all t ∈ I T , for all ϕ, ψ ∈ V u t , and for all u ∈ U ad I T , there exist M 1 > 0 and M 2 > 0, such that

3.11
For any t ∈ I T , ϕ ∈ V u t , and u ∈ U ad I T , set H t, u; ϕ Ω u t ϕp x, t dx, then H t, u; ϕ is a continuous linear functional on V u t .According to Riesz representation theorem 14 , there exists a unique element f t; u ∈ V u t , such that H t, u; ϕ f t; u , ϕ V u t ,V u t , ∀ϕ ∈ V u t .

3.12
From the above, for each u ∈ U ad I T , the system HBP can be written as a parabolic evolution equation denoted by PEE : where t ∈ I T , u ∈ U ad I T , and V u t dt < ∞} with the inner product defined as w, w T 0 w t , w t V t dt.
3.14 Obviously, the above process is reversible, the system PEE is thus equivalent to HBP .The definition of weak solution of the system PEE is given in the following.From the above theorem, we can also know that the system HBP has a unique solution y t, x; u ∈ L 2 I T , Ω u t ; H u t ∩ C I T , Ω u t ; H u t .

Identification Model
We define the performance criterion as where T t, x; u is the temperature function obtained from the system 2.6 , T t, x is the observed temperature.Let f 0 u y t, x; u − y t, x 2 , then 4.1 is expressed as The goal of this study is to make the temperature T t, x; u obtained from the system 2.6 fits the measured data as far as possible.Then, the identification model of the system HBP is expressed as min J u

SIP
where S U ad Q T {y t, x; u | y t, x; u is the solution of system HBP corresponding to u ∈ U ad I T }, and J u is defined by 4.1 .
Obviously, from 4.1 -4.2 , we can conclude that the problem SIP is equivalent to the following problem

SIPO
where S U ad Q T {T t, x; u | T t, x; u is the solution of system 2.6 corresponding to u ∈ U ad I T }.

Strong Continuity of the Weak Solution on Parameters
Theorem 4.1.Suppose that assumptions (A1)-(A5) hold, then the mapping u → y t, x; u is strongly continuous.
Proof.For any given parameter u 0 ∈ U ad I T , let {u n } ⊂ U ad I T be a feasible parameter sequence, such that u n − u 0 → 0 as n → ∞, and y n and y 0 the solutions of PEE corresponding to u n and u 0 , respectively.Set w n y n − y 0 , p n p t, x; u n , p 0 p t, x; u 0 , then we obtain the following system

4.4
Multiply both sides of the first equation of 4.4 by w n and integrate them over Ω u t × 0, t , then we have Using Gronwall inequality, we obtain From all the above and u n −u 0 → 0, we get |p n −p 0 | → 0, the mapping u → y t, x; u is thus strongly continuous.

Existence of Optimal Parameter
Theorem 4.2.Suppose that assumptions (A1)-(A5) hold, then there exists at least one optimal parameter u * ∈ U ad I T satisfies the identification problem (SIP).
Proof.Obviously, f 0 u is continuous on V u t , and, from 4.3 , J u I T Ω u t f 0 u dx dt ≥ 0. By Theorem 4.1, the mapping u → y t, x; u is continuous.Hence, the mapping u → J u is continuous on U ad I T .Since U ad I T is a nonempty, bounded, and closed set, there exists u * ∈ U ad I T , such that for all u ∈ U ad I T , J u * ≤ J u , e.g., u * ∈ U ad I T is an optimal parameter, thus the desired result is obtained.

Necessary Optimality Condition
Let u * ∈ U ad I T be the optimal parameter of the problem SIP , following the convex compactness of U ad I T and the continuity of J u on U ad I T , we can prove that J u • is Gateaux differentiable at u * ∈ U ad I T and its Gateaux derivative DJ u * • exists.Hence, we can get the following necessary conditions for optimality.Theorem 4.3.Suppose that assumptions (A1)-(A5) are valid, and let u * ∈ U ad I T be the optimal parameter of the system (SIP), then u * • satisfies the following inequality:

Optimization of the Parameter Identification Model
In this section, we aim at constructing an optimal algorithm to solve the parameter identification problem SIPO .On the domain Ω i t , we decompose 2.6 as with the initial condition and the boundary and penetrating conditions
According to the partial differential theory, we can obtain the following theorem.
Theorem 5.1.Suppose that assumptions (A1)-(A5) are valid, then, for any u ∈ U ad I T , the subsystem IBP i has a unique weak solution T i x, t; u ∈ C Ω i t , I T , U ad ; R which depends continuously on u ∈ U ad I T .
Let M, N be the numbers of the measured spatial and temporal spots, respectively, x k the measured depth, t j the measured time, T x k , t j ; u the calculated temperature of sea ice at depth x k and time t j from the system 2.6 , T x k , t j the measured temperature data at depth x k and time t j , then the practical identification problem can be expressed as

SIPOD i
Obviously, the performance function of SIPOD can be described as From Theorem 4.2, we can get that there also exists at least one optimal parameter satisfying the subproblem SIPOD i . Set T x k , t j ; u , i 1, 2, 3, T x k , t j , i 1, 2, 3, then, the identification problem is expressed as and the subidentification problem is min J bi u SIPOD bi Theorem 5.2.If u * is the optimal parameter of SIPOD, then u * is also that of SIPOD b .

Numerical Optimization Algorithm
In this section, the semi-implicit finite difference scheme 15 is employed to discretize, and the nonoverlapping Schwarz alternating direction method to solve 5.1 .And related to the measured temperature data {T z k , t j }, a numerical optimization algorithm is constructed with the following steps.
Step 1. Select starting points u 0 u 0 1 , u 0 2 ∈ U ad I T , n directions e 1 , e 2 , . . ., e n , starting step length Δu Δu 1 , Δu 2 , Δu 1 > 0, Δu 2 > 0, acceleration factor α, accuracy ε > 0, and the maximum iteration number k max , set v 0 u 0 , k j 0; Step 2. Calculate the numerical solution T x, t; v j , by nonoverlapping Schwarz alternating direction method and semi-implicit difference scheme. Step Step 4. If j < n, set j j 1, go to Step 3; for the case j n, if J b v n 1 ≤ J b u k , go to Step 5, else, go to Step 6. Step and j 0, go to Step 2.

Numerical Results
Using the above optimal algorithm, we can obtain the temperature distribution T x k , t j ; u of sea ice in the Arctic.The data set was measured by a monitor buoy installed on floe ice in the Arctic.
The comparison between the calculated and measured sea ice temperatures from November 1, 2003 to February 29, 2004 was shown in Figure 2, with the test temperatures data being smoothed for the unomitted transparent errors in some of the original data.The measured and simulated data in Figure 2 both indicate a general increasing trend of the temperature of sea ice, while a decreasing increment with increasing depth of sea ice.Meanwhile, because the data had been measured during the freezing period of sea ice, the temperature generally decreases with increasing time.The calculated temperature fits well with the measured data; furthermore, the better fit between the calculated and measured temperatures is corresponding to the deeper sea ice.The temperature deviation between the simulated and measured temperatures is 0.361, which is a much small error for the temperature of sea ice.It can be concluded that the numerical results by our method reflect the actual variation of the sea ice temperature distribution and approach the measured data well.

Conclusions
In this study, we aim to simulate the temperature of sea ice with thicknesses of snow and sea ice as identified parameters.We established a parabolic distributed system of the temperature field of sea ice with a time-dependent domain in the Arctic sea ice and proposed its parameter identification problem with thicknesses of snow and sea ice as identified variables and the deviation between the calculated and measured sea ice temperatures as the performance criterion.And we proved the unique existence of the weak solution of the thermodynamic system, discussed the existence of the optimal parameter, and derived the one-order necessary condition.Finally, based on the nonoverlapping domain decomposition method and semi-implicit difference scheme, we constructed an optimization algorithm to simulate the sea ice temperature in the Arctic.The results Figure 2 show that the temperature of sea ice increases with increasing depth of sea ice, while, because the measurement had been made during the freezing period of sea ice, it decreases with increasing time.
The simulated temperatures of sea ice fitted well with the measured data, and the deviation between the simulated and measured temperature of sea ice was slight.It indicates that the numerical results by our method reflect the actual variation of the sea ice temperature distribution in space and time.Meanwhile, it is illuminated that the algorithm proposed in the present study is valid and more suitable for the deeper sea ice.
Nevertheless, it is important to note that our mathematical framework was verified only using in situ data from the Arctic, and the variabilities of the physical parameters of the coupled sea ice system, such as density, specific heat, were not considered, more measurements would be required for a true verification of this method in the future work.However, we believe that our method is a promising approach worthy of further studies under different environmental conditions in the Arctic as well as in the Antarctic.

Figure 1 :
Figure 1: Sketch of the sea ice thermodynamic model.

Figure 2 :
Figure 2: Comparison of calculated and measured temperatures.
is the admissible parameter set, Γ is a known nonempty, bounded, closed convex subset on R 2 .And let Ω u t Ω 1u t ∪ Ω 2u t ∪ Ω 3u t be the spatial domain of the distributed parameter system HBP at the time t, where Ω iu t Ω i t , i 1, 2, 3.For any t ∈ I T and u ∈ U ad I T , let H u t L 2 Ω u t ; R be a separable Hilbert space, Journal of Applied Mathematics is continuous, and V u t is dense in H u t .The inner product and the induced norm on H u t V u t , respectively, then V u t ⊂ H u t .Let •, • H u t and •, • V u t denote the inner product on H u t and V u t , respectively, H u t and V u t the dual spaces of H u t and V u t , respectively, and •, • V u t ,V u t the dual paring between V u t and V u t , then V u t , H u t , V u t is a Gelftriple space with V u t → H u t ≡ H u t → V u t , namely, the embedding V u t → H u t 6 L 2 I T , Ω u t ; H u t ∩ C I T , Ω u t ; H u t which depends continuously on f and y 0 .By the definition of V u t , we know that ψ ∈ V u t satisfies the Dirichlet condition ψ x 0, x ∈ ∂Ω u t /∂Ω • u t for any given t ∈ I T .The system PEE thus has a unique weak solution which depends continuously on f and y 0 15 .
a •, u; y •, •; u , ξ f •; u , ξ H u t , 3.15 where ξ ∈ V u t and u ∈ U ad I T .Theorem 3.2.Suppose that assumptions (A1)-(A5), Properties 1 and 3 hold.For all t ∈ I T , let Ω u t x u t, Ω 0 is bounded, and there exists a countable base {ϕ 1 t , ϕ 2 t , . ..} on V u t , such that ϕ it t ∈ V u t , i 1, 2, . . . .Then, the system (PEE) has a unique weak solution y t, x; u ∈ {T i x, t; u | T i x, t; u is the solution of the subsystem IBP i corresponding to u ∈ U ad I T }, i 1, 2, 3, and for the numerical implementation, the problem SIPOD can be decomposed to the following three connected subproblems