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
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.
A thermodynamic model for the Arctic sea ice was firstly developed in  (hereafter MU71), in which the authors investigated the correspondence between the annual cycle pattern acquired from the simulation and the empirical data of . 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”  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,  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  by taking into account the internal brine pocket melting on surface ablation and the vertically varying salinity profile. Their thermodynamic model was demonstrated by  using a global climate model with a Lagrangian ice thickness distribution. Combining these two models,  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 , 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  through a satellite called ”ICESat” during 2003-2008, and compared with the data in  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,  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 . In , 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 , 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  and estimators with boundary measurements in [9, 10]. As further extensions,  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  investigated an input-to-state stability of the control of Stefan problem in  with respect to an unknown heat loss at the interface.
Our conference paper  applied the state estimation for the Stefan problem developed in  to the thermodynamic model of Arctic sea ice in . The present paper improves the design of  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.
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 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.
The sea ice thickness is positive and upper bounded, i.e. there exists such that
The velocity of the thickness is bounded, i.e., there exists a positive constant such that
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.
where . Next, we define the estimation error states as
where the negative sign is added to be consistent with the model developed in  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.
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.
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:
where is to be determined. Taking the first and second spatial derivatives of the transformation (3.2), we get
Taking integration by parts and substituting the boundary condition (38), we get
Substituting and into (3.2), we get
Moreover, substituting into (3.2) yields
and the observer gains must satisfy
and the function must satisfy
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
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
Next, we consider
Let be the Lyapunov functional defined by
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
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.
4.1 Input Parameters
The input parameters are taken from  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 , 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 . 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 .
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