Arctic Sea Ice State Estimation From Thermodynamic PDE Model

Arctic Sea Ice State Estimation
From Thermodynamic PDE Model

[    [

Recent rapid loss of the Arctic sea ice motivates the study of the Arctic sea ice thickness. Global climate model that describes the ice’s thickness evolution requires an accurate spatial temperature profile of the Arctic sea ice. However, measuring the complete temperature profile is not feasible within and throughout the Arctic icecap. Instead, measuring the ice’s thickness is doable with the acquisition of data from submarine and satellite devices. In this paper, we develop a backstepping observer algorithm to estimate the temperature profile for the Arctic sea ice model via available measurements of sea ice thickness and sea ice surface temperature. The observer is designed in a rigorous manner to drive the temperature profile estimation error to zero, for a salinity-free sea ice model. Moreover, the proposed observer is used to estimate the temperature profile of the original sea ice model with salinity via numerical simulation. In comparison with the straightforward open-loop algorithm, the simulation results illustrate that our observer design achieves ten times faster convergence of the estimated temperature.

UCSD]Shumon Koga, UCSD]Miroslav Krstic

Department of Mechanical and Aerospace Engineering, University of California, San Diego, La Jolla, CA 92093-0411 USA 

Key words:  Arctic sea ice, State estimation, Moving boundaries, Backstepping


1 Introduction

1.1 Motivation

The Arctic sea ice has been studied intensively in the field of climate and geoscience. One of the main reasons is due to ice-albedo feedback which influences climate dynamics through the high reflectivity of sea ice. The other reason is the rapid decline of the Arctic sea ice extent in the recent decade shown in several observations. These observations motivate the investigation of future sea ice amount. Several studies have developed a computational model of the Arctic sea ice and performed numerical simulations of the model with initial sea ice temperature profile. However, the spatially distributed temperature in sea ice is difficult to recover in realtime using a limited number of thermal sensors. Hence, the online estimation of the sea ice temperature profile based on some available measurements is crucial for the prediction of the sea ice thickness.

1.2 Literature

A thermodynamic model for the Arctic sea ice was firstly developed in [23] (hereafter MU71), in which the authors investigated the correspondence between the annual cycle pattern acquired from the simulation and the empirical data of [31]. The model involves a temperature diffusion equation evolving on a spatial domain defined as the sea ice thickness. Due to melting or freezing phenomena, the aforementioned spatial domain is time-varying. Such a model is called “Stefan problem” [6] which is described by a parabolic partial differential equation (PDE) with a state-dependent moving boundary driven by a Neumann boundary value.

Refined models of MU71 have been suggested in literature. For instance, [27] proposed a numerical model to achieve faster and accurate computation of MU71 by discretizing the temperature profile into some layers and neglecting salinity effect. The energy conserving model of MU71 was introduced in [2] by taking into account the internal brine pocket melting on surface ablation and the vertically varying salinity profile. Their thermodynamic model was demonstrated by [3] using a global climate model with a Lagrangian ice thickness distribution. Combining these two models, [32] developed an energy conserving three-layer model of sea ice by treating the upper half of the ice as a variable heat capacity layer.

Remote sensing techniques have been employed to obtain the Arctic sea ice data in several studies. In [7], authors suggested an algorithm to calculate sea ice surface temperature using the satellite measured brightness temperatures, which provided an excellent measurement of the actual surface temperature of the sea ice during the Arctic cold period. The Arctic sea ice thickness data were acquired in [20] through a satellite called ”ICESat” during 2003-2008, and compared with the data in [26] observed by submarine during 1958-2000. More recent data describing the evolution of the sea ice thickness have been collected between 2010 and 2014 from the satellite called ”CryoSat-2” [22, 21].

On the other hand, state estimation has been studied as a specific type of data assimilation which utilizes the numerical model along with measured value. For finite dimensional systems associated with noisy uncertain measurement, a systematic well known approach is Kalman Filter. Another well known method is Luenberger type state observer which reconstructs the state variable from partially measured variables. For the application to sea ice, [4] developed an adjoint based method as an iterative state and parameter estimation for the coupled sea ice-ocean in the Labrador Sea and Baffin Bay to minimize an uncertainty-weighted model-data misfit in a least square sense as suggested in [33, 34] using Massachusetts Institute of Technology general circulation model (MITgcm) developed in [24]. In [5], the same methodology was applied to reconstruct the global ocean and ice concentration. Their sea ice model is based on the zero-layer approximation of numerical model in [27], which is a crude model lacking an internal heat storage and promoting fast melting.

State estimation for infinite dimensional systems is much more challenging than the estimation of finite dimensional systems. The infinite dimensional systems are often modeled by means of PDEs which describes a time evolution of spatially distributed state variables. A systematic estimation method for PDEs was developed via ”backsteppping design” which employs the state observer structure for PDEs, see [16, 18, 29, 30, 28] for instance. Recently, the backstepping method for PDEs has been employed to the Stefan problem for the design of boundary controllers in [8] and estimators with boundary measurements in [9, 10]. As further extensions, [11] designed a state feedback control for the Stefan problem under the material’s convection, [12, 15] designed a delay compensated control for the Stefan problem, and [14] investigated an input-to-state stability of the control of Stefan problem in [8] with respect to an unknown heat loss at the interface.

1.3 Contributions

Our conference paper [13] applied the state estimation for the Stefan problem developed in [10] to the thermodynamic model of Arctic sea ice in [23]. The present paper improves the design of [13] by removing the availability of the temperature gradient at the ice-ocean interface which is nearly impossible to measure. The designed algorithm provides the temperature profile estimation in Arctic sea ice via backstepping method. The observer utilizes the available measurements of the sea ice thickness and the ice surface temperature. A simplified model for the Arctic sea ice is formulated by neglecting the salinity effect and the exponential convergence of its temperature estimation is ensured by using the designed observer. The simulation study for the proposed observer along with the original model is performed in order to investigate the convergence performance numerically. The simulation results illustrate that the estimated temperature profile converges uniformly to the actual sea ice temperature ten times faster than the straightforward open-loop estimator.

1.4 Organization

This paper is organized as follows: Section 2 describes the thermodynamic model of the Arctic sea ice in MU71 and introduces its simplification. Section 3 develops the state estimation design for the simplified model via backstepping PDE observer and shows the exponential stability of the estimation error system. Section 4 illustrates the simulation of the designed observer applied to the original model of MU71. The paper ends with the concluding remarks and future direction in Section 5.

2 Thermodynamic Model of Sea Ice

The thermodynamic model of MU71 describes the time evolution of sea ice temperature profile in the vertical axis along its thickness which also evolves in time due to accumulation or ablation caused by energy balance.

2.1 Snow Covered Season Model

Fig. 1: Schematic of the vertical model of the Arctic sea ice.

Fig. 1 provides a schematic of the Arctic sea ice model. During the seasons other than the summer (July and August), the sea ice is covered by snow, and the surface position of the snow also evolves in time. Let , denote the temperature profile of snow and sea ice, and and denote the thickness of snow and sea ice. The total incoming heat flux from the atmosphere is denoted by , and the heat flux from the ocean is denoted by . The Arctic sea ice model suggested by MU71 gives the governing equations of a Stefan-type free boundary problem formulated as


where , , , , , , , , and are solar radiation penetrating the ice, Stefan-Bolzman constant, thermal conductivity of snow, density of snow, heat capacity of pure ice, thermal conductivity of pure ice, density of pure ice, melting point of surface snow, and melting point of bottom sea ice. The total heat flux from the air includes the following terms


where , , , , and denote the incoming solar short-wave radiation, the long-wave radiation from the atmosphere and clouds, the flux of sensible heat, the latent heat in the adjacent air, and the surface albedo, respectively. The heat capacity and thermal conductivity of the sea ice are affected by the salinity as


where denotes the salinity in the sea ice. and represent the weight parameters. The thermodynamic model (2.1)-(7) allows us to predict the future thickness of snow and sea ice and temperature profile given the accurate initial temperature profile and thickness. However, from the practical point of view, it is not feasible to obtain the complete temperature profile due to a limited number of thermal sensors. To deal with the problem, the estimation algorithm is designed so that the state estimation converges to the actual state starting from an initial estimate.

2.2 Simplification of the Model

Before considering the state estimation design, first we impose a simplification on the system (2.1)-(7). The effect of the salinity profile on the physical parameters is assumed to be sufficiently small so that it can be negligible, i.e.


Therefore, the heat equation of the sea ice temperature (2.1) is rewritten as


where the diffusion coefficient is defined as . Next, we impose the following assumptions.

Assumption 1

The sea ice thickness is positive and upper bounded, i.e. there exists such that

Assumption 2

The velocity of the thickness is bounded, i.e., there exists a positive constant such that


We formulate the observer structure for sea ice temperature estimation based on the simplified sea ice model composed of (12) and (6).

3 State Estimation Design

In this section, we derive the estimation algorithm utilizing some available measurements and show the exponential convergence of the designed estimation to the simplified sea ice temperature. The ice thickness and surface temperature are measured in several studies [7, 20, 22, 26].

3.1 Observer Structure

Suppose that the sea ice thickness and the ice surface temperature are obtained as available measurements and , i.e.


The state estimate of the sea ice temperature is governed by a copy of the plant (12) and (6)-(7) plus the error injection of , namely, as follows:


where . Next, we define the estimation error states as


where the negative sign is added to be consistent with the model developed in [10] for the liquid phase. Subtraction of the observer system (3.1)-(3.1) from the system (12) and (6)-(7) yields the closed-system of estimation error as


Our goal is to design the observer gains , , , so that the temperature error converges to zero, which enables the state estimate to track the actual sea ice temperature . The main theorem of this paper is stated as follows.

Theorem 3

Let Assumptions 1 and 2 hold. Consider the estimation error system (23)-(26) with the design of the observer gains


where , , and are positive free parameters, is defined by


and denotes the modified Bessel function of the -th kind. Then, there exist positive constants and such that, for all , the norm


satisfies the following exponential decay


namely, the origin of the estimation error system is exponentially stable in the spatial norm.

Remark 1

The observer gains (3)-(30) include the thickness , so the gains are not precomputed offline, but are easily calculated online, along with the state estimation.

Remark 2

The estimation of the snow temperature profile is also achievable using the same observer structure as in Theorem 3 with measurements of the snow thickness and snow surface temperature. To avoid the lengthy statement, we do not provide it in this paper.

3.2 Gain Derivation via State Transformation

The backstepping is a well-known method to design the observer gains for PDEs [16, 17, 18, 19]. Hereafter, the partial derivative of a variable in and are denoted as the variable with a subscript of and , respectively. For the estimation error system (23)–(26), we apply the following invertible transformations:


which map the estimation error system (23)-(26) into the following target system:


where is to be determined. Taking the first and second spatial derivatives of the transformation (3.2), we get


Next, taking the time derivative of (3.2) along the solution of the target system (3.2)–(39), we get


Taking integration by parts and substituting the boundary condition (38), we get


Thus, by (3.2) and (43), we have


Substituting and into (3.2), we get


Moreover, substituting into (3.2) yields


Therefore, for the equations (23)–(26) to hold, the gain kernel functions must satisfy the following conditions:


and the observer gains must satisfy


and the function must satisfy


The solutions to (48)–(50) are uniquely given by


Then, by taking the derivatives of the obtained gain kernels, we get


where is defined by (31). Using (56)–(59), the conditions (3.2)–(54) are led to the explicit formulations of the observer gains given as (3)–(30). In the similar manner, the conditions for the gain kernel functions of the inverse transformation (3.2) are given by


and, the function is obtained by


The solutions to (60)–(62) are given by


Using the solutions (64) and (65), the function is obtained explicitly by (63). Then, the solution (63) also satisfy the condition (3.2). Hence, the transformation from to is invertible.

3.3 Stability Analysis

In this section, we prove the exponential stability of the origin of the the estimation error system (23)-(26) in the spatial norm. First, we show the exponential stability of the origin of the target system (3.2)-(39). We consider the following Lyapunov functional


where denotes the spatial norm of , defined by . Taking the time derivative of (66) together with the solution of (3.2)-(39) yields


Applying Young’s and Cauchy-Schwarz inequalities to the last term in (67) with the help of Assumption 2, we have


where . Applying (68) to (67), we obtain the following inequality:


Next, we consider


Taking the time derivative of (70) along with (39) leads to


Let be the Lyapunov functional defined by


Combining (69) and (71), the time derivative of (72) is shown to satisfy the following inequality:


Hence, by choosing the gain parameter to satisfy


the inequality (73) is led to


Applying comparison principle to the differential inequality (75), we get


Hence, the target system (3.2)-(39) is exponentially stable at the origin in the spatial norm. Due to the invertibility of the transformations (3.2) and (3.2), there exist positive constants and such that the following inequalities hold


Hence, by (76) and (3.3), we arrive at


which completes the proof of Theorem 3. Note that the designed backstepping observer achieves faster convergence with a possibility of causing overshoot since the overshoot coefficient is a monotonically increasing function in the observer gains’ parameters .

4 Numerical Simulation

While we have focused on the simplified PDE (12) to derive a rigorous proof of the proposed state estimation design (3.1)-(3.1) with observer gains given by (3)–(30), simulation studies are performed by applying the estimation design to the original thermodynamic model (2.1)-(7) including salinity.

(a) Dynamic behavior of snow and sea ice.
(b) Time evolution of temperature profile in sea ice
Fig. 2: Simulation tests of the plant (2.1)–(7) with input parameters on annual cycle. Both (a) and (b) are in good agreement with the simulation results in [23].

4.1 Input Parameters

Unit W/m W/m W/m W/m
Jan. 0 168 19.0 0
Feb. 0 166 12.3 -0.323
Mar. 30.7 166 11.6 -0.484 0.83
Apr. 160 187 4.68 -1.45 0.81
May. 286 244 -7.26 -7.43 0.82
Jun. 310 291 -6.30 -11.3 0.78
Jul. 220 308 -4.84 -10.3 0.64
Aug. 145 302 -6.46 -10.7 0.69
Sep. 59.7 266 -2.74 -6.30 0.84
Oct. 6.46 224 1.61 -3.07 0.85
Nov. 0 181 9.04 -0.161
Dec. 0 176 12.8 -0.161
Table 1: Average monthly values for the energy fluxes from the atmosphere
Symbol Unit Value
kg/m 330
W/m/C 0.31
kg/m 917
J/kg/C 2110
W/m/C 2.034
kJ C/kg 18.0
W/m 0.117
W/m 1.59
/m 1.5
C -0.1
C -1.8
Table 2: Physical parameters of snow and sea ice

The input parameters are taken from [23] in SI units and Table 1 shows the monthly averaged values of heat fluxes coming from the atmosphere for each months. Table 2 shows the physical parameters of snow and pure ice (without salinity), weight parameters and of salinity effect on heat capacity and thermal conductivity, radiative energy , and melting temperatures at surface and bottom . Following the energy conserving model in [2], the salinity profile is described by


where , , and .

4.2 Simulation Test of MU71

Using the given data, firstly the simulation of (2.1)-(7) is performed and showed in Fig. 2 to recover the evolution of and in the annual season as in [23]. The dynamic behavior of the snow surface and the bottom of sea ice are shown in Fig. 2 (a), and the time evolution of the temperature profile in sea ice is illustrated in Fig. 2 (b). It can be seen that both of Fig. 2 (a) and (b) have a good agreement with the simulation results shown in [23].

(a) Open-loop estimation, i.e., and aaaa for .
(b) The proposed estimation with the observer gains given in (3)–(30).
(c) The proposed estimation with larger value of aaaa than the value used in (b). The overshoot is observed.
(d) The proposed estimation with smaller value of than the value used in (b). The convergence speed gets slower.
Fig. 3: Simulation results of the plant (2.1)–(7) and the estimator (3.1)-(3.1) using parameters on January. The designed backstepping observer achieves faster convergence to the actual state than the straightforward open-loop estimation.
(a) Temperature profiles on January 1st.
(b) Temperature profiles on January 2nd.
(c) Temperature profiles on January 3rd.
(d) The time evolution of .
Fig. 4: Simulation result of the plant (2.1)–(7) and the estimator (3.1)-(3.1) with parameters used in Fig. 3 (b).

4.3 Simulation Results of State Estimation

The simulation results of temperature estimation computed by (3.1)-(3.1) along with the available measurements obtained by the online calculation of (2.1)-(7) are shown in Fig. 3. Here the initial temperature profile is set to be piecewise linear in space plus the sine wave perturbation, formulated as


where which is obtained by solving fourth order algebraic equation from (2.1) and the input data, and is set as [C]. The initial temperature estimation is supposed to be a quadratic function in space, given by