Simulation and data assimilation in an idealized coupled atmosphere–ocean–sea ice floe model with cloud effects

Sea ice plays a crucial role in the climate system, particularly in the Marginal Ice Zone (MIZ), a transitional area consisting of fragmented ice between the open ocean and consolidated pack ice. As the MIZ expands, understanding its dynamics becomes essential for predicting climate change impacts. However, the role of clouds in these processes has been largely overlooked. This paper addresses that gap by developing an idealized coupled atmosphere–ocean–ice model incorporating cloud and precipitation effects, tackling both forward (simulation) and inverse (data assimilation) problems. Sea ice dynamics are modeled using the discrete-element method, which simulates floes driven by atmospheric and oceanic forces. The ocean is represented by a two-layer quasi-geostrophic (QG) model, capturing mesoscale eddies and ice–ocean drag. The atmosphere is modeled using a two-layer saturated precipitating QG system, accounting for variable evaporation over sea surfaces and ice. Cloud cover affects radiation, influencing ice melting. The idealized coupled modeling framework allows us to study the interactions between atmosphere, ocean, and sea ice floes. Specifically, it focuses on how clouds and precipitation affect energy balance, melting, and freezing processes. It also serves as a testbed for data assimilation, which allows the recovery of unobserved floe trajectories and ocean fields in cloud-induced uncertainties. Numerical results show that appropriate reduced-order models help improve data assimilation efficiency with partial observations, allowing the skillful inference of missing floe trajectories and lower atmospheric winds. These results imply the potential of integrating idealized models with data assimilation to improve our understanding of Arctic dynamics and predictions.

1 Introduction

Sea ice is a critical component of our climate system, serving both as a reflective shield that deflects solar radiation and as an insulator that regulates oceanic heat (Thomas, 2017; Weeks, 2010). By controlling the exchange of heat, moisture, and momentum between the ocean and atmosphere, sea ice plays a key role in shaping global weather patterns and broader climate systems. Therefore, understanding the dynamics of sea ice is essential for improving climate models and making accurate predictions about future climate conditions (Bigg, 2003; Meier et al., 2014; Gildor and Tziperman, 2001; Bhatt et al., 2014; Leppäranta, 2011; Maslowski et al., 2012; Thomson et al., 2018; Weeks and Ackley, 1986; Weiss, 2013).

The Marginal Ice Zone (MIZ), the transitional area between the open ocean and consolidated pack ice, is characterized by fragmented ice floes that interact dynamically with oceanic and atmospheric forces, playing a critical role in energy exchange, ocean circulation, and climate regulation (Dumont, 2022; Thomson et al., 2018). As the climate warms, the MIZ is expanding, increasing the fragmentation of ice and intensifying ocean–atmosphere interactions. This growth contributes to shifts in ocean currents, weather patterns, and the overall climate system. The expansion also accelerates the ice–albedo feedback, where more open water absorbs solar radiation, driving further ice melt. Given its growing influence, understanding and modeling the MIZ and its ice floe dynamics is essential for predicting climate change impacts and developing effective adaptation strategies (Manucharyan and Thompson, 2017; Strong and Rigor, 2013; Timmermans et al., 2018; Squire, 2020).

In Earth system models, sea ice is typically represented using continuum frameworks with viscous-plastic rheology (Hibler III, 1979; Hunke and Dukowicz, 1997; Tremblay and Mysak, 1997; Toyoda et al., 2019), which effectively captures large-scale dynamics but often has difficulties in accounting for brittle behavior and fine-scale fragmentation. These models work well at basin scales. However, they lack the resolution needed to simulate individual ice floes and their interactions. In contrast, the discrete-element method (DEM) focuses on individual ice floes in Lagrangian coordinates and allows a more detailed representation of local interactions between ice, ocean, and atmosphere (Cundall, 1988, 1979; Hart et al., 1988). DEM also reduces computational costs for simulating MIZ by eliminating the need for the advective transport schemes required in continuum models, while offering flexible spatial resolution, making it particularly useful for modeling the complex dynamics of the MIZ (Lindsay and Stern, 2004; Manucharyan and Montemuro, 2022; Damsgaard et al., 2018; Bouillon and Rampal, 2015; Rampal et al., 2016; Deng et al., 2024).

One aspect that has not received enough attention in previous studies is the effect of clouds on sea ice dynamics. On the one hand, clouds influence the thermodynamics of the atmosphere–ice system by modulating radiative fluxes and precipitation, which subsequently affect ice melting and growth processes (Shine et al., 1984; Huang et al., 2019; Liu et al., 2012; Huang et al., 2017; Morrison et al., 2019; Kay and Gettelman, 2009). Understanding these interactions is crucial for accurately predicting sea ice features under different climate scenarios. Including cloud–sea ice interactions in the modeling framework helps enhance our understanding of MIZ dynamics, particularly for studying the response of the DEM to atmospheric influences. On the other hand, significant challenges appear when observing ice floes in the presence of clouds. In situ measurements are often sparse (Brunette et al., 2022; Gabrielski et al., 2015; Hutchings et al., 2012; Itkin et al., 2017; Lei et al., 2020) and have limited spatiotemporal resolution (Cámara-Mor et al., 2010; Kwok, 2018). As a result, satellite imagery is widely used to monitor ice floe motion in the MIZ, which is then employed to infer ocean currents (Manucharyan et al., 2022; Lopez-Acosta et al., 2019; Chen et al., 2022b; Covington et al., 2022). However, ice floes can become obscured in satellite images due to intermittent cloud cover. Understanding the accuracy of recovering floe trajectories and the ocean fields with cloud cover facilitates the study of the MIZ. Recent advances in sea ice data assimilation encompass a broad spectrum of approaches. For example, Lisæter et al. (2003) and Massonnet et al. (2014) assimilated passive-microwave concentration and altimetry-derived thickness into coupled ice–ocean models with an ensemble Kalman filter, substantially reducing drift and thickness errors. Riedel and Anderson (2024) accounted for the bounded, non-Gaussian statistics of sea ice variables within the observation operator, which refines the posterior analyses of both ice and snow states. At the fully coupled level, Penny et al. (2019) introduced a strongly coupled data-assimilation (SCDA) framework that puts sea-surface and ice increments directly into the atmospheric analysis, further improving the short-term forecasts in the marginal-ice zone. With traditional Eulerian approaches, Chen et al. (2022a) and Deng et al. (2025) developed an efficient Lagrangian scheme that reconstructs mesoscale currents and vorticity from a limited set of tracked floes, even if only partial trajectories are observed due to clouds. Nevertheless, current data-assimilation frameworks for fully coupled atmosphere–ocean–ice models still lack a consistent treatment of cloud and precipitation effects.

This paper works toward filling these gaps by developing an idealized coupled atmosphere–ocean–ice model that incorporates the effects of clouds and precipitation. This model serves several important purposes. It allows the study of the fundamental physics governing interactions between the atmosphere, ocean, and ice floes, which gives a comprehensive understanding of how these components influence each other. It also provides an idealized modeling framework for analyzing the effect of clouds on ice floes. It specifically examines how cloud cover and precipitation modify the energy balance and influence the processes of ice melting and freezing. Moreover, it functions as a testbed for evaluating the accuracy of inferring missing observations of ice floe trajectories and the underlying ocean fields in the presence of clouds through data assimilation (DA). The primary goals of this study are thus twofold, addressing both forward (model simulation) and inverse (DA) problems: to develop and analyze the idealized coupled atmosphere–ocean–ice model that incorporates the effects of clouds and offers insights into the fundamental interactions within the MIZ, and to develop and test an efficient DA scheme using reduced-order models aimed at recovering unobserved variables and fields despite limited and uncertain observations.

The remainder of the paper is organized as follows: Sect. 2 presents the coupled atmosphere–ocean–ice model. Section 3 discusses the development of cheap surrogate forecast models for studying DA. Section 4 provides numerical simulation results that demonstrate the model dynamics, while Sect. 5 discusses the DA results. The paper is concluded in Sect. 6.

2 The idealized coupled atmosphere–ocean–sea ice floe system

2.1 Overview

To incorporate the effects of clouds and precipitation, we develop a coupled atmosphere–ocean–ice system that provides an understanding of the interactions among these components.

In this framework, sea ice dynamics are modeled using the DEM, which represents individual ice floes as circular elements with specific sizes and masses (Cundall, 1988; Hart et al., 1988). The DEM effectively captures the interactions and shape-preserving behaviors of floes, with their motion driven by atmospheric and oceanic forces. These forces are quantified through surface integrals over the floes.

Ocean dynamics are modeled using a two-layer quasi-geostrophic (QG) model, known as the Phillips model (Vallis, 2017; Salmon, 1998). This model effectively simulates eddies resulting from baroclinic instabilities, which are crucial for accurately representing oceanic conditions. The model is configured to reflect an Arctic Ocean regime (Qi and Majda, 2016), capturing key interactions between the ocean and the ice floes.

The atmospheric component employs a two-layer saturated precipitating quasi-geostrophic (PQG) model (Smith and Stechmann, 2017; Edwards et al., 2020a, b; Hu et al., 2021). This model addresses the atmospheric forces acting on the ice floes and incorporates the precipitation dynamics that contribute to their growth. The source of water vapor is represented by a parameter quantifying evaporation within the saturated PQG framework (Edwards et al., 2020a, b), varying between the sea surface and the ice floe surface. This variation allows the model to simulate how changes in ice floe distribution inversely impact atmospheric thermodynamics. In addition, radiation is incorporated to model the melting of ice floes, with the magnitude adjusted based on cloud thickness within the saturated PQG model.

In the following, Sect. 2.2 presents the DEM model; Sect. 2.3 outlines the two-layer QG model; Sect. 2.4 introduces and discusses the saturated PQG model; and Sect. 2.5 explores the thermodynamic processes across the different models. Figure 1a illustrates the coupled atmosphere–ice–ocean system along with its interacting components, while Fig. 1b provides a sectional overview of these components within the coupled atmosphere–ice–ocean system, illustrating their interactions.

Figure 1 (a) Illustration of the coupled atmosphere–ice–ocean system. (b) Sectional breakdown of the coupled atmosphere–ice–ocean system.

2.2 Ice floe dynamics: the DEM model

2.2.1 Governing equations

The DEM is employed to model the motion of sea ice floes, which are simplified as rigid circular bodies (Cundall, 1979, 1988; Hart et al., 1988; Chen et al., 2021). Each floe is defined by its position, ${\mathit{x}}^l=(x^l,y^l)$, and angular displacement, Ω^l, where $l=\mathrm{1},\mathrm{2},\mathrm{\dots },L$ indexes the individual floes. The governing equations for the motion of each floe are as follows:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}& m^l{\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}{\mathit{x}}^l}{\mathrm{d}{t}^{\mathrm{2}}}}=\underset{A}{\iint }{\mathbf{F}}^l\mathrm{d}A+{\mathbf{F}}_{\text{contact}}^l,\\ & I^l{\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}{\mathrm{\Omega }}^l}{\mathrm{d}{t}^{\mathrm{2}}}}=\underset{A}{\iint }{\mathcal{T}}^l\mathrm{d}A+{\mathcal{T}}_{\text{contact}}^l,\end{array}\end{array}$$

where the position x^l is defined at the center of mass of the lth floe. The second-order time derivative $\frac{{\mathrm{d}}^{\mathrm{2}}{\mathit{x}}^l}{\mathrm{d}{t}^{\mathrm{2}}}$ represents the acceleration of the floe, influenced by the contact force with other floes ${\mathbf{F}}_{\text{contact}}^l$ and the total external force F^l, integrated over the floe's area A. Similarly, the angular acceleration $\frac{{\mathrm{d}}^{\mathrm{2}}{\mathrm{\Omega }}^l}{\mathrm{d}{t}^{\mathrm{2}}}$ results from the torque due to contacts with other floes ${\mathcal{T}}_{\text{contact}}^l$ and the external torque 𝒯^l, also integrated over A. Here, t denotes time, m^l is the mass of floe l, and I^l is its moment of inertia. Periodic boundary conditions are imposed in both the x and y directions to ensure continuity in the domain.

By defining the velocity at the floe's center of mass as ${\mathit{u}}^l=(u^l,v^l)$ and the angular velocity as ω^l, the equations in Eq. (1) can be reformulated into a set of first-order differential equations:

$$\begin{array}{}\text{(2)}& {\displaystyle }\mathrm{d}{\mathit{x}}^l={\mathit{v}}^l\mathrm{d}t,\text{(3)}& {\displaystyle }\mathrm{d}{\mathbf{\Omega }}^l={\mathit{\omega }}^l\mathrm{d}t,\text{(4)}& {\displaystyle }\begin{array}{rl}\mathrm{d}{\mathit{v}}^l& ={\displaystyle \frac{\mathrm{1}}{m^l}}\left(\underset{\text{Contact forces}}{\underbrace{\sum _j\left({\mathit{f}}_{\mathrm{n}}^{lj}+{\mathit{f}}_{\mathrm{t}}^{lj}\right)}}+\underset{\text{Ocean drag force}}{\underbrace{{\mathcal{D}}_{\mathrm{o}}\left({\mathit{u}}_{\mathrm{o}}\left({\mathit{x}}^l\right)-{\mathit{v}}^l\right)}}\right.\\ & \left.+\underset{\text{Atmosphere drag force}}{\underbrace{{\mathcal{D}}_{\mathrm{a}}\left({\mathit{u}}_{\mathrm{a}}\left({\mathit{x}}^l\right)-{\mathit{v}}^l\right)}}\right)\mathrm{d}t,\end{array}\text{(5)}& {\displaystyle }\begin{array}{rl}\mathrm{d}{\mathit{\omega }}^l& ={\displaystyle \frac{\mathrm{1}}{I^l}}\left(\underset{\text{Contact torque}}{\underbrace{\sum _j\left(r^l{\mathit{n}}^{lj}\times {\mathit{f}}_t^{lj}\right)\cdot \widehat{\mathit{z}}}}\right.\\ & +\underset{\text{Ocean drag torque}}{\underbrace{{\mathcal{T}}_{\mathrm{o}}\left({\displaystyle \frac{\mathrm{\nabla }\times {\mathit{u}}_{\mathrm{o}}\left({\mathit{x}}^l\right)}{\mathrm{2}}}-{\mathit{\omega }}^l\widehat{\mathit{z}}\right)}}\\ & \left.+\underset{\text{Atmosphere drag torque}}{\underbrace{{\mathcal{T}}_{\mathrm{a}}\left({\displaystyle \frac{\mathrm{\nabla }\times {\mathit{u}}_{\mathrm{a}}\left({\mathit{x}}^l\right)}{\mathrm{2}}}-{\mathit{\omega }}^l\widehat{\mathit{z}}\right)}}\right)\mathrm{d}t.\end{array}\end{array}$$

The notations for the variables and parameters in Eqs. (2)–(5) are listed in Table 1.

Table 1 Notation in the DEM model.

The details of the ocean drag force and torque are covered in Sect. 2.2.2, while the atmospheric drag force and torque are detailed in Sect. 2.2.3. The details of the contact forces are provided in Appendix A. While the current coupling model does not include the Coriolis force on ice floe movement, it is still a valid assumption as this research focuses on short-term, localized ice floe dynamics where collision and drag forces are dominant (Thorndike and Colony, 1982; Steele et al., 1997; Thorndike, 1986).

2.2.2 Oceanic forcing

In this section we present the mathematical formulations for the drag force and torque exerted on an ice floe by ocean currents. These expressions are crucial for understanding how the ocean's movement affects the floe's velocity and angular velocity.

Drag force

The velocity of the lth ice floe, denoted as v^l, is affected by the ocean's drag force, which follows a quadratic drag law. The drag force exerted on the ice floe is expressed as

$$\begin{array}{}\text{(6)}& {\displaystyle }{\mathcal{D}}_{\mathrm{o}}^l\left({\mathit{u}}_{\mathrm{o}}\left({\mathit{x}}^l\right)-{\mathit{v}}^l\right)=\stackrel{\mathrm{̃}}{\mathit{\alpha }}({\mathit{u}}_{\mathrm{o}}-{\mathit{v}}^l)\left|{\mathit{u}}_{\mathrm{o}}-{\mathit{v}}^l\right|,\end{array}$$

where u_o(x^l) represents the sea surface current velocity at the centroid of the cylinder floe, i.e., x^l. The coefficient $\stackrel{\mathrm{̃}}{\mathit{\alpha }}$ is defined as

$$\begin{array}{}\text{(7)}& {\displaystyle }\stackrel{\mathrm{̃}}{\mathit{\alpha }}=d_{\mathrm{o}}{\mathit{\rho }}_{\mathrm{o}}\mathit{\pi }(r^l{)}^{\mathrm{2}},\end{array}$$

where d_o is the ocean drag coefficient, ρ_o is the density of ocean water, and r^l is the radius of the ice floe.

Drag torque

The torque due to ocean drag, denoted as 𝒯_o, influences the angular velocity ω^l of the ice floe. Assuming a quadratic drag law, the governing equation for the torque is expressed as

$$\begin{array}{}\text{(8)}& {\displaystyle }{\mathcal{T}}_{\mathrm{o}}\left({\displaystyle \frac{\mathrm{\nabla }\times {\mathit{u}}_{\mathrm{o}}}{\mathrm{2}}}-{\mathit{\omega }}^l\right)={\stackrel{\mathrm{̃}}{\mathit{\beta }}}_{\mathrm{o}}\left({\displaystyle \frac{\mathrm{\nabla }\times {\mathit{u}}_{\mathrm{o}}}{\mathrm{2}}}-{\mathit{\omega }}^l\right)\left|{\displaystyle \frac{\mathrm{\nabla }\times {\mathit{u}}_{\mathrm{o}}}{\mathrm{2}}}-{\mathit{\omega }}^l\right|.\end{array}$$

Here, $\frac{\mathrm{\nabla }\times {\mathit{u}}_{\mathrm{o}}}{\mathrm{2}}$ represents half of the curl of the ocean surface velocity, corresponding to the angular velocity of the ocean. The coefficient ${\stackrel{\mathrm{̃}}{\mathit{\beta }}}_{\mathrm{o}}$ is defined as

$$\begin{array}{}\text{(9)}& {\displaystyle }\stackrel{\mathrm{̃}}{\mathit{\beta }}=d_{\mathrm{o}}{\mathit{\rho }}_{\mathrm{o}}\mathit{\pi }(r^l{)}^{\mathrm{4}},\end{array}$$

where d_o is the ocean drag coefficient and ρ_o is the density of the ocean water, consistent with the parameters used in Eq. (7).

2.2.3 Atmospheric forcing

In this section we present the mathematical formulations for the drag force and torque exerted on an ice floe by atmospheric wind. These expressions are crucial for understanding how the atmosphere's movement affects the floe's velocity and angular velocity.

Drag force

The velocity of the ice floe, denoted by v^l, is influenced by the atmospheric drag force, which follows a quadratic drag law. The equation governing this relationship is

$$\begin{array}{}\text{(10)}& {\displaystyle }{\mathcal{D}}_{\mathrm{a}}^l\left({\mathit{u}}_{\mathrm{a}}\left({\mathit{x}}^l\right)-{\mathit{v}}^l\right)=\stackrel{\mathrm{̃}}{\mathit{\alpha }}({\mathit{u}}_{\mathrm{a}}-{\mathit{v}}^l)\left|{\mathit{u}}_{\mathrm{a}}-{\mathit{v}}^l\right|,\end{array}$$

where u_a represents the atmosphere near-surface wind velocity. The coefficient $\stackrel{\mathrm{̃}}{\mathit{\alpha }}$ is defined as

$$\begin{array}{}\text{(11)}& {\displaystyle }{\stackrel{\mathrm{̃}}{\mathit{\alpha }}}_{\mathrm{a}}=d_{\mathrm{a}}{\mathit{\rho }}_{\mathrm{a}}\mathit{\pi }(r^l{)}^{\mathrm{2}},\end{array}$$

where d_a is the atmosphere drag coefficient and ρ_a is the density of the air.

Drag torque

The governing equation for the angular velocity perturbed by the atmosphere, 𝒯_a, is a quadratic function, which yields

$$\begin{array}{}\text{(12)}& {\displaystyle }{\mathcal{T}}_{\mathrm{a}}\left({\displaystyle \frac{\mathrm{\nabla }\times {\mathit{u}}_{\mathrm{a}}}{\mathrm{2}}}-{\mathit{\omega }}^l\right)=\stackrel{\mathrm{̃}}{\mathit{\beta }}\left({\displaystyle \frac{\mathrm{\nabla }\times {\mathit{u}}_{\mathrm{a}}}{\mathrm{2}}}-{\mathit{\omega }}^l\right)\left|{\displaystyle \frac{\mathrm{\nabla }\times {\mathit{u}}_{\mathrm{a}}}{\mathrm{2}}}-{\mathit{\omega }}^l\right|,\end{array}$$

where ${\stackrel{\mathrm{̃}}{\mathit{\beta }}}_{\mathrm{a}}$ is the coefficient of drag, defined as

$$\begin{array}{}\text{(13)}& {\displaystyle }{\stackrel{\mathrm{̃}}{\mathit{\beta }}}_{\mathrm{a}}=d_{\mathrm{a}}{\mathit{\rho }}_{\mathrm{a}}\mathit{\pi }(r^l{)}^{\mathrm{4}}.\end{array}$$

Here, d_a is the atmosphere drag coefficient and ρ_a is the density of the air, consistent with the parameters used in Eq. (11).

2.3 Ocean dynamics: a two-layer quasi-geostrophic model

The two-layer QG model operates in a rotating reference frame, with two layers of equal depth bounded by a rigid lid at the top and a flat bottom. The governing equations of this model are expressed in terms of barotropic and baroclinic modes for potential vorticity (PV) anomalies, with periodic boundary conditions imposed in both the x and y directions (Qi and Majda, 2016; Vallis, 2017; Salmon, 1998). The model is governed by the following equations:

$$\begin{array}{}\text{(14)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial q_{\mathrm{1}}^{\mathrm{o}}}{\partial t}}& +J({\mathit{\psi }}_{\mathrm{1}}^{\mathrm{o}},q_{\mathrm{1}})+\mathit{\beta }{\displaystyle \frac{\partial {\mathit{\psi }}_{\mathrm{1}}^{\mathrm{o}}}{\partial x}}+U_{\mathrm{o}}{\displaystyle \frac{\partial }{\partial x}}\left(\mathrm{\Delta }{\mathit{\psi }}_{\mathrm{1}}^{\mathrm{o}}+(k_{\mathrm{d}}^{\mathrm{o}}{)}^{\mathrm{2}}{\mathit{\psi }}_{\mathrm{1}}^{\mathrm{o}}\right)\\ & =-{\mathit{\nu }}_{\mathrm{o}}{\mathrm{\Delta }}^{\mathrm{4}}{q}_{\mathrm{1}}^{\mathrm{o}},\end{array}\text{(15)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial q_{\mathrm{2}}^{\mathrm{o}}}{\partial t}}& +J({\mathit{\psi }}_{\mathrm{2}}^{\mathrm{o}},q_{\mathrm{2}})+\mathit{\beta }{\displaystyle \frac{\partial {\mathit{\psi }}_{\mathrm{2}}^{\mathrm{o}}}{\partial x}}-U_{\mathrm{o}}{\displaystyle \frac{\partial }{\partial x}}\left(\mathrm{\Delta }{\mathit{\psi }}_{\mathrm{2}}^{\mathrm{o}}+(k_{\mathrm{d}}^{\mathrm{o}}{)}^{\mathrm{2}}{\mathit{\psi }}_{\mathrm{2}}^{\mathrm{o}}\right)\\ & ={\mathit{\kappa }}_{\mathrm{o}}\mathrm{\Delta }{\mathit{\psi }}_{\mathrm{2}}^{\mathrm{o}}-{\mathit{\nu }}_{\mathrm{o}}{\mathrm{\Delta }}^{\mathrm{4}}{q}_{\mathrm{2}}^{\mathrm{o}},\end{array}\end{array}$$

where the subscript $(\cdot {)}_i^{\mathrm{o}}$, $i=\mathrm{1},\mathrm{2}$, denotes the oceanic layers. The term $q_i^{\mathrm{o}}$ denotes the PV and ${\mathit{\psi }}_i^{\mathrm{o}}$ denotes the streamfunction. The term $k_{\mathrm{d}}^{\mathrm{o}}$ represents the baroclinic deformation wavenumber corresponding to the Rossby radius of deformation L_d. Additionally, $J(A,B)=A_x{B}_y-A_y{B}_x$ denotes the Jacobian operator. A large-scale vertical shear $(U_{\mathrm{o}},-U_{\mathrm{o}})$, with the same strength but opposite directions, is assumed in the background to induce baroclinic instability. In the dissipation terms on the right-hand sides of the equations, besides hyperviscosity ${\mathit{\nu }}_{\mathrm{o}}{\mathrm{\Delta }}^{\mathrm{4}}{q}_i^{\mathrm{o}}$, only Ekman friction ${\mathit{\kappa }}_{\mathrm{o}}\mathrm{\Delta }{\mathit{\psi }}_{\mathrm{2}}^{\mathrm{o}}$ is used, with κ_o indicating the strength of the friction applied to the surface layer of the ocean flow. Additionally, the relationships between the PV and the streamfunction are described by the following equations:

$$\begin{array}{}\text{(16)}& {\displaystyle }& {\displaystyle }{q}_{\mathrm{1}}^{\mathrm{o}}=\mathrm{\Delta }{\mathit{\psi }}_{\mathrm{1}}^{\mathrm{o}}-{\displaystyle \frac{(k_{\mathrm{d}}^{\mathrm{o}}{)}^{\mathrm{2}}}{\mathrm{2}}}({\mathit{\psi }}_{\mathrm{1}}^{\mathrm{o}}-{\mathit{\psi }}_{\mathrm{2}}^{\mathrm{o}}),\text{(17)}& {\displaystyle }& {\displaystyle }{q}_{\mathrm{2}}^{\mathrm{o}}=\mathrm{\Delta }{\mathit{\psi }}_{\mathrm{2}}^{\mathrm{o}}+{\displaystyle \frac{(k_{\mathrm{d}}^{\mathrm{o}}{)}^{\mathrm{2}}}{\mathrm{2}}}({\mathit{\psi }}_{\mathrm{1}}^{\mathrm{o}}-{\mathit{\psi }}_{\mathrm{2}}^{\mathrm{o}}).\end{array}$$

Here, q_(⋅) denotes the potential vorticity at the (⋅) layer, following the notation used in Qi and Majda (2016) and Chen et al. (2021, 2022b). Its precise definition may vary slightly depending on the component (e.g., atmosphere or ocean) and is clarified in the corresponding context.

2.4 Atmosphere dynamics: a two-layer saturated precipitating quasi-geostrophic model

The PQG model, as detailed in Smith and Stechmann (2017), Edwards et al. (2020b), and Hu et al. (2021), captures synoptic-scale dynamics, particularly in extratropical regions. This model extends traditional quasi-geostrophic dynamics by incorporating moisture-related thermodynamic processes, such as water vapor dynamics, cloud formation, phase transitions, and precipitation. This integration provides a more realistic representation of the large-scale meteorological patterns commonly observed in these areas.

In this study, we focus on the structure and statistics of total water within fully saturated domains. By employing a fully saturated or convective setup, we simplify the analysis and simulation, avoiding the complexities introduced by phase changes or convective thresholds that could complicate the coupling model. The two-layer fully saturated PQG model with periodic boundary conditions imposed in both the x and y directions yields the following equations:

$$\begin{array}{}\text{(18)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial q_{\mathrm{1}}^a}{\partial t}}& +J({\mathit{\psi }}_{\mathrm{1}}^a,q_{\mathrm{1}}^a)+\mathit{\beta }{\displaystyle \frac{\partial {\mathit{\psi }}_{\mathrm{1}}^a}{\partial x}}+U_{\mathrm{a}}{\displaystyle \frac{\partial }{\partial x}}\left(\mathrm{\Delta }{\mathit{\psi }}_{\mathrm{1}}^a+(k_{\mathrm{d}}^{\mathrm{a}}{)}^{\mathrm{2}}{\mathit{\psi }}_{\mathrm{1}}^a\right)\\ & =-{\mathit{\nu }}_{\mathrm{a}}{\mathrm{\Delta }}^{\mathrm{4}}{q}_{\mathrm{1}}^a,\end{array}\text{(19)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial q_{\mathrm{2}}^{\mathrm{a}}}{\partial t}}& +J({\mathit{\psi }}_{\mathrm{2}}^{\mathrm{a}},q_{\mathrm{2}}^{\mathrm{a}})+\mathit{\beta }{\displaystyle \frac{\partial {\mathit{\psi }}_{\mathrm{2}}^{\mathrm{a}}}{\partial x}}-U_{\mathrm{a}}{\displaystyle \frac{\partial }{\partial x}}\left(\mathrm{\Delta }{\mathit{\psi }}_{\mathrm{2}}^{\mathrm{a}}+(k_{\mathrm{d}}^{\mathrm{a}}{)}^{\mathrm{2}}{\mathit{\psi }}_{\mathrm{2}}^{\mathrm{a}}\right)\\ & ={\mathit{\kappa }}_{\mathrm{a}}\mathrm{\Delta }{\mathit{\psi }}_{\mathrm{2}}^{\mathrm{a}}-{\mathit{\nu }}_{\mathrm{a}}{\mathrm{\Delta }}^{\mathrm{4}}{q}_{\mathrm{2}}^{\mathrm{a}},\end{array}\text{(20)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial M}{\partial t}}& +J({\mathit{\psi }}_m,M)+v_m{\displaystyle \frac{\partial M_{\mathrm{bg}}}{\partial y}}\\ & =-{\displaystyle \frac{V_{\mathrm{p}}}{\mathrm{\Delta }z}}{q}_{\mathrm{t},\mathrm{m}}+E-{\mathit{\nu }}_{\mathrm{a}}{\mathrm{\Delta }}^{\mathrm{4}}M.\end{array}\end{array}$$

In this model, the subscript $(\cdot {)}_i^{\mathrm{a}}$, $i=\mathrm{1},\mathrm{2}$, denotes the atmospheric layers. The term $q_i^{\mathrm{a}}$ represents the PV, ${\mathit{\psi }}_i^{\mathrm{a}}$ the streamfunction, and M the balanced moisture variable. Also, V_p is the precipitation fall speed and E is the evaporation rate. The relationships between vorticity, streamfunction, and PV follow as outlined in Eqs. (16)–(17). The baroclinic deformation wavenumber $k_{\mathrm{d}}^{\mathrm{a}}$ corresponds to the Rossby radius of deformation L_d. The total water mixing ratio, q_t,m, is related to the equivalent potential temperature, θ_e,m, as follows:

$$\begin{array}{}\text{(21)}& {\displaystyle }{q}_{\mathrm{t},\mathrm{m}}=M-G_{\mathrm{M}}{\mathit{\theta }}_{\mathrm{e},\mathrm{m}},\end{array}$$

where G_M is the ratio of the background vertical gradients of total water mixing ratio, and θ_e,m the equivalent potential temperature. The equation for θ_e,m is given by

$$\begin{array}{}\text{(22)}& {\displaystyle }{\mathit{\theta }}_{\mathrm{e},\mathrm{m}}={\displaystyle \frac{\stackrel{\mathrm{̃}}{L}}{L_{\mathrm{ds}}}}{\displaystyle \frac{{\mathit{\psi }}_{\mathrm{2}}^{\mathrm{a}}-{\mathit{\psi }}_{\mathrm{1}}^{\mathrm{a}}}{\mathrm{\Delta }z}}.\end{array}$$

The background PV, PV_i,bg, and balanced moisture, M_bg, are defined as

$$\begin{array}{}\text{(23)}& {\displaystyle }P{V}_{i,\mathrm{bg}}& {\displaystyle }=(-\mathrm{1}{)}^i{\left({\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }z}}\right)}^{\mathrm{2}}{\left({\displaystyle \frac{\stackrel{\mathrm{̃}}{L}}{L_{\mathrm{ds}}}}\right)}^{\mathrm{2}}(\mathrm{2}{U}_{\mathrm{a}}y),\text{(24)}& {\displaystyle }{M}_{\mathrm{bg}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }z}}{\displaystyle \frac{\stackrel{\mathrm{̃}}{L}}{L_{\mathrm{ds}}}}\left(\mathrm{2}{U}_{\mathrm{a}}y\right),\end{array}$$

where the vertical shear $(U_{\mathrm{a}},-U_{\mathrm{a}})$, with the same strength but opposite directions, is assumed in the background to induce baroclinic instability in the atmosphere. Additionally, the characteristic length scale $\stackrel{\mathrm{̃}}{L}$ and the saturated deformation length scale L_ds satisfy the following relationship:

$$\begin{array}{}\text{(25)}& {\displaystyle }(k_{\mathrm{d}}^{\mathrm{a}}{)}^{\mathrm{2}}=\mathrm{8}{\left({\displaystyle \frac{\stackrel{\mathrm{̃}}{L}}{L_{\mathrm{du}}}}\right)}^{\mathrm{2}}={\displaystyle \frac{\mathrm{8}}{\mathrm{1}+G_{\mathrm{M}}}}{\left({\displaystyle \frac{\stackrel{\mathrm{̃}}{L}}{L_{\mathrm{ds}}}}\right)}^{\mathrm{2}}.\end{array}$$

Remark 1.

The background zonal velocities U^a and U^o in the atmospheric and oceanic QG models are implicitly included in the computation of drag forces and torques on the ice floes. These background zonal velocities affect the flow fields, thereby influencing the ice floe dynamics and trajectories in the coupled system.

2.5 Thermodynamic interactions between atmosphere, ocean, and sea ice components

2.5.1 Floe freezing and melting

In the atmospheric model, the magnitude of precipitation is primarily regulated by the variable q_t, which is derived from the balanced moisture variable M and the streamfunctions ${\mathit{\psi }}_{\mathrm{1}}^{\mathrm{a}}$ and ${\mathit{\psi }}_{\mathrm{2}}^{\mathrm{a}}$. The evolution Eq. (20) governs the transport of moisture as well as the falling of precipitation, represented by $-\frac{V_{\mathrm{p}}}{\mathrm{\Delta }z}{q}_{\mathrm{t},\mathrm{m}}$, and evaporation, denoted by E.

In the initial coupling between ice floes and the atmosphere, several factors are taken into account. First, the precipitation fall speed V_p is assumed to be constant, with the amount of precipitation at each location determined by the variable q_t,m. Second, the evaporation rate E varies between oceanic and ice surfaces (Omstedt et al., 1997; Bintanja and Selten, 2014). Specifically, evaporation is described by

$$\begin{array}{}\text{(26)}& {\displaystyle }E(x,y)=\left\{\begin{array}{ll}{E}_{\mathrm{i}}& \text{if }(x,y)\in \mathbb{Q},\\ E_{\mathrm{o}}& \text{if }(x,y)\in {\mathbb{R}}^{\mathrm{2}}\setminus \mathbb{Q},\end{array}\right.\end{array}$$

where ℚ represents regions covered by ice floes, and ℝ² denotes the double-periodic domain of the ocean. In the ice–atmosphere coupled model, the total water content q_t is significantly influenced by spatial variations in the evaporation rate across ice-covered and open ocean regions. Over the open ocean, evaporation typically occurs at a higher and more consistent rate due to the relatively warm surface temperature and the absence of insulating ice cover. This facilitates a steady exchange of moisture between the ocean surface and the atmosphere. In contrast, the evaporation rate over sea ice is markedly lower, as ice is a barrier that reduces direct contact between the water surface and the atmosphere. This difference in evaporation rates is further complicated by the presence of ice floes, which introduce localized variations in the evaporation rate. Specifically, larger ice floes can significantly suppress evaporation in their vicinity, while smaller floes may have a negligible impact. Consequently, the model is required to account for the heterogeneous distribution of evaporation rates to accurately simulate total water content and its dynamics within the coupled sea ice–atmosphere system. This difference is crucial for understanding the overall moisture balance, cloud formation, and precipitation patterns in polar regions. In particular, the evaporation rate can be parameterized as a function of floe size (radius r_l) and the distance from the center of the floe:

$$\begin{array}{}\text{(27)}& {\displaystyle }E(x,y){\displaystyle }=E_{\mathrm{o}}-\sum _{l=\mathrm{1}}^L{\stackrel{\mathrm{̃}}{E}}_{\text{ice}}(x,y;x_l,y_l,r_l),\end{array}$$

where

$$\begin{array}{}\text{(28)}& {\displaystyle }{\stackrel{\mathrm{̃}}{E}}_{\text{ice}}(x,y;x_l,y_l,r_l){\displaystyle }=(r_l-r_{\mathrm{thr}}{)}^{+}\cdot E_{\text{ice}}(x,y;x_l,y_l,r_l).\end{array}$$

The influence of the lth ice floe on the evaporation rate at location (x,y) is modeled as

$$\begin{array}{}\text{(29)}& \begin{array}{rl}& E_{\text{ice}}(x,y;x_l,y_l,r_l)=a_l\cdot {\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{2}\mathit{\pi }{r}_l^{\mathrm{2}}}}}\\ & \times \exp \left(-{\displaystyle \frac{(x-x_l{)}^{\mathrm{2}}+(y-y_l{)}^{\mathrm{2}}}{\mathrm{2}{r}_l^{\mathrm{2}}}}\right),\end{array}\end{array}$$

where a_l is a scaling parameter given by $a_l=\frac{\mathrm{1}}{L_{\text{domain}}}\frac{r_l-r_{\mathrm{thr}}}{r_{\mathrm{thr}}}.$ Here, r_l is the radius of the lth ice floe, and (x_l,y_l) denotes its location. The standard deviation of the Gaussian distribution, representing the spatial influence of the ice floe, is σ(r_l)=r_l. To account for the negligible impact of small ice floes, we introduce the threshold r_thr in the ramp function $(r_l-r_{\mathrm{thr}}{)}^{+}$, defined as

$$\begin{array}{}\text{(30)}& (r_l-r_{\mathrm{thr}}{)}^{+}=\left\{\begin{array}{ll}{r}_l-r_{\mathrm{thr}}& \text{if }{r}_l>r_{\mathrm{thr}},\\ \mathrm{0}& \text{if }{r}_l\le r_{\mathrm{thr}},\end{array}\right.\end{array}$$

where r_thr=20 km is the chosen threshold radius. Ice floes with a radius smaller than r_thr are assumed to have a negligible impact on the evaporation rate distribution. The term $(r_l-r_{\mathrm{thr}}{)}^{+}$ in Eq. (28) ensures that only the portion of the ice floe's radius exceeding this threshold contributes to the reduction in evaporation. As a result, the overall evaporation rate E(x,y) is determined by subtracting the cumulative effect of all significant ice floes from the constant base evaporation rate over the ocean. The influence of each ice floe is represented by a Gaussian function centered at its location, with its effect scaled by the floe's radius beyond the threshold.

Additionally, precipitation/snow that falls onto the surface of an ice floe contributes to an increase in its height, depth, or mass (Massom et al., 2001; Provost et al., 2017). This increase is uniformly distributed across each floe. The rate of increase in the depth of the floe l over time Δt is given by

$$\begin{array}{}\text{(31)}& {\displaystyle \frac{\mathrm{d}{h}^l}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{w}}\mathit{\pi }{r}_j^{\mathrm{2}}}}\underset{{\mathrm{\Omega }}_j}{\iint }{\displaystyle \frac{V_{\mathrm{p}}}{\mathrm{\Delta }z}}{q}_{\mathrm{t},\mathrm{m}}(x,y)\mathrm{d}x\mathrm{d}y,\end{array}$$

where Ω_l is the subdomain of the lth ice floe and ρ_w is the density of water. It is important to note that with the assumption of ice floes being circular in shape, the mass m^l of each floe can be calculated using the formula

$$\begin{array}{}\text{(32)}& {\displaystyle }{m}^l={\mathit{\rho }}_{\mathrm{ice}}\mathit{\pi }(r^l{)}^{\mathrm{2}}{h}^l,\end{array}$$

where ρ_ice represents the density of ice, and r^l denotes the radius of the floe.

It is noted that the assumption that precipitation increases ice floe thickness is based on the idealization that all precipitation is interpreted as snowfall. This assumption is most valid during colder months and in regions where snow predominates, and may not hold during summer or fall in the MIZ, as highlighted in Boisvert et al. (2023). As a first step toward investigating the influence of precipitation on ice floes, we adopt this assumption to explore its impact within an idealized modeling framework.

2.5.2 Radiation and transfer of radiant energy

The term “insolation” refers to the amount of incoming solar energy (Berger, 1978). To calculate the energy absorbed by sea ice from sunlight, the intercepted energy is multiplied by one minus the albedo value (Miller et al., 2010; Wang et al., 2016). The albedo quantifies the fraction of light reflected away from the ice, so one minus the albedo represents the fraction of light energy absorbed. The total energy absorbed by the lth sea ice over time Δt is given by

$$\begin{array}{}\text{(33)}& {\displaystyle }{E}_l^{\mathrm{ice}}=\mathrm{\Delta }t\underset{{\mathrm{\Omega }}_l}{\iint }\mathit{\gamma }\left(\mathit{x}\right)E_{\mathrm{s}}(\mathrm{1}-\mathit{\alpha })\mathrm{d}x\mathrm{d}y.\end{array}$$

Here, E_s represents the solar insolation, known as the solar constant. While this value may differ in polar winter conditions, it is treated as constant in this study as a simplification appropriate for the short simulation period. α is the albedo of the ice. The variable γ represents the fraction of radiation that penetrates the cloud layer, which varies inversely with the total water content; more total water results in a lower γ value.

To assess how radiation impacts the size of each ice floe, we assume that any change in the ice floe's volume is uniformly distributed across the floe. The formula for the reduction in the depth of the lth floe over time Δt can be expressed as

$$\begin{array}{}\text{(34)}& {\displaystyle \frac{\mathrm{d}{h}^l}{\mathrm{d}t}}=-{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{ice}}\mathit{\pi }(r^l{)}^{\mathrm{2}}}}{\displaystyle \frac{E_l^{\mathrm{ice}}}{C_{\mathrm{ice}}}},\end{array}$$

where C_ice represents the specific heat capacity of the ice.

3 Data assimilation of the coupled system

The coupled model developed above can serve as a testbed for evaluating the accuracy of inferring missing observations of ice floe trajectories and the underlying ocean fields in the presence of clouds through DA. This topic is crucial not only for understanding dynamical coupling and inference capabilities but also as a prerequisite for effective forecasting.

In the presence of cloud cover, DA encounters the challenge of missing observations, raising the question “How can we recover unobserved variables and fields when observations are absent?” However, due to the complexity of the coupled model, using it directly as the forecast model in ensemble DA proves computationally expensive. For instance, the full coupled atmosphere–ice–ocean model is implemented as a sequential MATLAB script without parallelization. A single forward simulation over one model day (i.e., 86 400 s) requires approximately 3.70 CPU hours, corresponding to a wall-clock time of 0.82 h on an Apple M1 Max processor with 32 GB of RAM. It is therefore not feasible to use this full-order model directly in data assimilation settings due to its high computational cost. To address this issue, we develop a low-cost surrogate model for non-linear DA with partial observations. The surrogate model simplifies the DEM system by replacing contact forces with white noise and utilizing reduced-order models for the atmosphere and ocean in spectral space (Chen, 2023; Majda and Harlim, 2012; Majda et al., 2019; Chen et al., 2022b). Such a strategy significantly reduces computational costs while preserving the essential features of the original model.

Despite the lack of adequate polar observations, satellite data on ice floes and wind conditions can assist in recovering unobserved fields (Haugen, 2014). We employ the local ensemble transform Kalman filter (LETKF; Hunt et al., 2007; Bishop et al., 2001) to integrate model forecasts and partial observations through Bayesian updating (Evensen, 2003; Houtekamer and Mitchell, 2005).

Remark 2.

In the presence of sea ice floes, the direct coupling between the atmosphere and ocean is relatively weak compared to other components of the system, such as the atmosphere–ice floe interaction. Moreover, the atmospheric and oceanic models have different timescales – specifically, the ocean evolves on a much slower timescale than the atmosphere – and the current study focuses on short-term simulations. As a simplification, we choose to neglect the atmosphere–ocean coupling, which remains a reasonable assumption within the scope of the present modeling framework.

3.1 Cheap surrogate forecast models for the coupled system

DA involves two key steps: forecasting and analysis. The forecasting step utilizes a forecast model to obtain predicted statistics, which does not necessarily need to be the true underlying dynamics. In practice, the actual system is never fully known, and comprehensive models can be costly to run, particularly for ensemble forecasts. Consequently, using appropriate, inexpensive surrogate or reduced-order models (ROMs) that capture the essential dynamical and statistical features of the original system is often crucial for facilitating practical DA (Majda, 2016; Farrell and Ioannou, 1993; Berner et al., 2017; Branicki et al., 2018; Majda and Chen, 2018; Li and Stechmann, 2020; Harlim and Majda, 2008; Chen and Majda, 2018; Kang and Harlim, 2012).

In the coupled atmosphere–ocean–sea ice model developed above, the most computationally intensive components, such as detailed floe–floe and ocean–atmosphere interactions, are replaced with suitable surrogates that approximate their behavior with significantly lower complexity. The challenging task of modeling contact forces among sea ice floes is simplified by representing these forces as white noise. By alleviating the computational burden, ROMs enable more frequent updates to the model state using new observational data, thereby enhancing the accuracy and reliability of predictions.

While it would be ideal to perform a baseline data assimilation experiment using the full coupled atmosphere–ocean–ice model, such a setup is computationally prohibitive. A single 66 d forecast with the full model requires approximately 244 CPU hours, and an ensemble of 50 members, which is a relatively small ensemble size, would require over 12 000 CPU hours per assimilation cycle. Furthermore, the full model operates on a much finer spatial grid than the reduced-order surrogate, resulting in a substantially higher-dimensional state variable and significantly higher computational costs.

3.1.1 Surrogate forecast model for sea ice motions

The surrogate forecast model utilized in DA starts with a simplified DEM model for the sea ice motions:

$$\begin{array}{}\text{(35)}& {\displaystyle }\mathrm{d}{\mathit{x}}^l={\mathit{v}}^l\mathrm{d}t,\text{(36)}& {\displaystyle }\begin{array}{rl}\mathrm{d}{\mathit{v}}^l& ={\displaystyle \frac{\mathrm{1}}{m^l}}\left(\underset{\text{Ocean drag force}}{\underbrace{{\mathcal{D}}_{\mathrm{o}}^l\left({\mathit{u}}_{\mathrm{o}}-{\mathit{v}}^l\right)}}+\underset{\text{Atmosphere drag force}}{\underbrace{{\mathcal{D}}_{\mathrm{a}}^l\left({\mathit{u}}_{\mathrm{a}}-{\mathit{v}}^l\right)}}\right)\mathrm{d}t\\ & +{\mathit{\sigma }}_{v}d{\mathbf{W}}_v^l\left(t\right),\end{array}\end{array}$$

where x^l is the position of the ice floe, v^l is the velocity of the ice floe contributed by the ocean drag force ${\mathcal{D}}_{\mathrm{o}}^l({\mathit{u}}_{\mathrm{o}}-{\mathit{v}}^l)$ and the atmospheric drag force ${\mathcal{D}}_{\mathrm{a}}^l({\mathit{u}}_{\mathrm{a}}-{\mathit{v}}^l)$, m^l is the mass of the ice floe, σ_v represents the intensity of the stochastic term, and ${\mathbf{W}}_v^l\left(t\right)$ is a Wiener process representing the stochastic component. The stochastic forcing is used to effectively approximate the instantaneous floe–floe interactions in the forecast system, significantly reducing the computational cost while at the same time providing a physically motivated source of model-error variance that maintains adequate ensemble spread. Equation (35) describes the change in position over time, while Eq. (36) describes the change in velocity, incorporating both ocean and atmospheric drag forces, as well as a stochastic term to account for random contact forces.

3.1.2 Surrogate forecast model for the atmosphere dynamics

Recall that the atmospheric wind fields are modeled using the two-layer saturated PQG framework. However, directly applying this two-layer PQG model in DA proves to be computationally expensive. To overcome this challenge, we propose constructing stochastic surrogate models for a limited number of spectral modes of the streamfunction. Given that the upper atmospheric variable is observed, it is essential for the surrogate models to maintain the connection between the upper and lower layers. To this end, we employ both barotropic and baroclinic formulations for the surrogate model that automatically couple the dynamics of the two layers.

The spectral representation of the streamfunctions at two layers is given by

$$\begin{array}{}\text{(37)}& {\mathit{\psi }}_{\mathrm{1}}=\sum _{\mathit{k}}\widehat{{\mathit{\psi }}_{\mathrm{1},\mathit{k}}}\left(t\right)e^{i\mathit{k}\cdot \mathit{x}},{\mathit{\psi }}_{\mathrm{2}}=\sum _{\mathit{k}}{\widehat{\mathit{\psi }}}_{\mathrm{2},\mathit{k}}\left(t\right)e^{i\mathit{k}\cdot \mathit{x}},\end{array}$$

where ψ₁ and ψ₂ are the streamfunctions at the first and second layers, respectively. The summation is over wavenumbers k, $\widehat{{\mathit{\psi }}_{\mathrm{1},\mathit{k}}}\left(t\right)$ and ${\widehat{\mathit{\psi }}}_{\mathrm{2},\mathit{k}}\left(t\right)$ are the time-dependent spectral coefficients for the respective layers, i is the imaginary unit, and x represents the spatial coordinates. On the other hand, the barotropic and baroclinic formulations for each Fourier wavenumber yield

$$\begin{array}{}\text{(38)}& {\widehat{\mathit{\psi }}}_{\mathrm{bt},\mathit{k}}={\displaystyle \frac{\widehat{{\mathit{\psi }}_{\mathrm{1},\mathit{k}}}+{\widehat{\mathit{\psi }}}_{\mathrm{2},\mathit{k}}}{\mathrm{2}}},{\widehat{\mathit{\psi }}}_{\mathrm{bc},\mathit{k}}={\displaystyle \frac{\widehat{{\mathit{\psi }}_{\mathrm{1},\mathit{k}}}-{\widehat{\mathit{\psi }}}_{\mathrm{2},\mathit{k}}}{\mathrm{2}}},\end{array}$$

where ${\widehat{\mathit{\psi }}}_{\mathrm{bt},\mathit{k}}$ and ${\widehat{\mathit{\psi }}}_{\mathrm{bc},\mathit{k}}$ are the barotropic and baroclinic streamfunctions for the wavenumber k, respectively. The subscript “bt” therefore labels the barotropic mode, while “bc” labels the baroclinic component. Linear stochastic surrogate models can be constructed for the barotropic and baroclinic streamfunctions for each Fourier wavenumber (Chen, 2023; Majda and Harlim, 2012):

$$\begin{array}{}\text{(39)}& {\displaystyle }\begin{array}{rl}\mathrm{d}{\widehat{\mathit{\psi }}}_{\mathrm{bt},\mathit{k}}& =(-{\mathit{\gamma }}_{\mathrm{bt},\mathit{k}}+i{\mathit{\omega }}_{\mathrm{bt},\mathit{k}}){\widehat{\mathit{\psi }}}_{\mathrm{bt},\mathit{k}}\mathrm{d}t+f_{\mathrm{bt},\mathit{k}}\mathrm{d}t\\ & +{\mathit{\sigma }}_{\mathrm{bt},\mathit{k}}\mathrm{d}{W}_{\mathrm{bt},\mathit{k}}\left(t\right),\end{array}\text{(40)}& {\displaystyle }\begin{array}{rl}\mathrm{d}{\widehat{\mathit{\psi }}}_{\mathrm{bc},\mathit{k}}& =(-{\mathit{\gamma }}_{\mathrm{bc},\mathit{k}}+i{\mathit{\omega }}_{\mathrm{bc},\mathit{k}}){\widehat{\mathit{\psi }}}_{\mathrm{bc},\mathit{k}}\mathrm{d}t+f_{\mathrm{bc},\mathit{k}}\mathrm{d}t\\ & +{\mathit{\sigma }}_{\mathrm{bc},\mathit{k}}\mathrm{d}{W}_{\mathrm{bc},\mathit{k}}\left(t\right),\end{array}\end{array}$$

where ${\mathit{\gamma }}_{\mathrm{bt},\mathit{k}},{\mathit{\omega }}_{\mathrm{bt},\mathit{k}},f_{\mathrm{bt},\mathit{k}},{\mathit{\sigma }}_{\mathrm{bt},\mathit{k}}$ and ${\mathit{\gamma }}_{\mathrm{bc},\mathit{k}},{\mathit{\omega }}_{\mathrm{bc},\mathit{k}},f_{\mathrm{bc},\mathit{k}},{\mathit{\sigma }}_{\mathrm{bc},\mathit{k}}$ are parameters obtained by calculating the statistical quantities of the time series, while W_bt,k(t) and W_bc,k(t) are independent Wiener processes. In particular, γ_bt,k and ω_bt,k are estimated using the cross-correlation function of the time series of ${\widehat{\mathit{\psi }}}_{\mathrm{bt},\mathit{k}}$:

$$\begin{array}{}\text{(41)}& {\displaystyle }{\mathrm{XC}}_{\mathrm{bt},\mathit{k}}\left(t\right)=\sin \left({\mathit{\omega }}_{\mathrm{bt},\mathit{k}}t\right)e^{-{\mathit{\gamma }}_{\mathrm{bt},\mathit{k}}t},\end{array}$$

where XC_bt,k(t) represents the cross-correlation function for the barotropic component at wavenumber k. The other two parameters can be approximated using the following formulas:

$$\begin{array}{}\text{(42)}& {\displaystyle }& {\displaystyle }{T}_{\mathrm{bt},\mathit{k}}={\displaystyle \frac{{\mathit{\gamma }}_{\mathrm{bt},\mathit{k}}}{{\mathit{\omega }}_{\mathrm{bt},\mathit{k}}^{\mathrm{2}}+{\mathit{\gamma }}_{\mathrm{bt},\mathit{k}}^{\mathrm{2}}}},{\mathit{\theta }}_{\mathrm{bt},k}={\displaystyle \frac{{\mathit{\omega }}_{\mathrm{bt},\mathit{k}}}{{\mathit{\omega }}_{\mathrm{bt},\mathit{k}}^{\mathrm{2}}+{\mathit{\gamma }}_{\mathrm{bt},\mathit{k}}^{\mathrm{2}}}},\text{(43)}& {\displaystyle }& {\displaystyle }{f}_{\mathrm{bt},\mathit{k}}={\displaystyle \frac{m_{\mathrm{bt},\mathit{k}}(T_{\mathrm{bt},\mathit{k}}-i{\mathit{\theta }}_{\mathrm{bt},\mathit{k}})}{T_{\mathrm{bt},\mathit{k}}^{\mathrm{2}}+{\mathit{\theta }}_{\mathrm{bt},\mathit{k}}^{\mathrm{2}}}},\text{(44)}& {\displaystyle }& {\displaystyle }{\mathit{\sigma }}_{\mathrm{bt},\mathit{k}}=\sqrt{{\displaystyle \frac{\mathrm{2}{E}_{\mathrm{bt},\mathit{k}}{T}_{\mathrm{bt},\mathit{k}}}{T_{\mathrm{bt},\mathit{k}}^{\mathrm{2}}+{\mathit{\theta }}_{\mathrm{bt},\mathit{k}}^{\mathrm{2}}}}},\end{array}$$

where m_bt,k is the mean, E_bt,k is the variance, and T_bt,k and θ_bt,k are the real and imaginary parts of the decorrelation time, respectively. Concurrently, the streamfunctions in the two layers from the surrogate model, ψ₁ and ψ₂, can be recovered using the following relations:

$$\begin{array}{}\text{(45)}& {\displaystyle }\widehat{{\mathit{\psi }}_{\mathrm{1},\mathit{k}}}={\widehat{\mathit{\psi }}}_{\mathrm{bt},\mathit{k}}+{\widehat{\mathit{\psi }}}_{\mathrm{bc},\mathit{k}},{\widehat{\mathit{\psi }}}_{\mathrm{2},\mathit{k}}={\widehat{\mathit{\psi }}}_{\mathrm{bt},\mathit{k}}-{\widehat{\mathit{\psi }}}_{\mathrm{bc},\mathit{k}}.\end{array}$$

The velocities in the two layers of the atmosphere are given by

$$\begin{array}{}\text{(46)}& {\displaystyle }{\mathit{u}}_j^{\mathrm{a}}=(u_j^{\mathrm{a}},v_j^{\mathrm{a}})=\left({\displaystyle \frac{\partial {\stackrel{\mathrm{̃}}{\mathit{\psi }}}_j}{\partial y}},-{\displaystyle \frac{\partial {\stackrel{\mathrm{̃}}{\mathit{\psi }}}_j}{\partial x}}\right),\end{array}$$

where $j=\mathrm{1},\mathrm{2}$ denotes the layer. Here, ${\mathit{u}}_j^{\mathrm{a}}$ represents the velocity vector, with $u_j^{\mathrm{a}}$ and $v_j^{\mathrm{a}}$ being the velocity components in the x and y directions, respectively. The surrogate model streamfunction for layer j is given by

$$\begin{array}{}\text{(47)}& {\displaystyle }{\stackrel{\mathrm{̃}}{\mathit{\psi }}}_j=\sum _{\mathit{k}\in {\mathbf{K}}_{\mathrm{r}}}{\widehat{\mathit{\psi }}}_{j,\mathit{k}}\left(t\right)e^{i\mathit{k}\cdot \mathit{x}},\end{array}$$

where k denotes the wavenumber vector, ${\widehat{\mathit{\psi }}}_{j,\mathit{k}}\left(t\right)$ is the time-dependent spectral coefficient for layer j, and K_r represents the set of retained spectral modes.

3.1.3 Surrogate forecast model for the ocean dynamics

Recall that the ocean fields are generated by the two-layer QG model. Only the surface layer of the ocean fields is used in coupling the atmosphere and the ice floes. To simplify the model, a stochastic surrogate model is constructed solely for the near-surface ocean streamfunction, represented as

$$\begin{array}{}\text{(48)}& \begin{array}{rl}\mathrm{d}{\widehat{\mathit{\psi }}}_{\mathrm{o},\mathit{k}}& =(-{\mathit{\gamma }}_{\mathrm{o},\mathit{k}}+i{\mathit{\omega }}_{\mathrm{o},\mathit{k}}){\widehat{\mathit{\psi }}}_{\mathrm{o},\mathit{k}}\mathrm{d}t+f_{\mathrm{o},\mathit{k}}\mathrm{d}t\\ & +{\mathit{\sigma }}_{\mathrm{o},\mathit{k}}\mathrm{d}{W}_{\mathrm{o},\mathit{k}}\left(t\right),\end{array}\end{array}$$

where ${\widehat{\mathit{\psi }}}_{\mathrm{o},\mathit{k}}$ is the spectral coefficient of the ocean streamfunction at wavenumber k, γ_o,k is the damping coefficient, ω_o,k is the frequency, f_o,k is the forcing term, σ_o,k is the noise amplitude, and W_o,k(t) represents the Wiener process.

The parameters γ_o, ω_o, f_o, and σ can be approximated similarly using the formulas provided in Eqs. (41) and (44). The velocity of near-surface ocean current can be obtained similarly in Eqs. (46) and (47).

3.2 Uncertainty in observing ice floes

The presence of clouds poses a significant challenge in observing the location of sea ice floes, especially when using optical and infrared satellite imagery (Reiser et al., 2020; Hyun and Kim, 2017; Wright and Polashenski, 2018). Clouds can obscure the surface, which prevents sensors from capturing clear images of the ice floes underneath. Such an issue causes gaps in observational data, where the exact position and extent of the floes cannot be accurately determined. Note that even when clouds are partially transparent, the scattering and absorption of light can distort observed images, leading to biases in the inferred location and size of the floes. In addition, clouds can form and dissipate relatively rapidly, which complicates the temporal consistency of observations. Consequently, cloud cover introduces a substantial source of uncertainty in observing and monitoring sea ice floes, especially in regions with frequent or persistent cloudiness.

The observability of ice floes varies significantly with their size. Due to their extensive surface area, large ice floes are generally easier to detect in satellite imagery, even in cloudy conditions. These large-sized floes may still be partially visible through breaks in the clouds or cloud-affected imagery, allowing for some degree of position and movement tracking. However, the accuracy of this tracking can be affected since cloud-induced distortions may obscure edges or lead to errors in determining the exact boundaries of floes. In contrast, small ice floes are more susceptible to being completely obscured by clouds. Due to their limited size, these floes can be entirely covered by clouds, resulting in missing data that complicates the assimilation of accurate ice dynamics into models. Identifying small floes is also hampered by the limited spatial resolution of many satellite sensors. Even minor cloud-induced distortions can make small floes undetectable or misclassified, significantly underestimating their presence and distribution. The difficulty in observing small floes gives a challenge for accurately modeling sea ice behavior, as these smaller elements can play a critical role in the overall dynamics and thermodynamics of sea ice, particularly in processes such as melt pond formation and interactions with ocean and atmospheric conditions.

To summarize, while large floes can often be partially observed even in cloudy conditions, the observation of small floes is much more challenging. This leads to large uncertainties in identifying their location, extent, and influence on the sea ice system. This disparity emphasizes the need for developing advanced observational techniques and DA methods to bridge the gaps and reduce the inaccuracies caused by cloud cover.

To represent observational uncertainty in DA, we use the total water content q_t(x,t) as a controlling factor. Above each floe, we calculate the mean total water content [q_t(x,t)] at time t:

$$\begin{array}{}\text{(49)}& {\displaystyle }\left[q_{\mathrm{t}}\right(\mathit{x},t\left)\right]={\displaystyle \frac{\mathrm{1}}{\left|{\mathrm{\Omega }}_l\right|}}\underset{{\mathrm{\Omega }}_l}{\int }{q}_{\mathrm{t}}(\mathit{x},t)\mathrm{d}\mathit{x}.\end{array}$$

This mean value encapsulates the spatial distribution of water content above each ice floe, serving as an approximation for the uncertainty in observations. Variations in [q_t(x,t)] from one floe to another can indicate the degree of uncertainty inherent in the observational data, as it reflects the heterogeneity in the physical characteristics of the ice. In particular, we set a threshold, ${\stackrel{\mathrm{̃}}{q}}_{\mathrm{t}}$, such that the observational uncertainty ${\mathit{\sigma }}_l^{\mathrm{obs}}$ is given, for the lth floe, by

$$\begin{array}{}\text{(50)}& \begin{array}{rl}& {\mathit{\sigma }}_l^{\mathrm{obs}}\left({\mathit{x}}_l\right(t\left)\right)=\\ & \left\{\begin{array}{ll}\mathrm{5}\times {\mathrm{10}}^{\mathrm{2}}\text{m}& \text{if }\left[q_{\mathrm{t}}\right(\mathit{x},t\left)\right]<{\stackrel{\mathrm{̃}}{q}}_t\left(\text{small observation uncertainty}\right),\\ \mathrm{2}{r}_l& \text{if }\left[q_{\mathrm{t}}\right(\mathit{x},t\left)\right]\ge {\stackrel{\mathrm{̃}}{q}}_t\left(\text{large observation uncertainty}\right),\end{array}\right.\end{array}\end{array}$$

where r_l denotes the radius of the lth floe and x_l its trajectory. It is important to note that when the mean total water content over the floe is high, which indicates significant cloud cover and can be classified as unobserved, its position can still be approximated. However, these estimates are often highly inaccurate. Consequently, in the data-assimilation setting, we assign floes classified as unobserved a markedly inflated observational uncertainty, taken here as twice their radius.

In summary, the total water content q_t interacts with ice floes in two different ways. First, the spatial coverage of an ice floe over the ocean affects the local evaporation, which in turn modifies the distribution of q_t over ice and ocean. Second, the mean total water content [q_t] above each floe is used to parameterize cloud-related observation uncertainty in floe localization, with higher [q_t] indicating thicker clouds and thus higher uncertainty.

Remark 3.

An ice floe is classified as “unobserved” whenever the mean total water content over the floe exceeds the threshold, $\left[q_{\mathrm{t}}\right]\ge {\stackrel{\mathrm{̃}}{q}}_{\mathrm{t}}$ – cf. Eq. (50). In that case we inflate the observation-error standard deviation to ${\mathit{\sigma }}_l^{\mathrm{obs}}=\mathrm{2}{r}_l$, with r_l the floe radius. Because this value is of the same order as the floe's diameter, the associated observation exerts negligible influence on the analysis. Among all the test cases, the smallest floe radius is $r_l^{min}=\mathrm{8}\times {\mathrm{10}}^{\mathrm{3}}\text{m}$, so an unobserved floe takes at least $\mathrm{2}{r}_l=\mathrm{1.6}\times {\mathrm{10}}^{\mathrm{4}}\text{m}$ of uncertainty, clearly distinguished from the observed case with ${\mathit{\sigma }}_l^{\mathrm{obs}}=\mathrm{5}\times {\mathrm{10}}^{\mathrm{2}}\text{m}$.

3.3 Local ensemble transform Kalman filter

The LETKF (Hunt et al., 2007; Bishop et al., 2001) functions similarly to standard ensemble Kalman filters but implements filtering within localized domains. Specifically, the LETKF improves the process by partitioning the global data set into smaller, overlapping regions. Each region is updated independently using local observations. It effectively reduces long-range error correlations and enhances computational efficiency. This localization is essential for large-scale applications. It makes LETKF more scalable and accurate when dealing with spatially extended systems.

LETKF applies filtering exclusively in physical space for both the trajectories of ice floes, which are in the Lagrangian frame of reference, and the streamfunctions of the atmosphere and ocean that are in the Eulerian frame of reference. This unified approach maintains consistency across different types of geophysical data, potentially simplifying the assimilation process and directly addressing spatial correlations.

4 Numerical simulation results of the coupled atmosphere–ocean–sea ice model

4.1 Setup

The domain size considered here is 400 km×400 km in the marginal ice zone. The size distribution of ice floes follows a power law, represented as $p\left(r\right)=a{k}^a/r^{a+\mathrm{1}}$ (Stern et al., 2018), where the constants k and a are parameters of the model. The numerical integration time step, Δt, is set to 58.2 s. Three different distributions of floe sizes are considered in this study and each distribution comprises 48 floes, designed to simulate different ice coverage and collision frequency:

• regime I: large-sized ice floes nearly cover the entire ocean domain, resulting in frequent collisions;

• regime II: medium-sized ice floes cover approximately half of the ocean domain, leading to less frequent collisions; and

• regime III: small-sized ice floes occupy only a small portion of the ocean domain, where collisions are rare.

These different scenarios are designed to explore the dynamics of ice floes under different coverage conditions and assess their impacts on the atmospheric model. The detailed parameters of the DEM model are listed in Table 2.

Table 2 Parameters in the ice floe model.

The coupled atmosphere–ocean-ice model employs different numerical schemes for each component. The atmospheric and oceanic equations are discretized in space using a spectral method and integrated in time with an adaptive third-order Runge–Kutta scheme, following Qi and Majda (2016) and Edwards et al. (2020a, b). The sea ice component is simulated using the DEM, with time integration using a forward Euler scheme to resolve floe dynamics dominated by contact and drag forces.

The model parameters for the two-layer QG model are summarized in Table 3. These parameters are specifically tailored to represent conditions typical of the high-latitude ocean regime (Qi and Majda, 2016).

Table 3 Parameters in QG equations for high-latitude ocean regime.

The parameters relevant to the atmospheric regime within the framework of the saturated PQG equations are detailed in Table 4 (Edwards et al., 2020a, b; Smith and Stechmann, 2017; Hu et al., 2021). The atmospheric and oceanic systems exhibit multiscale characteristics both temporally and spatially. Specifically, typical atmospheric wind velocities range from 8 to 10 m s⁻¹ (equivalent to 800 km d⁻¹). They are significantly faster than ocean current speeds, which are around 0.1 m s⁻¹ (corresponding to 10 km d⁻¹). In addition, atmospheric winds exhibit rapid temporal changes compared to the relatively slower temporal variations observed in ocean currents. In addition, in this work, the terms “total water” and “total water mixing ratio” are used interchangeably, following the convention adopted in Edwards et al. (2020a, b), to maintain consistency and clarity. It is worth noting that the dimensional time steps for the atmospheric and ocean models are obtained by rescaling their respective non-dimensional counterparts using their different characteristic scales appropriate to each component under the Arctic regime; for example, the atmospheric model's time step of 58.2 s corresponds to a non-dimensional value of $\mathrm{\Delta }t=\mathrm{5}\times {\mathrm{10}}^{-\mathrm{4}}$, based on the characteristic timescale T employed in the PQG model.

Table 4 Parameters in saturated PQG equations for the high-latitude atmosphere regime.

4.2 Simulated atmospheric and ocean fields

Figure 2 shows the time evolution of the PVs (PV₁ and PV₂) in the atmosphere from t=404.4 to 808.8 h and 1213.3 h. These snapshots illustrate the development of the PV structures and their time evolution. They indicate the dynamic interactions within the coupled atmosphere–ocean–sea ice system.

Figure 2 PV₁ and PV₂ for atmosphere at different time instances.

Figure 3 shows the upper-layer PV field of the ocean. It reveals distinct patterns in the evolution of the PV structures with both large- and small-scale localized features. At time t=404.4 h, the PV field looks relatively uniform, with moderate gradients indicating steady circulation. When the time arrives at t=808.8 h, the PV distribution becomes more complicated, with sharper gradients and more eddies. The spatial patterns suggest an increased vorticity intensity, with the appearance of turbulent mixing and instabilities in certain regions. At t=1213.3 h, the PV field is further strengthened. It is more concentrated, indicating a highly dynamic state with stronger interactions among different oceanic layers. The gradient patterns and isolated vortices highlight the non-linear nature of the response in the ocean to external forcing and the internal energy transfers occurring within the coupled atmosphere–ice–ocean system.

Figure 3 PV₂ field for the upper-layer ocean at different time instances.

4.3 Simulated floe trajectories

The trajectories presented in Fig. 4 help understand the role of floe mass in modulating the sea ice motions. The floes have decreasing mass from top to bottom. The results here illustrate the relationship between floe size, mass, and external forces, such as wind stress and ocean currents.

The largest floe is shown in the top panel. The floe has a higher inertia and exhibits a smoother and more streamlined trajectory. The results in this panel suggest that, under similar forcing conditions, larger and heavier floes resist rapid changes in motion and respond more gradually to external perturbations due to their momentum. These floes dominate ice dynamics over longer timescales. They illustrate slower and more predictable moving patterns.

In contrast, the medium-sized floe in the middle panel shows a more irregular motion, with noticeable deviations from a linear path. This is because its reduced mass allows for more immediate responses to short-term variations in external forcing, such as fluctuations in atmospheric wind or ocean currents. A balance between inertia and responsiveness appears in the motion of this floe. It highlights that mid-sized floes can exhibit complex dynamics even under relatively steady forcing conditions.

The bottom panel shows the trajectory of a small-sized floe. Due to the low mass and inertia of the floe, it has the strongest irregular motion. Therefore, the motion of small-sized floes is highly susceptible to external forces. Such floes contribute to the more chaotic aspects of sea ice movement, especially in regions with strong atmospheric or oceanic variability.

Figure 4 Trajectories of the largest floe in three different regimes (I, II, and III, from top to bottom), all starting from the same initial location, marked by a blue rectangle.

4.4 Simulated precipitation in the sea ice–atmosphere coupled model

From a modeling perspective, accurately capturing the differences in evaporation rates between ice-covered and open ocean areas is essential for simulating the coupled atmosphere–sea ice system. The model must account for the spatial variability introduced by ice floes, including their size, distribution, and movement, as these factors significantly influence local evaporation rates.

We calculate the spatially averaged total water content for each individual ice floe as

$$\begin{array}{}\text{(51)}& {\displaystyle }\left[q_{\mathrm{t}}{]}_l\right(t)={\displaystyle \frac{\mathrm{1}}{\left|{\mathrm{\Omega }}_l\right|}}\underset{{\mathrm{\Omega }}_l}{\int }{q}_{\mathrm{t}}(x,y,t)\mathrm{d}\mathit{x}.\end{array}$$

Here, [q_t]_l(t) represents the average total water content over the area of the lth ice floe at time t. The integral is taken over the domain Ω_l, which denotes the spatial extent of the lth floe, with $\left|{\mathrm{\Omega }}_l\right|$ representing the area of this ice floe.

Figure 5 Floe distribution at different times for three different regimes (top to bottom: regime I, regime II, and regime III), with different quantities: first column: floe distribution at t=808.8 h; second and third columns: total water at t=808.8 and 1213.3 h; fourth column: averaged total water at t=808.8 h.

By integrating over the floe's area, we obtain a representative measure of the total water content, which can be used to analyze the floe's contribution to the overall hydrological balance in the sea ice–atmosphere system.

The first three columns of Fig. 5 display the ice floes and total water content across different regimes. In the last column of Fig. 5, the time-averaged spatially averaged total water content, $\stackrel{\mathrm{‾}}{\left[q_{\mathrm{t}}{]}_l\right(t)}$, is plotted as a function of different floe radii. The results demonstrate that $\stackrel{\mathrm{‾}}{\left[q_{\mathrm{t}}{]}_l\right(t)}$ decreases as the floe radius increases. This trend suggests that larger ice floes tend to have lower average water content due to reduced moisture exchange.

5 Numerical DA results

5.1 Setup

The floe locations and upper atmospheric data are the only observational inputs used in the DA process for the coupled system. The forecast time step, Δt, is set to 58.2 s, allowing the model to resolve finer temporal details. Observations are available every 1500 numerical integration time steps, which corresponds to approximately Δt^obs=24.2 h. This frequency aligns with the acquisition of satellite images, which occurs roughly every 24 h. The observational uncertainty for the ice floe locations, when there is no cloud cover, is set to 0.5 km. In contrast, the noise in the streamfunction is quantified as 20 % of the standard deviation of the streamfunction at a given grid point, denoted as std(ψ_x). The total duration of the DA spans 1601.5 h, or approximately 66.73 d. While the truth is computed from a high-resolution model with 128×128 grids, the reduced-order forecast model in spectral space contains 16×16 modes to enhance computational efficiency. The ensemble size in the LETKF is 300. The localization radius is set to 200 km, which limits the influence of observations to a localized region, thereby reducing spurious correlations in the DA. The parameters used in the DA are summarized in Table 5. It is worth noting that the atmospheric model uses a time step of Δt=58.2 s, while the DEM sea ice model is updated every 80 atmospheric steps, giving Δt≈1.29 h. Observations for data assimilation are assumed to be available every 1500 atmospheric steps, or Δt^obs≈24.2 h.

Table 5 Parameters in the DA.

In the numerical tests, we consider two different settings for the observations of ice floes: (1) plentiful observations and (2) sparse observations. In case (1), most of the floes (≥ 70 %) in the domain are observed during the observation period, whereas in case (2), only a few floes (about 30 %) are observed during the same period.

5.2 Recovered atmosphere and ocean fields

5.2.1 Recovered atmosphere fields

Figure 6 presents the performance of DA in recovering the streamfunctions of the atmosphere flow fields, which are crucial in describing the large-scale circulation of the atmosphere, at two different layers: the upper-layer streamfunction, ψ₂ (top panel), and the near-surface layer streamfunction, ψ₁ (bottom panel). The comparison is at t=970.61 h.

Despite the noisy observations, the posterior mean of the upper-layer streamfunction, ψ₂, closely matches the truth. The recovered near-surface streamfunction, ψ₁, shows slight deviations from the truth because ψ₁ is not directly observed. Its inference relies on a combination of observed floe trajectories and ψ₂, which introduces some uncertainty into the recovery of ψ₁. Nevertheless, despite the slight increase in uncertainty, the spatiotemporal patterns and amplitudes are largely recovered in the DA solution, demonstrating the effectiveness of DA in estimating unobserved states.

Figure 6 Recovered atmosphere streamfunction at the upper layer, ψ₂ (top panel), and near-surface layer, ψ₁ (bottom panel) at t=970.61 h: comparison among the truth, the posterior and prior means, and upper-layer observations.

5.2.2 Recovered upper-layer ocean field

Figure 7 shows the performance of DA in recovering the upper-layer ocean streamfunction. The figure compares the true ocean streamfunction (top panel) with the posterior mean (bottom panel) at several time instances.

The DA effectively restores the magnitude of the ocean fields and captures some of the large-scale patterns. The absence of small-scale features in the DA results is expected, as floe movements are primarily driven by atmospheric forces, with ocean drag playing a relatively minor role. Consequently, the ocean field has limited observability, resulting in less accurate state estimation.

Figure 7 Recovered upper-layer ocean streamfunction: comparison among the truth (a, b, c) and the posterior mean at different time instances, (d, e, f) and (g, h, i).

Figure 8 DA results for the ice floe trajectories. Panels (a)–(c) depict the floes' trajectories over time with the posterior mean (green), observations (red), and the true trajectory (blue).

5.3 Recovered floe trajectories

Figure 8 compares the true and estimated trajectories of three ice floes under sparse observational conditions in regime II over a period of 66.7 d. Similar results were observed in the other two regimes, so they are omitted here. In each panel, the posterior mean trajectory is shown in green, the observations in red, and the true trajectory in blue. The observational uncertainty, introduced via Eq. (50), shows that when cloud cover, ${\stackrel{\mathrm{̃}}{q}}_t$, remains below a certain threshold, the uncertainty stays small (around 500 m), allowing for accurate floe location recovery. However, when cloud cover exceeds this threshold, the observational uncertainty increases significantly, which can lead to substantial inaccuracies in the data assimilation of floe trajectories.

Figure 8a shows the trajectory of an ice floe with a radius of 22.73 km, which is relatively large for this regime (see Fig. 5). The results demonstrate a reasonably close match between the posterior mean and the true trajectory, despite sparse observational data. This is expected, as the value of $\stackrel{\mathrm{‾}}{\left[q_{\mathrm{t}}{]}_l\right(t)}$ is smaller for larger floes (as shown in the averaged total water content in the second row of Fig. 5), leading to lower observational uncertainty for floes with larger radii.

Panel (b) shows the trajectory of an ice floe with a radius of 22.03 km, slightly smaller than the one in panel (a), which may lead to greater observational uncertainty at certain times. During these instances, the red observational points deviate significantly from the true trajectory. Nonetheless, the DA still successfully recovers the floe trajectory, despite some highly inaccurate observations.

Finally, panel (c) illustrates the trajectory of an ice floe with a radius of 12.23 km, significantly smaller than those in panels (a) and (b). In this case, ${\stackrel{\mathrm{̃}}{q}}_t$ remains relatively high throughout, indicating larger observational uncertainties. Consequently, the red observation points deviate significantly from the true locations at most observational time instances. Although the DA recovered trajectory is less accurate compared to the previous two cases, it still captures the general trajectory pattern relative to the true path.

Table 6 Normalized RMSE of the posterior (third column) and prior (fourth column) DA results in regime I.

Table 7 Normalized RMSE of the posterior (third column) and prior (fourth column) DA results in regime II.

Table 8 Normalized RMSE of the posterior (third column) and prior (fourth column) DA results in regime III.

5.4 Skill scores

Finally, we present the normalized root-mean-square error (RMSE) of the posterior mean estimate as a skill score for DA across different regimes. The state variables in this DA study consist of (1) the locations of the ice floes, and (2) the streamfunctions of the upper and near-surface atmosphere, as well as the ocean. Since the ice floe positions are Lagrangian quantities, while the streamfunctions are represented on Eulerian grids, the normalized RMSE skill scores are defined differently for the floe trajectories x and the streamfunction ψ:

$$\begin{array}{}\text{(52)}& \begin{array}{rl}& {\mathrm{RMSE}}_{\mathit{x}}={\displaystyle \frac{\sqrt{\frac{\mathrm{1}}{M}{\sum }_{k=\mathrm{1}}^M\left|{\mathit{x}}^{\text{truth}}\right(t_k)-{\mathit{x}}^{\text{posterior}}(t_k){|}^{\mathrm{2}}}}{\sqrt{\frac{\mathrm{1}}{M}{\sum }_{k=\mathrm{1}}^M\left|{\mathit{x}}^{\text{truth}}\right(t_k\left)\right(t_k){|}^{\mathrm{2}}}}},\\ & {\mathrm{RMSE}}_{\mathit{\psi }}=\\ & {\displaystyle \frac{\sqrt{\frac{\mathrm{1}}{NM}{\sum }_{j=\mathrm{1}}^N{\sum }_{k=\mathrm{1}}^M\left|{\mathit{\psi }}^{\text{truth}}\right({\mathit{x}}_j,t_k)-{\mathit{\psi }}^{\text{posterior}}({\mathit{x}}_j,t_k){|}^{\mathrm{2}}}}{\sqrt{\frac{\mathrm{1}}{NM}{\sum }_{j=\mathrm{1}}^N{\sum }_{k=\mathrm{1}}^M\left|{\mathit{\psi }}^{\text{truth}}\right(t_k\left)\right({\mathit{x}}_j,t_k){|}^{\mathrm{2}}}}},\end{array}\end{array}$$

where N is the number of grid points and M is the number of observations, (⋅)^truth(t_k) is the true position of the floes at time t_k, and (⋅)^posterior(t_k) represents the model's posterior mean estimate of the floe position at time t_k. When the normalized RMSE exceeds 1, it indicates a loss of skill in the estimation, as the error in the posterior mean reaches the level of the equilibrium standard deviation.

The DA skill scores across the three regimes are shown in Tables 6–8. For the fully observed atmospheric state variable, ${\mathit{\psi }}_{\mathrm{2}}^{\mathrm{a}}$, and the unobserved variable, ${\mathit{\psi }}_{\mathrm{1}}^{\mathrm{a}}$, the DA demonstrates strong performance. This is expected, as ${\mathit{\psi }}_{\mathrm{1}}^{\mathrm{a}}$ and ${\mathit{\psi }}_{\mathrm{2}}^{\mathrm{a}}$ are closely linked by barotropic and baroclinic dynamics, allowing for relatively accurate recovery of ${\mathit{\psi }}_{\mathrm{1}}^{\mathrm{a}}$ even though it is not directly observed. In contrast, the estimation of the ocean streamfunction, ψ^o, is less accurate. This is because in the coupled model, the floe movements are primarily driven by atmospheric forces, while ocean drag has a relatively weak influence. As a result, the ocean field has limited observability, leading to less accurate state estimation. On the other hand, the estimation of floe trajectories is generally accurate, with accuracy slightly decreasing from high-concentration (regime I) to low-concentration regimes (regime III). This is due to the lower-concentration regime experiencing more cloud cover (Fig. 5), which amplifies observational uncertainties.

6 Conclusion

In this paper, we developed an idealized coupled atmosphere–ocean–ice model to investigate the effects of clouds on sea ice dynamics in the MIZ. Our model integrates the DEM to simulate the movement and interaction of individual ice floes, combined with a two-layer QG ocean model and a two-layer PQG atmospheric model that includes saturated precipitation processes. This framework facilitates studying the interactions between atmospheric, ocean, and sea ice floes. It specifically focuses on cloud-induced radiative and precipitation effects on sea ice evolution. The paper addresses both forward (model simulation) and inverse (DA) problems. For the former, we study the interactions between different model components; for the latter, we focus on recovering unobserved floe trajectories obscured by cloud cover and inferring ocean and atmospheric fields using limited observations.

The results from this study show the significant influence of clouds on sea ice dynamics. They also help in understanding the non-trivial interactions within the coupled atmosphere–ocean–ice system. Integrating idealized modeling with advanced DA techniques forms a powerful tool for enhancing the prediction of Arctic climate processes, particularly in the MIZ, where small-scale interactions are critical.

There are several topics for future research work. One is about further refining the model to incorporate more sophisticated thermodynamic processes, and the other is to test its performance under varying climate scenarios. Moreover, the continued development of DA schemes that can handle more comprehensive observational data sets (Curry et al., 1996) will be essential for advancing our understanding of sea ice dynamics and for improving the accuracy of Earth system models. In addition, it would also be worthwhile to explore how the statistical properties of ice floe trajectories vary with floe size. The current model does not include the Coriolis force, as it focuses on short-term, localized ice floe dynamics where collision and drag forces are dominant. As a future direction, we plan to incorporate the Coriolis effect to evaluate its impact on larger-scale and longer-term ice floe motion. As a future direction, we aim to extend the current framework by incorporating two-way coupling between the ocean and ice components, as well as a fully coupled atmosphere–ocean interaction, to better capture feedback mechanisms and improve the accuracy of ocean state estimation in data assimilation settings. Because existing in situ and remote sensors sample only limited portions of the stratified upper ocean, its full three-dimensional state remains only partially observable in the present experiments; a quantitative assessment of this limitation is left for future study.

Appendix A: Contact forces in the coupling model

A1 Normal contact force

In the mathematical framework described in Eqs. (36) and (5), the normal contact force between the lth and jth ice floes is governed by Hooke's law of linear elasticity. Specifically, the force ${\mathit{f}}_{\mathit{n}}^{lj}$ is expressed as

$$\begin{array}{}\text{(A1)}& {\displaystyle }{\mathit{f}}_{\mathit{n}}^{lj}=c^{lj}{E}^{lj}\mathit{\delta }{n}^{lj}{\mathit{n}}^{lj},\end{array}$$

where n^lj denotes the normal vector originating from the center of floe l and pointing towards the center of floe j. The function ${\mathit{\delta }}_n^{lj}$, representing a delta function, ensures that the contact force between the jth and lth floes is non-zero only when they are in contact:

$$\begin{array}{}\text{(A2)}& {\displaystyle }{\mathit{\delta }}_n^{lj}=\left\{\begin{array}{ll}\mathrm{1}& \text{if }\left|d^{lj}\right|-(r^l+r^j)<\mathrm{0},\\ \mathrm{0}& \text{otherwise. }\end{array}\right.\end{array}$$

Here, r^l and r^j are the radii of the respective floes, and d^lj is the distance between their centers. The chord length c^lj, oriented transversely across the cross-sectional area, is calculated as

$$\begin{array}{}\text{(A3)}& {\displaystyle }{c}^{lj}={\displaystyle \frac{\mathrm{1}}{d^{lj}}}\sqrt{\mathrm{4}\left(d^{lj}{)}^{\mathrm{2}}\right(r^{max}{)}^{\mathrm{2}}-{\left((d^{lj}{)}^{\mathrm{2}}-(r^{min}{)}^{\mathrm{2}}+(r^{max}{)}^{\mathrm{2}}\right)}^{\mathrm{2}}},\end{array}$$

where E^lj represents Young's modulus of elasticity. The parameters c^lj and d^lj, crucial for understanding the mechanical interactions, are illustrated in Fig. A1.

A2 Tangent contact force

The tangential force at the contact interface between two ice floes is proportional to their relative tangential velocities. The direction of this force, denoted by t^lj, is perpendicular to the normal direction. This force, ${\mathit{f}}_{\mathit{t}}^{lj}$, is defined by

$$\begin{array}{}\text{(A4)}& {\displaystyle }{\mathit{f}}_{\mathit{t}}^{lj}=c^{lj}{G}^{lj}\mathrm{\Delta }{v}_t^{lj}{\mathit{t}}^{lj},\end{array}$$

where c^lj refers to the chord length previously defined in Eq. (A3), G^lj is the shear modulus, and $\mathrm{\Delta }{v}_t^{lj}$ is the difference in velocity along the tangential direction between the lth and jth ice floes.

Figure A1 Illustration of geometrical quantities for computing the tangential and normal contact force.