# Thermal Diffusivity Identification of Distributed Parameter Systems to Sea Ice

Authors: Liqiong Shi; Zhijun Li; Enmin Feng; Yila Bai; Yu Yang
Journal: Journal of Applied Mathematics
 (2013)

---

## Abstract

A method of optimal control is presented as a numerical tool for solving the sea ice heat transfer problem governed by a parabolic partial differential equation. Taken the deviation between the calculated ice temperature and the measurements as the performance criterion, an optimal control model of distributed parameter systems with specific constraints of thermal properties of sea ice was proposed to determine the thermal diffusivity of sea ice. Based on sea ice physical processes, the parameterization of the thermal diffusivity was derived through field data. The simulation results illustrated that the identified parameterization of the thermal diffusivity is reasonably effective in sea ice thermodynamics. The direct relation between the thermal diffusivity of sea ice and ice porosity is physically significant and can considerably reduce the computational errors. The successful application of this method also explained that the optimal control model of distributed parameter systems in conjunction with the engineering background has great potential in dealing with practical problems.


---

## Body

## 1. Introduction

In recent years, optimal control has been regarded as a powerful method with the ever-increasing demands in science, technology, and engineering [1, 2]. Many examples have an infinite number of degrees of freedom and are often modeled as infinite-dimensional distributed parameter systems (DPSs) such as thermal processes [3], spatially distributed chemical reactions [4], fluid processes [5], and flexible beams or plates [6]. Their states can vary both temporally and spatially. In the view of the complexity of these practical problems, distributed parameter systems, usually described by partial differential equations, have received a special attention and are widely used for modeling in a number of applications, relative to lumped parameter systems.Consequently, considerable literature on the control theory of distributed parameter systems has been reported since the 1960s [7–9], for example, the solution, optimality conditions, and the analysis of an optimal control problem of DPSs. However, the applications are still limited, because no general system has all specific characteristics of every process. In many cases, it is difficult to obtain an exact partial differential equation of the process only from the physical and chemical laws, and thus, the system identification is often employed to estimate the unknown DPSs from data [10]. In this paper, the method of system identification is employed to the heat transfer problem of sea ice. Well known thermodynamic sea ice processes, especially the freeze and melt cycle in the polar regions, play an important role in the global climate system via their influence on the horizontal extent and the thickness of sea ice. The thermodynamics of sea ice, controlled by the parabolic heat conduction equation, is a typical example of distributed parameter systems.In fact, the optimal control has been used to deal with ice problems [11–15]. Based on the sea ice temperature distribution, Lv et al. constructed an optimal control model with the thermodynamic parameters (density, specific heat, and thermal conductivity) of snow, ice, and ocean mixed layer as the control variable [11]; Fang et al. obtained the characteristics of the ice thicknesses and the ice temperatures by establishing a distributed parameter coupled system with the ice thickness as a movable inner boundary of the space domain [12]; Tan et al. presented an optimal model with the thicknesses of snow and sea ice as parametric variables [13]. They all emphasized theoretical research on sea ice physical parameters, referring to the thermodynamic coupled system, the existence of the optimal control, and the necessary conditions for optimality, which have provided an academic foundation for the simulation and forecast of sea ice thermodynamics. Furthermore, Bai et al. [14, 15] derived ice-water drag coefficients of an isolated ice floe and the thermal diffusivity of river ice with the method of parameter identification. Their work has realized an application of the optimal control of distributed parameter systems to ice parameters.The physical parameters are crucial to exactly understand the sea ice thermal behavior. Much previous research in sea ice thermodynamics focused on the thermal conductivity, the specific heat, the latent heat, and so on, through the three determining factors—ice temperature, ice salinity, and ice density [16, 17]. Sea ice porosity which is defined as the ratio of brine and gas volume to the total volume, is also controlled by sea ice temperature, salinity, and density together [18, 19]. With the deepening of research, sea ice mechanical properties became evaluated by using ice porosity, replacing ice temperature, salinity, and density [20–22], but the thermal properties of sea ice are seldom directly studied on the basis of ice porosity. Thus, we intend to examine the thermal diffusivity of sea ice through ice porosity to fully understand the thermal process in sea ice by the method of optimal control of DPSs.In this paper, we proceed from the one-dimensional heat conduction equation description of sea ice. Based on the field measurements from Antarctica in 2005 and 2006, an optimal control model of distributed parameter systems was put forward for identifying the thermal diffusivity of sea ice by taking thermal characteristics of sea ice and the related research results as specific constraints, and the deviation between the calculated ice temperature and the measurements was used as the performance criterion. The goal is to obtain a general parameterization of the thermal diffusivity and ice porosity for Antarctic sea ice. This identified parameterization of the thermal diffusivity of sea ice and ice porosity shows more direct relation between the thermal diffusivity and ice porosity. Meanwhile, the parameterization for the two parameters not only makes their physical meaning more clear, but also simplifies the complicated calculation by using ice thermal conductivity, specific heat capacity, and ice density.


## 2. Problem Analysis

Physically, the heat transfer process within ice can be described by the classical one-dimensional heat conduction equation [23]. Usually, a heat source term is added in the heat conduction equation as considering the incoming solar radiation. However, the heat source term could be ignored here, because we only discuss heat transfer in the middle ice layer, where the influence from the solar radiation is limited. The solar radiation is absent especially in polar winter. The distributed parameter systems with the initial condition and the boundary conditions can be given by
(1)∂∂tT(x,t)=λ(ν)∂2∂x2T(x,t),(x,t)∈Q,T(x,t)|t=0=T0(x),x∈[x1,x2],T(x,t)|x=x1=T1(t),t∈(0,tf],T(x,t)|x=x2=T2(t),t∈(0,tf],
where λ(ν) is the thermal diffusivity of sea ice, ν is ice porosity, T is sea ice temperature, x is the spatial variable with x∈[x1,x2], x1 and x2 are the upper and lower boundaries, t is the time variable with t∈(0,tf], and tf is the final time for measurement. Let Ω=[x1,x2], I=(0,tf], Q={(x,t)∈R2:x∈[x1,x2],andt∈(0,tf]}=Ω×I. The initial and the boundary conditions are all given continuous functions through interpolating and fitting the measured ice temperature. Because λ(ν), T0(x), T1(t), and T2(t) are all continuous and differentiable, satisfying T0(x1)=T1(0), T0(x2)=T2(0), the system (1) exists a unique solution.The thermal diffusivity of ice is defined as the ratio of thermal conductivity to the product of ice density and specific heat. In a sense, thermal diffusivity is the measure of thermal inertia [24]. It is the most directly observed thermal property, associated with the rate of ice temperature changes. The higher the thermal diffusivity, the faster the propagation of heat into the medium. In other word, substances with high thermal diffusivity rapidly adjust their temperature to that of their surroundings, because they conduct heat quickly in comparison to their volumetric heat capacity or thermal bulk. In the previous system, the thermal diffusivity of sea ice λ(ν) is an unknown function of ice porosity to be estimated by the system identification [10]. Therefore, it requires an optimal control model to determine the mathematical relationship λ(ν) between the thermal diffusivity and ice porosity, which not only describes the thermal behavior, but also represents the thermal properties of sea ice. According to the previous research results in sea ice thermodynamics and the limitation of the available field measurements, the following issues should be discussed.The formulae for the effective specific heat capacity and thermal conductivity of sea ice by Schwerdtfeger [17] are employed to calculate the thermal diffusivity of sea ice by the definition of thermal diffusivity and the in situ available data of measured ice temperature, salinity, and density. Figure 1 shows how the thermal diffusivity of sea ice varies with ice temperature, salinity, and density. According to the calculation, the minimum and the maximum of the thermal diffusivity for the study area in Antarctica are 0.58×10-7m2·s-1 and 10.48×10-7m2·s-1. And also the statistical relationship of the porosity by Cox and Weeks [25] is used to estimate the minimum and the maximum of ice porosity, which are 11‰ and 250‰, respectively. In order to derive the representative results, we restrict the variation of the thermal diffusivity and ice porosity by expanding 10% of their minimum and the maximum calculated through the measured data. Consequently, the lower limits and the upper limits of the thermal diffusivity and the corresponding ice porosity in our model could be represented as
(2)λ∈[λl,λu]=[0.52×10-7,11.53×10-7],ν∈[νl,νu]=[10,275],
where the units of the thermal diffusivity λ and ice porosity ν are m2·s-1 and ‰, respectively. The subscripts l and u denote the lower limit and the upper limit.Variation of the thermal diffusivity of sea ice with ice temperature and salinity (a) and with ice temperature and density (b), calculated by previous formulae.
(a)
(b)It is important to note that the thermal diffusivity of sea ice decreases with the increase of ice porosity, at least within the range of ice porosity from 0‰  to 275‰. Mathematically, the thermal diffusivity of sea ice is a monotone function of ice porosity within its range of 0‰  to 275‰, which can be expressed as(3)dλ(ν)dν<0,ν∈(0,275].Another consideration is the extreme case. As ice porosity approaches 0, sea ice becomes to be pure ice, and the thermal diffusivity should be the value of pure ice. This case should not be neglected, because it originates from the concept of ice porosity. This property of sea ice is represented as(4)limv→0λ(ν)=10.8×10-7,
where the units of the thermal diffusivity λ and ice porosity ν are m2·s-1 and ‰, respectively.For the form of the parameterization of thermal diffusivity of sea ice and ice porosity is difficult to decide, formulations from literatures which are used to estimate the relations between sea ice properties and ice porosity (Table1) are tested and modified. With these formulations, the optimal control problem for examining the parameterization of the thermal diffusivity varied with ice porosity is converted into the parameter identification of determining parameter u=(k,a,b). Let
(5)Uad≔{u=(k,a,b)∣k∈I7={1,2,…,7},a∈[al,au],b∈[bl,bu]}.
For allu∈Uad, in the system (1), λ(ν) could be expressed as λ(ν;u). Parameter k∈I7={1,2,…,7} denotes the number of seven different forms of the parameterizations in Table 1, which illustrates the trend of the thermal diffusivity with ice porosity. According to the discussion of the extreme previous case, parameter a should be 10.8×10-7m2·s-1, which represents the thermal diffusivity of pure ice, while ice porosity equals 0. And parameter b∈[bl,bu] decides the rate of change in the thermal diffusivity curve. The related background of sea ice thermodynamics is presented in [30].Table 1
Formulations of sea ice properties and ice porosity from literatures and modified.
Number
Formulation
Reference and origin
1
λ
(
ν
)
=
a
(
1
-
b
ν
)
-
1
[26], used for the relation of sea ice tensile strength and brine volume
2
λ
(
ν
)
=
a
(
1
-
b
ν
)
4
[27], used for the relation of the elastic modulus and the porosity
3
λ
(
ν
)
=
a
(
1
-
b
ν
)
-
2
[27], used for the relation of the effect stress and the porosity
4
λ
(
ν
)
=
a
·
e
b
ν
[28], used for the elastic modulus and the porosity
5
λ
(
ν
)
=
a
(
1
+
ν
)
b
[19], used for the relation of the mechanical indicator and the porosity
6
λ
(
ν
)
=
a
(
1
-
b
ν
)
2
[29], used for the relation of maximum intensity and the porosity
7
λ
(
ν
)
=
a
(
1
-
b
ν
)
-
1
Newly modified form based on references


## 3. Thermal Diffusivity Identification of Sea Ice

Optimal control is commonly used for uncertain and complex inverse problems. On the basis of the previous comprehensive analysis, and field measurements in Antarctica, including ice temperature, ice salinity, and ice density, we could propose a succinct optimal control model to determine the parameterization of the thermal diffusivity of sea ice [31]:
(6)minf(u)=∫I∫Ω(T(x,t;u)-T¯(x,t))2dxdts.t.T(x,t;u)∈SUadu∈Uadλl≤λ(ν;u)≤λu,v∈[νl,νu]∂λ(ν;u)∂ν<0,ν∈(0,275],limv→0λ(ν;u)=10.8×10-7,
where T(x,t;u) is the calculated ice temperature from (1) and T¯(x,t) is a fitted function of the measured ice temperature. For the system (1), λ(ν), T0(x), and T1(t), T2(t) are all continuous and differentiable and satisfy T0(x1)=T1(0),T0(x2)=T2(0); thus, the system (1) exists a unique solution. Let T(x,t;u) denotes the solution of the system (1), then for allu∈Uad,T(x,t;u) is continuous concerning u∈Uad. Let
(7)SUad:={system(1)correspondingtou∈UadT(x,t;u)∣T(x,t;u)isthesolutionofthesystem(1)correspondingtou∈Uad}.LetV=[al,au]×[bl,bu]⊂R2, then V is a compact set. According to the system (1), the mapping u∈I7×V→T(x,t;u)∈SUad is continuous; thus, we can get that SUad is a compact set. Obviously, f(u) is continuous, and the set determined by the constraints in model (6) is compact. We derive ∃u*∈I7×V, such that for all u∈I7×V,f(u*)≤f(u). In other words, there exists an optimal parameter for the identification problem. With respect to properties and the necessary conditions for optimality and so forth, we will not repeat the detailed course of the proofs because of the limited length, and similar theoretical research has been reported in the literature [11–13].After overall discussion of the constraints in the model, for any givenk∈I7, parameter a=10.8×10-7m2·s-1, parameter b could be estimated by the range and the monotonicity of the thermal diffusivity.In the process of calculation, a simple numerical method—semi-implicit difference scheme—is applied to discretize the parabolic differential equation, for the convenience of calculation, and its good performance. This method is eplicit and unconditionally stable. LetΔt and Δx be time step and space step, respectively, nN=tf/Δt and iN=(x2-x1)/Δx denote time mesh points number and space mesh points number. We can get(8)∂T(x,t)∂t|in+(1/2)=λ(ν;u)∂2T(x,t)∂x2|in+(1/2),
(9a)Tin+1-TinΔt+o(Δt2)=λ(ν;u)(Ti+1n-TinΔx2-Tin+1-Ti-1n+1Δx2)=+o(Δx2+ΔtΔx+Δx2),i=1,2,…,iN-1,n=1,3,5,…,NN,(9b)Tin+1-TinΔt+o(Δt2)=λ(ν;u)(Ti+1n+1-Tin+1Δx2-Tin-Ti-1nΔx2)=i+o(Δx2+ΔtΔx+Δx2),i=iN-1,iN-2,…,1,n=2,4,6,…,NN-1.Therefore, the solution to the system (1) is simplified by(10a)Tin+1=λ(ν;u)Δt(Ti+1n+Ti-1n+1)+(Δx2-λ(ν;u)Δt)TinΔx2+λ(ν;u)Δti=1,2,…,iN-1,n=1,3,5,…,NN,(10b)Tin+1=λ(ν;u)Δt(Ti+1n+1+Ti-1n)+(Δx2-λ(ν;u)Δt)TinΔx2+λ(ν;u)Δti=iN-1,iN-2,…,1,n=2,4,6,…NN-1,where NN={nN,ifnN/2∉ZnN-1,otherwise.Considering the reality of the problem being solved, the simple genetic algorithm (GA) is employed in this paper. Genetic algorithms belong to the evolutionary algorithms, which generate solutions to optimization problems using techniques inspired by natural evolution and serve as an efficient option in science and engineering for solving practical problems [32]. Identifying the optimal settings for crossover and mutation probabilities is very important, because the variations of crossover and mutation operators influence the convergence and performance of the algorithm significantly [33]. Typical crossover probability and mutation probabilities are in the ranges of 0.5–1.0 and 0.005–0.05, respectively [34]. Therefore, we tried to test the different combinations of crossover probability and mutation probabilities in the typical ranges to determine the optimal crossover and mutation probabilities for seven formulations in Table 1 for the sake of improving the convergence performances of GA. The basic procedures are simply given as follows:Step 1.
Set the generation counterg←0, the maximum generation G for the evolution. Randomly generate an initial population P(0)={pi(u),(i=1,2,…,N)} including N individuals pi(u),(i=1,2,…,N).Step 2.
Calculate the thermal diffusivityλ(ν;u) for every individual pi(u) in P(k)={pi(u),(i=1,2,…,N)}, respectively, and judge their values. If λ(ν;u)∈[λl,λu] is satisfied, go to Step 4. Or remove the individual pi(u) and go to Step 3 to regenerate a new one.Step 3.
Generate a new individualpi(u) instead of the old one for the generation P(k) and go to Step 2. The goal for this is to make sure the thermal diffusivity of sea ice obtained in our study is in reasonable range.Step 4.
Evaluate the fitness valuef(u) of each generated individual by the objective function in the model (6). Apparently, the fitness value is the standard for judging the optimal solution.Step 5.
Produce the next generationP(k+1) according to the prior generation P(k) based on the crossover operation with a prespecified crossover probability and the mutation operation with a prespecified mutation probability. Here, the crossover and mutation probabilities are decided through a series of repeat tests and calculations to improve the performance of the algorithm. And the probabilities are very different for the different field measurements and seven formulations in Table 1.Step 6.
Ifg≤G, let g←g+1 and go to Step 2. Or terminate the calculation, and the individual with the minimum of the fitness value is as the optimal solution for the identification problem.According to the previous procedures and the measurements from seven ice samples collected in Antarctica in 2006, we derived the parameterizations of thermal diffusivity and ice porosity. Because of the various sampling dates, the obtained relations have various degrees of distinction, which are showed in Figure2. The optimal parameterization of the thermal diffusivity and ice porosity for the Antarctic sea ice is fitted and expressed as
(11)λ(ν)=10.8×10-7(1+ν)-0.302,v∈(0,275],
where the unit of the thermal diffusivity of sea ice λ is m2·s-1, and the unit of ice porosity is ‰.Figure 2
The optimal parameterization of the thermal diffusivity and ice porosity by fitting identified data from different ice samples at Antarctic Zhongshan Station on different dates in 2006.In Figure2, for the seven ice samples, the difference of parameter b causes the slight differences of the identified parameterizations. The reason for this is that the ranges of ice depth selected during identification are very various, although the ice porosity has reflected the predominant factors influencing properties of sea ice: ice temperature, ice salinity, and ice density. Most often, the initial rapid freezing on the upper surface of the ice forms granular crystals, and then, subsequent growth in the lower part of ice develops columnar crystals, which could affect the heat transfer within ice. Thus, a simple analysis was done between the proportion of the length of the granular crystals to the total length pgc for the ice simples used for our identification, and the values of parameter b and their correlation can be seen in Figure 3. It is obvious that the parameter b is highly dependent on ice crystal types, which illustrates that ice crystal has an important influence on the thermal diffusivity of sea ice and should be regarded as a factor to evaluate the thermal diffusivity in later work.Figure 3
The relation between the parameterb and the proportion of the length of the granular crystals to the total length pgc for the seven ice samples from Antarctic Zhongshan Station on different dates in 2006.


## 4. Numerical Simulations

To ensure the effectiveness of the optimal control model of distributed parameter systems and the validity of the parameterization of the thermal diffusivity of sea ice in this paper, we simulated vertical ice temperature profiles using the general parameterization of the thermal diffusivity varied with ice porosity (11) and compared them with the measured ice temperature from Chinese Antarctic Zhongshan Station on seven different dates in different months in 2005. Figure 4(a) shows the calculated and measured vertical temperature profiles at 13:00 on different dates. The simulated ice temperature profiles seemed to agree well with the measurements. For comparison, we also applied Schwerdtfeger’s formulae and the heat conduction equation to simulate the vertical temperature profiles. Figure 4(b) shows the results.The profiles of the measured and simulated ice temperature for 7 ice samples at 13:00 at Antarctic Zhongshan Station in 2005, which are calculated by the parameterization of the thermal diffusivity varied with ice porosity (a) and Schwerdtfeger’s formulae (b).
(a)
(b)In this paper, we introduce the average error(AE) to quantify the simulated results, which is defined as:
(12)AE=∑t∈I∑x∈Ω|T(x,t)-T¯(x,t)|∑t∈I∑x∈Ω|T¯(x,t)|.From Figure4, we can conclude that the trends of ice temperature deviations are approximately the same. The AE between simulated and measured ice temperature is 1.41% (Figure 4(a)) and 1.71% (Figure 4(b)), respectively. This indicates that the identified parameterization of the thermal diffusivity of sea ice (11) is applicable in sea ice thermodynamics and that the optimal control method we discussed is also effective.


## 5. Conclusions

Optimal control theory of distributed parameter systems plays an important role in physics, mechanics, engineering, chemical industry, and other fields of practical activity, because of the complexity of the system processes. An optimal control model of distributed parameter systems was proposed in this paper. Instead of words, mathematical expressions were taken as specific constraints of the model by analyzing thermal properties of sea ice and summarizing the previous research results, which linked mathematics with the essence of practical problems. This model can not only depict the heat transfer process within sea ice, but can also fully capture the specific features of sea ice.The general parameterization of the thermal diffusivity of sea ice in Antarctica was identified by applying the optimal control method:λ(ν)=10.8×10-7(1+ν)-0.302,ν∈(0,275].Using the parameterization, the thermal diffusivity of sea ice could be understood more directly through ice porosity, and simply calculated by ice porosity rather than by thermal conductivity, specific heat capacity, and ice density. Furthermore, this parameterization of the thermal diffusivity has more reasonable physical significance. And the simple analysis of the difference of the parameterizations for the different ice samples shows that ice crystal types have an effect on the thermal diffusivity of sea ice and should be discussed further.The utility of the optimal control model has been demonstrated by the simulated results. The average errors from the identified parameterization of sea ice thermal diffusivity and Schwerdtfeger’s formulae revealed the advantages and growing potentials of the optimal control combined with the background of practical problems.However, this paper focused on the application of the optimal control method to deal with the thermal diffusivity of sea ice; further efforts should be made to expand the practical application of the method.



---

*Source: 760378-2013-09-09.xml*