# Simulation and Domain Identification of Sea Ice Thermodynamic System

Authors: Bing Tan; Peng Lu; Zhijun Li; Enmin Feng; Juan Wang
Journal: Journal of Applied Mathematics
 (2012)

---

## Abstract

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.


---

## Body

## 1. 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 three-layer 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 time-dependent 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 Section2, 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.


## 2. The Coupled Thermodynamic System of Sea Ice

Consider the thermodynamic system of sea ice coupled by snow, ice, and sea water layers (Figure1), 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 IT:=(0,T).Figure 1
Sketch of the sea ice thermodynamic model.Letx 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 h1(t), h2(t), and h3(t) be the thicknesses of snow, ice, and sea water layers at the time t∈IT, respectively, then the total thickness at the time t can be expressed as
(2.1)L(t)=h1(t)+h2(t)+h3(t).Obviously, the thickness of sea water at the timet can be expressed as
(2.2)h3(t)=L(t)-h1(t)-h2(t).
According to the physical properties of sea ice, we have h1(t), h2(t)∈C1(IT;R).Leth1(0)=h10 and h2(0)=h20 be the initial thicknesses of snow and sea ice, respectively, which are known constants, then the initial thickness of sea water layer is
(2.3)h3(0)=h30=L(0)-h10-h20.Let the initial space domains of snow, ice, and sea water layers beΩ10=(0,h10], Ω20=(h10,h10+h20] and Ω30=(h10+h20,L(0)), respectively, then the initial space domain of the system can be described as
(2.4)Ω0=Ω10∪Ω20∪Ω30
obviously, Ω0 is a nonempty, bounded, and connected set.Additionally, letΩ1(t)=(0,h1(t)], Ω2(t)=(h1(t),h1(t)+h2(t)], and Ω3(t)=(h1(t)+h2(t),L(t)) be the space domains of snow, ice, and sea water layers at the time t∈IT, respectively, then the space domain of the snow-ice-water system at the time t∈IT can be expressed as
(2.5)Ω(t):=Ω1(t)∪Ω2(t)∪Ω3(t)=(0,L(t))⊂R.Apparently, the measured space-time domain denoted byQT=Ω(t)×IT 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(2.6)ρc⋅∂T(x,t)∂t=∂∂x(K⋅∂T(x,t)∂x)+g(x,t),(x,t)∈QT,T(x,0)=T0(x),x∈Ω¯0,∂T(x,t)∂t|x=0=Q1(t),t∈I¯T,T(x,t)|x=L(t)=T3(t),t∈I¯T,
where ρ=ρi is the density, c=ci the specific heat and K=Ki the thermal conductivity, T(x,t) the temperature of the snow-ice-water system at the depth x and the time t (unit: Kelvin), T0(x) the temperature at the depth x and the time t=0, and g(x,t)=gi(x,t) the heat source term. (x,t)∈Ωi×IT,  i=1,2,3, denotes the snow, ice, and sea water layer, respectively. Q1(t) and T3(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)
h1(t), h2(t)∈C2(IT;R) are bounded functions with positive values, and there exist hjl>0 and hju>0 satisfying hjl≤hju, such that hjl≤hj≤hju,j=1,2.Additionally, functions in the system (2.6) should be differential on x or t in some sense; thus, we assume that (A5)
T0(x)∈C2(Ωi∘;R), T3(t), Q1(t)∈C2(IT;R), gi(x,t)∈C2(Ωi∘×IT;R), i=1,2,3.According to the assumptions (A1) and (A2), we setα=ρc/K, and y(x,t)=T(x,t)-(x-L)Q1(t)-T3(t), then the system (2.6) can be described as the following homogeneous boundary value problem denoted by (HBP):
(2.7)∂y(x,t)∂t=α∂2y(x,t)∂x2+p(x,t),(x,t)∈QT,(2.8)y(x,0)=y0(x),x∈Ω¯0,(2.9)∂y(x,t)∂t|x=0=0,t∈I¯T,(2.10)y(x,t)|x=L=0,t∈I¯T,
where
(2.11)p(x,t)=1ρcg(x,t)-(x-L)Q1′(t)-T3′(t),y0(x)=T0(x)-(x-L)Q1(0)-T3(0).


## 3. Properties of the Coupled Thermodynamic System

For convenience, we letu1=h1(t), u2=h2(t), u=(u1,u2), and Uad(I¯T):={u(·)∈C2(IT;R+)∣u(t):I¯T→Γ is measurable, hjl≤uj≤hju,j=1,2}⊂C2(IT;R+) is the admissible parameter set, Γ is a known nonempty, bounded, closed convex subset on R2. 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 anyt∈I¯T and u∈Uad(I¯T), let Hu(t)=L2(Ωu(t);R) be a separable Hilbert space, Vu(t)={φ∣φ∈Hu′(t), ∂φ/∂x|x=0=φ|x=L=0} with norms denoted by |·|Hu(t) and ∥·∥Vu(t), respectively, then Vu(t)⊂Hu(t). Let (·,·)Hu(t) and 〈·,·〉Vu(t) denote the inner product on Hu(t) and Vu(t), respectively, Hu′(t) and Vu′(t) the dual spaces of Hu(t) and Vu(t), respectively, and 〈·,·〉Vu(t),Vu′(t) the dual paring between Vu(t) and Vu′(t), then (Vu(t),Hu(t),Vu′(t)) is a Gelf-triple space with Vu(t)↪Hu(t)≡Hu′(t)↪Vu′(t), namely, the embedding Vu(t)↪Hu(t) is continuous, and Vu(t) is dense in Hu(t). The inner product and the induced norm on Hu(t) are defined by
(3.1)(φ,ψ)Hu(t)=∫Ωu(t)‍φψdx,∀φ,ψ∈Hu(t),|φ|Hu(t)=(∫Ωu(t)‍φ2dx)1/2,∀φ∈Hu(t).
The inner product and the induced norm on Vu(t) are defined by
(3.2)〈φ,ψ〉Vu(t)=∫Ωu(t)‍(φψ+φxψx)dx,∀φ,ψ∈Vu(t),‖φ‖Vu(t)=(∫Ωu(t)‍(φ2+φx2)dx)1/2,∀φ∈Vu(t),
where φx and ψx are the derivatives of φ and ψ on the depth x, respectively.The dual paring betweenVu′(t) and Vu(t) is expressed as
(3.3)〈φ,ϕ〉Vu′(t),Vu(t)=∫Ωu(t)‍φϕdx.Letξ∈Vu(t), we multiply both sides of (2.7) by ξ and integrate them over Ωu(t):(3.4)∫Ωu(t)‍ξ∂y∂tdx=∫Ωu(t)‍ξα∂2y∂x2dx+∫Ωu(t)‍ξp(x,t)dx.Using the method of integration by parts and the homogeneous boundary conditions, (3.4) can be described as
(3.5)∫Ωu(t)‍ξ∂y∂tdx+∫Ωu(t)‍α∂ξ∂x∂y∂xdx=∫Ωu(t)‍ξp(x,t)dx.For anyt∈IT, we define a functional on Vu(t)×Vu(t) as
(3.6)a(t,u;φ,ψ)=∫Ωu(t)‍α∂φ∂x⋅∂ψ∂xdx.The properties of the functionala(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 Vu(t)×Vu(t).(2)
For allφ, ψ∈Vu(t), a(·,u;φ,ψ) is a measurable function on IT, and there exists c>0, such that |a(t,u;φ,ψ)|≤c·∥φ∥Vu(t)·∥ψ∥Vu(t).Proof.
(1) For allλ1, λ2∈R, for all φ1, φ2, ψ1, ψ2∈Vu(t), we have
(3.7)a(t,u;λ1φ1+λ2φ2,ψ)=∫Ωu(t)‍α∂λ1φ1+λ2φ2∂x⋅∂ψ∂xdx=λ1a(t,u;φ1,ψ)+λ2a(t,u;φ2,ψ),a(t,u;φ,λ1ψ1+λ2ψ2)=∫Ωu(t)‍α∂φ∂x⋅∂λ1ψ1+λ2ψ2∂xdx=λ1a(t,u;φ,ψ1)+λ2a(t,u;φ,ψ2).
So a(t,u;φ,ψ) is a bilinear functional on Vu(t)×Vu(t).
(2) Obviously, for allφ, ψ∈Vu(t), a(t,u;φ,ψ) is a measurable function on IT according to (3.6), and
(3.8)|a(t,u;φ,ψ)|=|∫Ωu(t)‍α∂φ∂x⋅∂ψ∂xdx|≤α(∫Ωu(t)‍φx2dx)1/2⋅(∫Ωu(t)‍ψx2dx)1/2≤α⋅‖φ‖Vu(t)⋅‖ψ‖Vu(t).
Taking c=α, then |a(t,u;φ,ψ)|≤c·∥φ∥Vu(t)·∥ψ∥Vu(t).Property 2.
Suppose that assumptions (A1)–(A5) hold, then, for the bilinear functionala(t,u;φ,ψ) defined by (3.6), there exist α0>0 and λ>0, such that |a(t,u;φ,φ)|+λ|φ|Hu(t)2≥α0∥φ∥Vu(t)2, for all φ∈Vu(t),∀t∈I¯T.Proof.
The conclusion follows directly from (3.6) and Garding inequality [14].From the Properties1 and 2 and Lax-Milgram theorem [15], there exists a bounded continuous linear operator A(t,u)∈L(Vu(t),Vu′(t)), such that
(3.9)a(t,u;φ,ψ)=〈A(t,u)φ,ψ〉Vu′(t),Vu(t),∀t∈I¯T,∀u∈Uad(I¯T),(3.10)A(t,u)φ=-α∂2φ∂x2,∀φ∈Vu(t),∀u∈Uad(I¯T).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 allt∈I¯T, for all φ,ψ∈Vu(t), and for all u∈Uad(I¯T), there exist M1>0 and M2>0, such that
(3.11)〈A(t,u)φ,φ〉Vu′(t),Vu(t)≥M1‖φ‖Vu(t)2,〈A(t,u)φ,ψ〉Vu′(t),Vu(t)≥M2‖φ‖Vu(t)2‖ψ‖Vu(t)2.
For anyt∈I¯T, φ∈Vu(t), and u∈Uad(I¯T), set H(t,u;φ)=∫Ωu(t)‍φp(x,t)dx, then H(t,u;φ) is a continuous linear functional on Vu(t). According to Riesz representation theorem [14], there exists a unique element f(t;u)∈Vu′(t), such that
(3.12)H(t,u;φ)=〈f(t;u),φ〉Vu′(t),Vu(t),∀φ∈Vu(t).
From the above, for eachu∈Uad(I¯T), the system (HBP) can be written as a parabolic evolution equation denoted by (PEE):
(3.13)yt+A(t,u)y=f(t;u)y(0)=y0,y0∈Hu(0),
where t∈I¯T, u∈Uad(I¯T), and f∈L2(IT;Vu′(·)), L2(IT;Vu(t))={w∣∫0T‍∥w(t)∥Vu(t)2dt<∞} with the inner product defined as
(3.14)(w,w¯)=∫0T‍〈w(t),w¯(t)〉Vu(t)dt.
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.Definition 3.1.
A functiony(t,x;u)∈L2(IT,Ωu(t);Hu(t))∩C(IT,Ωu(t);Hu′(t)) is a weak solution of (PEE), if y(t,x;u) satisfies
(3.15)〈dy(⋅,⋅;u)dt,ξ〉Vu′(t),Vu(t)+a(⋅,u;y(⋅,⋅;u),ξ)=(f(⋅;u),ξ)Hu(t),
where ξ∈Vu(t) and u∈Uad(I¯T).Theorem 3.2.
Suppose that assumptions (A1)–(A5), Properties1 and 3 hold. For all t∈I¯T, let Ωu(t)=xu(t,Ω0) is bounded, and there exists a countable base {φ1(t),φ2(t),…} on Vu(t), such that φit(t)∈Vu′(t), i=1,2,…. Then, the system (PEE) has a unique weak solution y(t,x;u)∈L2(IT,Ωu(t);Hu(t))∩C(IT,Ωu(t);Hu′(t)) which depends continuously on f and y0.Proof.
Let∂Ωu∘(t)={x∈Ωu(t):h(t,x,u)·n(x)=0}, here n(x) is the outward normal vector of ∂Ωu(t) at x. By the definition of Vu(t), we know that ψ∈Vu(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 y0 [15].
From the above theorem, we can also know that the system (HBP) has a unique solutiony(t,x;u)∈L2(IT,Ωu(t);Hu(t))∩C(IT,Ωu(t);Hu′(t)).


## 4. Parameter Identification of the Coupled System

### 4.1. Identification Model

We define the performance criterion as(4.1)J(u)=∫0T‍∫Ωu(t)‍(y(t,x;u)-y¯(t,x))2dxdt
with
(4.2)y(t,x;u)=T(t,x;u)-(x-L)Q1(t)-T3(t),y¯(t,x)=T¯(t,x)-(x-L)Q1(t)-T3(t),
where T(t,x;u) is the temperature function obtained from the system (2.6), T¯(t,x) is the observed temperature. Let f0(u)=(y(t,x;u)-y¯(t,x))2, then (4.1) is expressed as
(4.3)J(u)=∫0T‍∫Ωu(t)‍f0(u)dxdt.The goal of this study is to make the temperatureT(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
(SIP)min⁡J(u)s.t.y(t,x;u)∈SUad(QT)u∈Uad(I¯T),
where SUad(QT)={y(t,x;u)∣y(t,x;u) is the solution of system (HBP) corresponding to u∈Uad(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)min⁡J(u)=∫0T‍∫Ωu(t)‍(T(t,x;u)-T¯(t,x))2dxdts.t.T(t,x;u)∈SUad′(QT)u∈Uad(I¯T),
where SUad′(QT)={T(t,x;u)∣T(t,x;u) is the solution of system (2.6) corresponding to u∈Uad(I¯T)}.


### 4.2. Strong Continuity of the Weak Solution on Parameters

Theorem 4.1.
Suppose that assumptions (A1)–(A5) hold, then the mappingu→y(t,x;u) is strongly continuous.Proof.
For any given parameteru0∈Uad(I¯T), let {un}⊂Uad(I¯T) be a feasible parameter sequence, such that ∥un-u0∥→0 as n→∞, and yn and y0 the solutions of (PEE) corresponding to un and u0, respectively. Set wn=yn-y0, pn=p(t,x;un), p0=p(t,x;u0), then we obtain the following system
(4.4)∂wn∂t-α∂2wn∂x2=pn-p0w0(x,0)=0∂wn∂t|x=0=0wn|x=L=0.
Multiply both sides of the first equation of (4.4) by wn and integrate them over Ωu(t)×[0,t], then we have
(4.5)∫0t‍(wn,∂wn∂s)ds+∫0t‍α(s,un;wn,wn)ds=∫0t‍(wn,pn-p0)ds.
For the first term of (4.5) on the left side, we have
(4.6)∫0t‍(wn,∂wn∂s)ds=12(|wn(t)|2-|wn(0)|2)=12|wn(t)|2.
For the second term of (4.5) on the left side, according to Property 2, we have
(4.7)∫0t‍α(s,un;wn,wn)ds+λ∫0t‍|wn|2ds≥α0∫0t‍|wn|2ds,
where α0 and λ are the same with that of Property 2.
Using the elementary inequalityab≤(1/2ɛ)a2+(ɛ/2)b2 and taking ɛ=α0, then, for the right terms of (4.5), we have
(4.8)∫0t‍(wn,pn-p0)ds≤α02∫0t‍|wn|2ds+12α0∫0t‍|pn-p0|2ds.
LetYn(t)=|wn|2+α0∫0t‍|wn|2ds. Substituting (4.6)–(4.8) into (4.5), then we have
(4.9)Yn(t)≤1α0∫0t‍|pn-p0|2ds+(α0+2λ)∫0t‍|wn|2ds.
Using Gronwall inequality, we obtain(4.10)Yn(t)≤1α0eα0T∫0t‍|pn-p0|2ds.
From all the above and∥un-u0∥→0, we get |pn-p0|→0, the mapping u→y(t,x;u) is thus strongly continuous.


### 4.3. Existence of Optimal Parameter

Theorem 4.2.
Suppose that assumptions (A1)–(A5) hold, then there exists at least one optimal parameteru*∈Uad(I¯T) satisfies the identification problem (SIP).Proof.
Obviously,f0(u) is continuous on Vu(t), and, from (4.3), J(u)=∫IT‍∫Ωu(t)‍f0(u)dxdt≥0. By Theorem 4.1, the mapping u→y(t,x;u) is continuous. Hence, the mapping u→J(u) is continuous on Uad(I¯T). Since Uad(I¯T) is a nonempty, bounded, and closed set, there exists u*∈Uad(I¯T), such that for all u∈Uad(I¯T), J(u*)≤J(u), e.g., u*∈Uad(I¯T) is an optimal parameter, thus the desired result is obtained.


### 4.4. Necessary Optimality Condition

Letu*∈Uad(I¯T) be the optimal parameter of the problem (SIP), following the convex compactness of Uad(I¯T) and the continuity of J(u) on Uad(I¯T), we can prove that J(u(·)) is Gateaux differentiable at u*∈Uad(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 letu*∈Uad(I¯T) be the optimal parameter of the system (SIP), then u*(·) satisfies the following inequality:
(4.11)DJ(u*;u-u*)=lim⁡γ→0+J(u*+γ(u-u*))-J(u*)γ≥0,∀u∈Uad(I¯T).


## 4.1. Identification Model

We define the performance criterion as(4.1)J(u)=∫0T‍∫Ωu(t)‍(y(t,x;u)-y¯(t,x))2dxdt
with
(4.2)y(t,x;u)=T(t,x;u)-(x-L)Q1(t)-T3(t),y¯(t,x)=T¯(t,x)-(x-L)Q1(t)-T3(t),
where T(t,x;u) is the temperature function obtained from the system (2.6), T¯(t,x) is the observed temperature. Let f0(u)=(y(t,x;u)-y¯(t,x))2, then (4.1) is expressed as
(4.3)J(u)=∫0T‍∫Ωu(t)‍f0(u)dxdt.The goal of this study is to make the temperatureT(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
(SIP)min⁡J(u)s.t.y(t,x;u)∈SUad(QT)u∈Uad(I¯T),
where SUad(QT)={y(t,x;u)∣y(t,x;u) is the solution of system (HBP) corresponding to u∈Uad(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)min⁡J(u)=∫0T‍∫Ωu(t)‍(T(t,x;u)-T¯(t,x))2dxdts.t.T(t,x;u)∈SUad′(QT)u∈Uad(I¯T),
where SUad′(QT)={T(t,x;u)∣T(t,x;u) is the solution of system (2.6) corresponding to u∈Uad(I¯T)}.


## 4.2. Strong Continuity of the Weak Solution on Parameters

Theorem 4.1.
Suppose that assumptions (A1)–(A5) hold, then the mappingu→y(t,x;u) is strongly continuous.Proof.
For any given parameteru0∈Uad(I¯T), let {un}⊂Uad(I¯T) be a feasible parameter sequence, such that ∥un-u0∥→0 as n→∞, and yn and y0 the solutions of (PEE) corresponding to un and u0, respectively. Set wn=yn-y0, pn=p(t,x;un), p0=p(t,x;u0), then we obtain the following system
(4.4)∂wn∂t-α∂2wn∂x2=pn-p0w0(x,0)=0∂wn∂t|x=0=0wn|x=L=0.
Multiply both sides of the first equation of (4.4) by wn and integrate them over Ωu(t)×[0,t], then we have
(4.5)∫0t‍(wn,∂wn∂s)ds+∫0t‍α(s,un;wn,wn)ds=∫0t‍(wn,pn-p0)ds.
For the first term of (4.5) on the left side, we have
(4.6)∫0t‍(wn,∂wn∂s)ds=12(|wn(t)|2-|wn(0)|2)=12|wn(t)|2.
For the second term of (4.5) on the left side, according to Property 2, we have
(4.7)∫0t‍α(s,un;wn,wn)ds+λ∫0t‍|wn|2ds≥α0∫0t‍|wn|2ds,
where α0 and λ are the same with that of Property 2.
Using the elementary inequalityab≤(1/2ɛ)a2+(ɛ/2)b2 and taking ɛ=α0, then, for the right terms of (4.5), we have
(4.8)∫0t‍(wn,pn-p0)ds≤α02∫0t‍|wn|2ds+12α0∫0t‍|pn-p0|2ds.
LetYn(t)=|wn|2+α0∫0t‍|wn|2ds. Substituting (4.6)–(4.8) into (4.5), then we have
(4.9)Yn(t)≤1α0∫0t‍|pn-p0|2ds+(α0+2λ)∫0t‍|wn|2ds.
Using Gronwall inequality, we obtain(4.10)Yn(t)≤1α0eα0T∫0t‍|pn-p0|2ds.
From all the above and∥un-u0∥→0, we get |pn-p0|→0, the mapping u→y(t,x;u) is thus strongly continuous.


## 4.3. Existence of Optimal Parameter

Theorem 4.2.
Suppose that assumptions (A1)–(A5) hold, then there exists at least one optimal parameteru*∈Uad(I¯T) satisfies the identification problem (SIP).Proof.
Obviously,f0(u) is continuous on Vu(t), and, from (4.3), J(u)=∫IT‍∫Ωu(t)‍f0(u)dxdt≥0. By Theorem 4.1, the mapping u→y(t,x;u) is continuous. Hence, the mapping u→J(u) is continuous on Uad(I¯T). Since Uad(I¯T) is a nonempty, bounded, and closed set, there exists u*∈Uad(I¯T), such that for all u∈Uad(I¯T), J(u*)≤J(u), e.g., u*∈Uad(I¯T) is an optimal parameter, thus the desired result is obtained.


## 4.4. Necessary Optimality Condition

Letu*∈Uad(I¯T) be the optimal parameter of the problem (SIP), following the convex compactness of Uad(I¯T) and the continuity of J(u) on Uad(I¯T), we can prove that J(u(·)) is Gateaux differentiable at u*∈Uad(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 letu*∈Uad(I¯T) be the optimal parameter of the system (SIP), then u*(·) satisfies the following inequality:
(4.11)DJ(u*;u-u*)=lim⁡γ→0+J(u*+γ(u-u*))-J(u*)γ≥0,∀u∈Uad(I¯T).


## 5. Optimization and Numerical Algorithm

### 5.1. 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
(5.1)ρici⋅∂Ti(x,t)∂t=∂∂x(Ki⋅∂Ti(x,t)∂x)+gi(x,t),(x,t)∈Ωi(t)×IT,i=1,2,3,
with the initial condition
(5.2)Ti(x,0)=T0(x),x∈Ωi(t),i=1,2,3,
and the boundary and penetrating conditions
(5.3)∂T1(0,t)∂t=Q1(t)T1(h1(t),t)=T2(h1(t),t),t∈Iρ1c1k1⋅∂T1(x,t)∂x|x=h1(t)=ρ2c2k2⋅∂T2(x,t)∂x|x=h1(t)(5.4)T2(h1(t),t)=T1(h1(t),t)T2(h1(t)+h2(t),t)=T3(L,t),t∈Iρ2c2k2⋅∂T2(x,t)∂x|x=h1(t)+h2(t)-ρ3c3k3⋅∂T3(x,t)∂x|x=h1(t)+h2(t)=ρ2dh2(t)dt(5.5)T3(h1(t)+h2(t),t)=T2(h1(t)+h2(t),t)T3(L,t)=T3(t)t∈I.Denote the subsystem composed by (5.1), (5.2), and (5.3) as (IBP)1, by (5.1), (5.2), and (5.4) as (IBP)2, and (5.1), (5.2), and (5.5) as (IBP)3.According to the partial differential theory, we can obtain the following theorem.Theorem 5.1.
Suppose that assumptions (A1)–(A5) are valid, then, for anyu∈Uad(I¯T), the subsystem (IBP)i has a unique weak solution Ti(x,t;u)∈C(Ωi(t),IT,Uad;R) which depends continuously on u∈Uad(I¯T).LetM, N be the numbers of the measured spatial and temporal spots, respectively, xk the measured depth, tj the measured time, T(xk,tj;u) the calculated temperature of sea ice at depth xk and time tj from the system (2.6), T¯(xk,tj) the measured temperature data at depth xk and time tj, then the practical identification problem can be expressed as
(SIPOD)min⁡Jd(u)=∑k=1M‍∑j=1N‍(T(xk,tj;u)-T¯(xk,tj))2s.t.T(xk,tj;u)∈SUad′(QT)u∈Uad(I¯T).LetSiUad′(QT)={Ti(x,t;u)∣Ti(x,t;u) is the solution of the subsystem (IBP)i corresponding to u∈Uad(I¯T)},i=1,2,3, and for the numerical implementation, the problem (SIPOD) can be decomposed to the following three connected subproblems
(SIPODi)min⁡Jdi(u)=∑zk∈Ωi‍∑j=1N(T(xk,tj;u)-T¯(xk,tj))2‍s.t.Ti(xk,tj;u)∈SiUad′(QT),i=1,2,3u∈Uad(I¯T).Obviously, the performance function of (SIPOD) can be described as(5.6)Jd(u)=∑i=13‍Jdi(u).From Theorem4.2, we can get that there also exists at least one optimal parameter satisfying the subproblem (SIPODi).Set(5.7)TΩi(u)=∑xk∈Ωi‍∑j=1N‍T(xk,tj;u),i=1,2,3,TΩi=∑xk∈Ωi‍∑j=1N‍T¯(xk,tj),i=1,2,3,Jbi(u)=(TΩi(u)-TΩi)2,i=1,2,3,Jb(u)=∑i=13‍Jbi(u),
then, the identification problem is expressed as
(SIPODb)min⁡Jb(u)s.t.T(xk,tj;u)∈SUad′(QT)u∈Uad(I¯T),
and the subidentification problem is
(SIPODbi)min⁡Jbi(u)s.t.T(xk,tj;u)∈SiUad′(QT),i=1,2,3u∈Uad(I¯T).Theorem 5.2.
Ifu* is the optimal parameter of SIPOD, then u* is also that of (SIPODb).


### 5.2. 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(zk,tj)}, a numerical optimization algorithm is constructed with the following steps.Step 1.
Select starting pointsu0=(u10,u20)∈Uad(I¯T), n directions e1, e2,…,en, starting step length Δu=(Δu1,Δu2), Δu1>0, Δu2>0, acceleration factor α, accuracy ɛ>0, and the maximum iteration number kmax, set v0=u0, k=j=0;Step 2.
Calculate the numerical solutionT(x,t;vj), by nonoverlapping Schwarz alternating direction method and semi-implicit difference scheme.Step 3.
IfJb(vj+Δu·ej)<Jb(vj), set vj+1=vj+Δu·ej; for the case Jb(vj+Δu·ej)≥Jb(vj), if Jb(vj-Δu·ej)<Jb(vj), set vj+1=vj-Δu·ej, else, set vj+1=vj.Step 4.
Ifj<n, set j=j+1, go to Step 3; for the case j=n, if Jb(vn+1)≤Jb(uk), go to Step 5, else, go to Step 6.Step 5.
Setuk+1=vn+1, v0=uk+1+α(uk+1-uk), k=k+1, and j=0, go to Step 2.Step 6.
Ifk≥kmax, stop, set u*=uk, else, set v0=uk, uk+1=uk, k=k+1, and j=0, go to Step 2.


## 5.1. 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
(5.1)ρici⋅∂Ti(x,t)∂t=∂∂x(Ki⋅∂Ti(x,t)∂x)+gi(x,t),(x,t)∈Ωi(t)×IT,i=1,2,3,
with the initial condition
(5.2)Ti(x,0)=T0(x),x∈Ωi(t),i=1,2,3,
and the boundary and penetrating conditions
(5.3)∂T1(0,t)∂t=Q1(t)T1(h1(t),t)=T2(h1(t),t),t∈Iρ1c1k1⋅∂T1(x,t)∂x|x=h1(t)=ρ2c2k2⋅∂T2(x,t)∂x|x=h1(t)(5.4)T2(h1(t),t)=T1(h1(t),t)T2(h1(t)+h2(t),t)=T3(L,t),t∈Iρ2c2k2⋅∂T2(x,t)∂x|x=h1(t)+h2(t)-ρ3c3k3⋅∂T3(x,t)∂x|x=h1(t)+h2(t)=ρ2dh2(t)dt(5.5)T3(h1(t)+h2(t),t)=T2(h1(t)+h2(t),t)T3(L,t)=T3(t)t∈I.Denote the subsystem composed by (5.1), (5.2), and (5.3) as (IBP)1, by (5.1), (5.2), and (5.4) as (IBP)2, and (5.1), (5.2), and (5.5) as (IBP)3.According to the partial differential theory, we can obtain the following theorem.Theorem 5.1.
Suppose that assumptions (A1)–(A5) are valid, then, for anyu∈Uad(I¯T), the subsystem (IBP)i has a unique weak solution Ti(x,t;u)∈C(Ωi(t),IT,Uad;R) which depends continuously on u∈Uad(I¯T).LetM, N be the numbers of the measured spatial and temporal spots, respectively, xk the measured depth, tj the measured time, T(xk,tj;u) the calculated temperature of sea ice at depth xk and time tj from the system (2.6), T¯(xk,tj) the measured temperature data at depth xk and time tj, then the practical identification problem can be expressed as
(SIPOD)min⁡Jd(u)=∑k=1M‍∑j=1N‍(T(xk,tj;u)-T¯(xk,tj))2s.t.T(xk,tj;u)∈SUad′(QT)u∈Uad(I¯T).LetSiUad′(QT)={Ti(x,t;u)∣Ti(x,t;u) is the solution of the subsystem (IBP)i corresponding to u∈Uad(I¯T)},i=1,2,3, and for the numerical implementation, the problem (SIPOD) can be decomposed to the following three connected subproblems
(SIPODi)min⁡Jdi(u)=∑zk∈Ωi‍∑j=1N(T(xk,tj;u)-T¯(xk,tj))2‍s.t.Ti(xk,tj;u)∈SiUad′(QT),i=1,2,3u∈Uad(I¯T).Obviously, the performance function of (SIPOD) can be described as(5.6)Jd(u)=∑i=13‍Jdi(u).From Theorem4.2, we can get that there also exists at least one optimal parameter satisfying the subproblem (SIPODi).Set(5.7)TΩi(u)=∑xk∈Ωi‍∑j=1N‍T(xk,tj;u),i=1,2,3,TΩi=∑xk∈Ωi‍∑j=1N‍T¯(xk,tj),i=1,2,3,Jbi(u)=(TΩi(u)-TΩi)2,i=1,2,3,Jb(u)=∑i=13‍Jbi(u),
then, the identification problem is expressed as
(SIPODb)min⁡Jb(u)s.t.T(xk,tj;u)∈SUad′(QT)u∈Uad(I¯T),
and the subidentification problem is
(SIPODbi)min⁡Jbi(u)s.t.T(xk,tj;u)∈SiUad′(QT),i=1,2,3u∈Uad(I¯T).Theorem 5.2.
Ifu* is the optimal parameter of SIPOD, then u* is also that of (SIPODb).


## 5.2. 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(zk,tj)}, a numerical optimization algorithm is constructed with the following steps.Step 1.
Select starting pointsu0=(u10,u20)∈Uad(I¯T), n directions e1, e2,…,en, starting step length Δu=(Δu1,Δu2), Δu1>0, Δu2>0, acceleration factor α, accuracy ɛ>0, and the maximum iteration number kmax, set v0=u0, k=j=0;Step 2.
Calculate the numerical solutionT(x,t;vj), by nonoverlapping Schwarz alternating direction method and semi-implicit difference scheme.Step 3.
IfJb(vj+Δu·ej)<Jb(vj), set vj+1=vj+Δu·ej; for the case Jb(vj+Δu·ej)≥Jb(vj), if Jb(vj-Δu·ej)<Jb(vj), set vj+1=vj-Δu·ej, else, set vj+1=vj.Step 4.
Ifj<n, set j=j+1, go to Step 3; for the case j=n, if Jb(vn+1)≤Jb(uk), go to Step 5, else, go to Step 6.Step 5.
Setuk+1=vn+1, v0=uk+1+α(uk+1-uk), k=k+1, and j=0, go to Step 2.Step 6.
Ifk≥kmax, stop, set u*=uk, else, set v0=uk, uk+1=uk, k=k+1, and j=0, go to Step 2.


## 6. Numerical Results

Using the above optimal algorithm, we can obtain the temperature distributionT(xk,tj;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 Figure2, 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.Figure 2
Comparison of calculated and measured temperatures.


## 7. 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 (Figure2) 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.



---

*Source: 532870-2012-07-12.xml*