Satellite-Based Data Assimilation System for the Initialization of Arctic Sea Ice Concentration and Thickness Using CICE5

Lee, Jeong-Gil; Ham, Yoo-Geun

doi:10.3389/fclim.2022.797733

ORIGINAL RESEARCH article

Front. Clim., 01 March 2022

Sec. Predictions and Projections

Volume 4 - 2022 | https://doi.org/10.3389/fclim.2022.797733

This article is part of the Research Topic The State of the Art in Climate Predictions View all 9 articles

Satellite-Based Data Assimilation System for the Initialization of Arctic Sea Ice Concentration and Thickness Using CICE5

$\nJeong-Gil Lee$ Jeong-Gil Lee

Yoo-Geun Ham^*

Department of Oceanography, Chonnam National University, Gwangju, South Korea

The satellite-derived sea ice concentration (SIC) and thickness (SIT) observation over the Arctic region are assimilated by implementing the Ensemble Optimal Interpolation (EnOI) into the Community Ice CodE version 5.1.2 (CICE5) model. The assimilated observations are derived from Special Sensor Microwave Imager/Sounder (SSMIS) for the SIC, European Space Agency's (ESA) Soil Moisture and Ocean Salinity mission (SMOS) for the SIT of the thin ice, and ESA's CryoSat-2 satellite for the SIT of the thick ice. The SIC, and SIT observations are assimilated during 2000–2019, and 2011–2019, respectively. The quality of the reanalysis is evaluated by comparing with observation and modeled data. Three data assimilation experiments are conducted: noDA without data assimilation, Ver1 with SIC assimilation, and Ver2 with SIC and SIT assimilation. The climatological bias of the SIC in noDA was reduced in Ver1 by 29% in marginal ice zones during boreal winter, and 82% in pan-Arctic ocean during boreal summer. The quality of simulating the interannual variation of sea ice extent (SIE) is improved in all months particularly during boreal summer. The root-mean-square errors (RMSEs) of SIE anomaly in August are significantly reduced compared to noDA. However, the interannual variations of SIT is unrealistic in Ver1 which requires the additional assimilation of the SIT observation. The climatological bias of SIT over the Arctic was further reduced in Ver2 by 28% during boreal winter compared to that in Ver1. The interannual variability of SIT anomalies is realistically simulated in Ver2 by reducing the RMSEs of SIT anomalies by 33% in February, and 28% in August by compared to that in Ver1. The dominant interannual variation extracted by empirical orthogonal function (EOF) of SIT anomalies in Ver2 is better simulated than Ver1. The additional assimilation of SIT improves not only SIT, but also SIC. The climatological bias of SIE and the errors in leading EOF of SIC anomalies in Ver2 is further reduced compared to those in Ver1 during boreal winter. However, improvements led by assimilating additional SIT observation is not clear during boreal summer, possible due to the lack of available SIT observation during this season.

Introduction

Arctic sea ice is an important component of the Earth system in a warming climate that changes most dramatically during recent decades (Serreze et al., 2009; Screen and Simmonds, 2010; Perovich et al., 2017; Kwok, 2018). Arctic sea ice extent and thickness have consistently decreased rapidly in all seasons due to the anthropogenic forcing during last few decades (Kwok and Rothrock, 2009; Stroeve et al., 2012a). In addition to the contribution of anthropogenic forcing, internal variability dominating the Arctic summer circulation trend is also important cause of declining in September Arctic sea ice since 1979 (Ding et al., 2017). The rate of Arctic sea ice decline is unprecedented over at least the past 1,500 years according to reconstructed sea ice historical data (Kinnard et al., 2011), and an ice-free state is expected in the middle of twenty-first century according to climate model predictions under the high emission of the greenhouse gases (GHGs) scenarios (Wang and Overland, 2012; Notz and Stroeve, 2016; Guarino et al., 2020; Notz et al., 2020).

What is more, dramatic decrease of Arctic sea ice has increased interannual variability in the sea ice extent (SIE), particularly from summer to autumn. The interannual variability of the summer Arctic sea ice is non-linearly rising with a warmer climate, so that uncertainties for Arctic sea ice prediction also become large for recent decades (Goosse et al., 2009). Holland et al. (2008) showed the significant increase in standard deviation of September SIE anomaly in the global climate model projections, which eventually would contribute to decrease the predictability for Arctic sea ice in the future. On the other hand, Massonnet et al. (2018) suggested that Arctic sea ice volume will exhibit less year-to-year variability through the enhanced action of growth and melt processes as the ice thins in the future projections.

It has been recently recognized that Arctic sea ice plays an active role on the weather and climate system in the middle latitude in Northern Hemisphere as well as those in Arctic region (Zhang et al., 2012; Kim et al., 2014; Mori et al., 2014; Kug et al., 2015). The recent decrease in the sea ice affects the surface energy budget by the reducing the surface albedo and disappearing an insulating cap of the air-sea interfaces, which can increase the oceanic temperature to amplify the variation of the local near-surface temperature in Arctic (Stroeve et al., 2012b), then it induces the stationary Rossby wave-train to influence the climate in mid-latitudes particularly over the Eurasian and North America (Bader et al., 2011; Cohen et al., 2012; Mori et al., 2014; Kug et al., 2015), north Pacific such as the Aleutian low, and tropical central Pacific (Overland et al., 1999; Kennel and Yulaeva, 2020). That is, Arctic sea ice states are also considered as the cause of the weather and climate variability beyond the polar region (Deser et al., 2015). The possibility of mid-latitude effects of Arctic sea ice leads to increasing demand for the accurate predictions of the Arctic sea ice in recent decades, even though its influence on the mid-latitude climate significantly varies on the period analyzed, implying the necessity of the robust observational evidence (Overland et al., 2016; Blackport and Screen, 2020a,b; Cohen et al., 2020).

In addition, an accurate prediction of the Arctic sea ice leads to economic benefits by predicting the shipping route through the Arctic (i.e., Northern Sea Route, NSR) for commercial shipments between Europe/North America and Asia (Liu and Kronbak, 2010; Smith and Stephenson, 2013). Feasibility of the NSR is only possible when the climate model with high-resolution grid system is developed and the accurate predictions for the sea ice state (concentration and thickness) are guaranteed (Aksenov et al., 2017). Therefore, there has been a rising demand for Arctic sea ice predictions at sub-seasonal to interannual time scales (Eicken, 2013). In addition to the deterministic sea ice forecasts (Blanchard-Wrigglesworth et al., 2015; Wayand et al., 2019), as the forecast of Arctic sea ice is intrinsically uncertain due to the uncertainties originated from the chaotic nature of the sea ice component and irreducible model errors as with other climate system, probabilistic forecasts of Arctic sea ice have been explored recently to quantify the uncertainty in a forecast and to produce the optimal decision making (Dirkson et al., 2019; Gao et al., 2021).

The quality of Arctic sea ice prediction using the numerical model strongly depends on the degree of implementation of sea ice-related dynamical/physical processes in the dynamical model and the quality of the initial condition (Chevallier et al., 2013; Sigmond et al., 2013; Wang et al., 2013; Tietsche et al., 2014). In particular, Tietsche et al. (2014) insisted that shortage of the knowledge of Arctic sea ice initial conditions should impede the skillful predictions. To guarantee the skillful sea ice predictions, it is essential to improve the initial conditions by implementing data assimilation techniques using available observations (Lisæter et al., 2003).

Fortunately, the satellite-derived observational data of the sea ice concentration (SIC) have been successfully produced for over 30 years (Cavalieri et al., 1996). In addition, although retrieval methods are an early stage of development, pan-Arctic sea ice thickness (SIT) satellite observations have been produced since 2002 (Hendricks et al., 2018), and dataset with high spatial gridded resolution of 12.5 km and temporal resolution of 1 day have been also yielded from the 2010s (Tian-Kunze et al., 2014). Thus, several studies have carried out data assimilation of satellite-derived SIC and SIT into dynamical sea ice models using advanced assimilation techniques. Lisæter et al. (2003) used an ensemble Kalman filter (EnKF) to assimilate SIC data derived from Special Sensor Microwave/Imager (SSM/I) carried on board the satellites of the Defense Meteorological Satellite Program (DMSP). Lindsay and Zhang (2006) conducted a nudging scheme to assimilate observed monthly averaged SIC and sea ice velocity independently into the coupled ice-ocean model. Stark et al. (2008) assimilated the SSM/I SIC and sea ice drift data using an optimal interpolation scheme into the Met Office Forecasting Ocean Assimilation Model (FOAM). Caya et al. (2010) implemented the 3D-VAR to assimilate SIC into the Community Ice-Ocean Model (CIOM). The positive impact of assimilation on improving the simulation quality of the modeled SIC was well-established in aforementioned studies.

The next generation of researches emphasized the updates of the SIT in addition to the SIC for the accurate estimates of the growth/decaying speed of the sea ice, and the resultant sea ice openness regulating the ocean-atmosphere interactions are strongly dependent on the SIT. As a pioneer work, Tietsche et al. (2013) and Massonnet et al. (2013, 2015) showed that the quality of SIT in their reanalysis is improved by updating both SIC and SIT with the given SIC increment (i.e., multivariate updates). However, Massonnet et al. (2015) discussed that those results are needed to be further verified because of the inadequate observations of SIT and have large uncertainties in the estimation of SIT. In addition, their methodology has a clear limitation to improve the SIT simulation as the SIT updates is solely dependent on the SIC observation, and the degree of the update is based on the modeled relationship between SIC and SIT.

There are attempts directly assimilating SIT using satellites-based SIT observation (Yang et al., 2014; Chen et al., 2017; Blockley and Andrew Peterson, 2018; Mu et al., 2018; Xie et al., 2018). For example, Yang et al. (2014) shows that the data assimilation of SIT observation of SMOS into an ocean-sea ice model using local ensemble singular evolutive interpolated Kalman filter (SEIK) improves the forecast skill of SIT. Xie et al. (2018) assimilated SIT observation of merged measurements from CS2 and SMOS (CS2SMOS) into the TOPAZ4 system with Deterministic EnKF during boreal winter and showed the improvement in the SIC and sea ice drift in the central Arctic Ocean and Kara Sea.

However, the period of the experiments for the SIT assimilation is only for few months or years in previous researches, therefore, it was rather short to assess the degree of improvement in simulating long-term climatology and the interannual variability (Massonnet et al., 2015). This is, it has not been tried to reconstruct the Arctic sea ice states more than a decade by assimilating both satellite-based SIC and SIT observations.

In this study, we reconstructed the Arctic sea ice states for a decade (i.e., during 2010–2019) by assimilating satellite-derived SIC and SIT into the Los-Alamos sea ice model version 5.1.2 (CICE5) using Ensemble Optimal Interpolation (EnOI) scheme (Backeberg et al., 2014; Mignac et al., 2015; Kaurkin et al., 2016). In addition, the reanalysis quality by assimilating SIC observation is evaluated through the data assimilation experiment for two decades (i.e., during 2000–2019). This article is organized as follows. The description of the model and methods is given in Section Data and Methods. In Section Results, the impact of the data assimilation on the SIC and SIT states are assessed. Summary and Discussions are given in Section Summary and Discussions.

Data and Methods

The Sea Ice Model and Atmospheric/Oceanic Forcing

We used the Community Ice CodE (CICE) version 5.1.2 sea ice model, will be referred as CICE5 for simplicity, developed by the Los Alamos National Laboratory (LANL) (Hunke et al., 2015). The CICE model is a state-of-the-art sea ice model, which is originally developed as a stand-alone version, and it is now set up as a sea ice component of the National Center for Atmospheric Research (NCAR) Community Earth System Model (CESM) (Gent et al., 2011; Hurrell et al., 2013; Danabasoglu et al., 2020), and the others (Maclachlan et al., 2015). The CICE5 is a multicategory sea ice model to describe the ice thickness distribution with five thickness categories (fixed thickness boundaries are 0, 0.64, 1.39, 2.47, and 4.57 m by default). The sea ice in each category consists of seven vertical layers and one snow layer to resolve the sea ice temperature and salinity variations. The CICE5 has a thermodynamic process of mushy layer thermodynamics and a dynamic process of elastic-anistropic-plastic dynamics (Wilchinsky and Feltham, 2006; Turner et al., 2013). The most complex physics parameterization schemes such as melt pond and ridging is recently added in this version of the CICE. The horizontal resolution of this model is 1° approximately with 320 × 384 dimensions on a displaced pole grid and the time step is 1 h.

In this study, CICE5 is forced at the surface with atmospheric and oceanic fields. The National Center for Environmental Prediction Reanalysis 2 (NCEP R2) dataset is used for atmospheric forcing (Kanamitsu et al., 2002). The atmospheric forcing for CICE5 was passed on to the sea ice in flux form at the atmosphere-sea ice boundary (Lee et al., 2019). The near-surface atmospheric variables [6-hourly specific humidity (kg kg⁻¹), 2-m air temperature (K), 10-m U-wind and V-wind (m s⁻¹)] was used to calculate the turbulent fluxes using bulk formulae (Large and Yeager, 2009), and the near-surface air density (kg m⁻³) was calculated using the relationship delineated in Large and Yeager (2009). Monthly mean downward shortwave/longwave radiation (W m⁻²) and total precipitation rate (kg m⁻² s⁻¹) were used as atmospheric forcing.

The Optimal Interpolation Sea Surface Temperature (OISST) version 2 from National Oceanic and Atmospheric Administration (NOAA) (Reynolds et al., 2002) was used for daily sea surface temperature (SST) () forcing. The SST of mixed-layer slab ocean model in CICE5 is nudged toward the OISST data with 21-day Newtonian relaxation time-scale. The sea surface salinity was fixed constant value of 34 psu in this work.

Data for Assimilation and Validation

The satellite-observed daily SIC are obtained from the National Snow and Ice Data Center (NSIDC, http://nsidc.org/data/NSIDC-0079). They are derived from the Scanning Multichannel Microwave Radiometer on the Nimbus-7 satellite and from the Special Sensor Microwave Imager/Sounder (SSMIS) aboard Defense Meteorological Satellite (DMSP) F17 using the bootstrap algorithm (BT) version 3 (Comiso, 2017). The SSMIS BT data is on the polar stereographic projection with a spatial resolution of 25 km. The observational error value of 0.15 is uniformly given for accounting for measurement and representation errors (Janji et al., 2018), according to suggestion of Tonboe and Nielsen (2010). The OISST version 2 reanalysis data for daily SIC is also used for constraining the long-term control simulations to obtain the stationary ensemble perturbations as explained in Section Data Assimilation Scheme (Reynolds et al., 2002).

For SIT, the satellite-observation data have been derived from two different dataset of CryoSat-2 (CS2) and Soil Moisture Ocean Salinity (SMOS). The CS2 SIT is used here (http://nsidc.org/data/RDEFT4). The CS2, the European Space Agency (ESA) satellite mission, Level-4 SIT data based on radar altimeter measurements of sea ice freeboard is used for the area where the SIC is >0.7 (i.e., 70% of total area in the corresponding grid point) (Kurtz et al., 2014; Kurtz and Harbeck, 2017). The period for CS2 assimilation spans from January 1, 2011 to December 31, 2019. The SMOS is obtained from (https://icdc.cen.uni-hamburg.de/thredds/catalog/ftpthredds/smos_sea_ice_thickness/catalog.html). The L3C SMOS, the ESA mission, version 3.1 from January 1, 2011 to April 15, 2018 and version 3.2 from October 1, 2018 to December 31, 2019. Daily thickness data derived from L3B brightness temperatures using a single layer emissivity model is used for assimilating the modeled sea ice thinner than 1.0 m (Kaleschke et al., 2012, 2016; Tian-Kunze et al., 2014). The CS2 data is on the polar-stereographic grid system with a horizontal resolution of 25 km, and the SMOS data is on the NSIDC polar-stereographic grid system with a horizontal resolution of 12.5 km. The uncertainties of the observed SIT data are 1.5 m (0.3 m for random and + 1.2 m for systematic) for CS2 and 0.7 m (averaged of thin thickness) for SMOS, according to Chen et al. (2017). The satellite-based SIT observations are only available in the cold season from September (October) to May (April) for CS2 (SMOS and CS2SMOS), as it is difficult to discriminate between measurements from open ocean and sea ice due to melt ponds forming on the sea ice surface during the melt season (Tilling et al., 2018).

The verification for the quality of updated sea ice field by assimilating satellite data is conducted by utilizing observation and reanalysis data. The daily SIC derived from passive microwave sensor data (SSM/I) using the NASA Team algorithm (NT) is used to validate the reanalysis output (Cavalieri et al., 1996). The SSMIS NT data is used here (http://nsidc.org/data/NSIDC-0051). The Pan-Arctic Ice-Ocean Modeling and Assimilation System (PIOMAS) SIC and SIT are also used (Zhang and Rothrock, 2003). The weekly averaged SIT product version 2.3 of CS2SMOS derived from merging CS2 and SMOS using an optimal interpolation scheme is used to validate the updated thickness field (Ricker et al., 2017). The CS2SMOS data is used here (http://spaces.awi.de/display/CS2SMOS). Since CS2SMOS is created by merging CS2 and SMOS, there is some overlap with the assimilated data set. The grids of CS2SMOS are projected onto the 25 km EASE2 Grid. The model output, observation and reanalysis data are all re-gridded onto the 1 rectilinear latitude/longitude grid system from a displaced pole grid system used to solve the polar singularity problem for model and from Polar-Stereo graphic for satellite observations.

Data Assimilation Scheme

The satellite-derived sea ice data is assimilated using an Ensemble Optimal Interpolation (EnOI) scheme (Figure 1). The EnOI method approximates the background error covariance matrix (B matrix) by using a stationary ensemble member instead of a multiple realization of the model for the Ensemble Kalman Filter (EnKF; Evensen, 2003). Note that, because the mean state of Arctic sea ice is changing rapidly, the EnOI method with stationary B matrix used in this study has a limitation that it would cause a performance degradation compared to techniques considering flow-dependent features such as EnKF. While the EnKF is advantageous by using state-dependent ensemble perturbations, the use of EnKF is often limited as it requires multiple numbers of realization. Therefore, the EnOI is a cost-effective alternative to the EnKF to assimilating the complex dynamic models requiring high computation costs to operate (Oke et al., 2007). As a result, the data assimilation algorithms based on the stationary B matrix is still widely used to produce many SIC reanalysis products such as HadISST with OI (Rayner et al., 2003), NOAA OIv2 with OI (Reynolds et al., 2002), and ECMWF ORAS5 with 3D-Var (Zuo et al., 2019).

FIGURE 1

Figure 1. The schematic of the generation system for forecast error perturbations of $ε_{d + 1}^{n (y)}$ for calculating stationary background error covariance matrix (i.e., B matrix). The black (dotted) line represents the model long-time simulation nudged sea ice and ocean reanalysis with atmospheric forcing and the black (solid) line represents the 1-day forecast experiments forced atmospheric forcing data without restoring sea ice fields. The blue dot represents the initial conditions from the model long-time simulation. The red bidirectional arrow represents the 1-day forecast errors and the green left bracket represents the forecast perturbations using the bred vectors.

In our work, to obtain the stationary ensemble perturbations, two sets of model integrations are performed. First, long-term control simulations are obtained by nudging the daily OISST sea surface temperature (SST) data along with the daily OISIC SIC, and monthly PIOMAS SIT data (dotted black line in Figure 1). The restoring time scale are 21-day for SST, 6-h for SIC, and 5-day for SIT. Note that the small change in the nudging time-scale does not lead systematic changes in the nudged states. The state vector in a control simulation is denoted as $x_{d}^{s (y)}$ , where d is d^th day of a year, and y denotes a year between 1982 and 2019.

Secondly, the perturbed experiment is performed to obtain the ensemble perturbations $ε_{d + 1}^{n (y)}$ , which is the perturbation for the stationary B matrix. Initially (i.e., at January 1st, 1982), the initial condition by adding a random perturbation to a control simulation (i.e., $x_{d}^{s (y)} + {δ x}_{d}^{s (y)}$ ) are generated. Then, 1-day forecasted state $x_{d + 1}^{f (y)}$ is obtained by integration the model for 1-day [i.e., $x_{d + 1}^{f (y)} = M_{d \to d + 1} (x_{d}^{s (y)} + {δ x}_{d}^{s (y)})$ , where M_d→d+1 is the dynamical forecast model from day d to d+1]. Then, the difference from the control simulation at a corresponding day (i.e., $ε_{d + 1}^{n (y)} = x_{d + 1}^{f (y)} - x_{d + 1}^{s (y)}$ ) is defined as the perturbation (red bidirectional arrow in Figure 1). This difference is rescaled by 10% of its original amplitude, then, added to the initial condition at next day (i.e., ${δ x}_{d + 1}^{s (y)} = \frac{ε_{d + 1}^{n (y)}}{10}$ ) (Green left bracket in Figure 1). This procedure is repeated for every day from 1982 to 2019 to obtain the $ε_{d}^{n (y)}$ . As a result, the procedure to obtain the perturbations is similar to those for bred vectors (Toth and Kalnay, 1993, 1997).

The reason to utilize the bred vector approach is to prescribe fast-growing perturbations that can successfully mimics the forecast errors (Toth and Kalnay, 1997). The fast-growing errors dominate total forecast errors by definition, while the contribution of the non- or slow-growing errors for total forecast error becomes smaller in time. Therefore, the bred vector, which is known as a fastest growing perturbation, effectively represents forecast uncertainty with limited the number of ensemble members. Since the number of ensembles (38 in this study) is limited compared to the degrees of freedom of the model (about 10⁵ of CICE5), the ensemble perturbations generated by the procedure mimicking the breeding method grows faster than others to better represent total forecast errors with the limited number of samples.

The procedure for generating B(d) matrix (d denotes a number of day varying from 1 for 1st January to 365 for 31th December) using the perturbations is as follows. B(d) matrix has size of N_b × N_b (N_b is a dimension of the state vector) and $ε_{d}^{n (y)}$ denote the forecast error state vectors at day d in year y with size of N_b × 1. By assigning each year's perturbation sequentially at each column, ε_d for day d is generated whose size is N_b × N_y, where N_y is the number of years (i.e., 38). The stationary B(d) matrix is constructed as $B (d) \approx \frac{1}{N_{y}} [ε_{d}] {[ε_{d}]}^{T}$ . As B matrix is separately formulated by the calendar day, the B matrix accounts to the seasonal spatial differences in the perturbations. The B matrix for the each variable is formulated separately in a univariate sense: the observed information of one variable had no effect on the others during the assimilation procedures.

It should be noted that the generating algorithm of ensemble error perturbation used in this study does not take into account the model uncertainty originated from the physical parameters or the boundary conditions. That is, only the uncertainty led by the instability in the system was considered in the background error covariance matrix B, which can underestimate the amplitude of the error in the background states.

Then, in order to prevent non-physical update of analysis by the observation at unrealistically remote region, the localization scheme of Local Analysis (LA) is carried out (Sakov and Bertino, 2011). A localization radius of 800 km is applied as in Kimmritz et al. (2018) and Massonnet et al. (2015). With LA, the updates at each grid point are applied using observations within a specified distance. In addition to this, LA uses a quasi-Gaussian weight function as a localization function (Gaspari and Cohn, 1999), so observation and background error covariance in observation space are downweighted with increasing distance from the updated grid point like a Localized Ensemble Transform Kalman Filter (LETKF) (Hunt et al., 2007). Post-processing to conserve physical ranges for such variables is applied to keep the analysis states as physically reasonable values. That is, in the case with analysis value of the SIC with negative, and value >1, they are changed to 0, and 1, respectively (Kimmritz et al., 2018). The minimum value of the SIT analysis is set to 0.

In this study, the CICE5 has five categories of the SIT. A fractional sea ice area of each category is a prognostic variable, however, the total SIC, which is denoted by the sum of fractional ice areas for each category, is a diagnostic variable. As the satellite-derived SIC data only provide the aggregated concentration, the SIC increment is calculated using the difference between the modeled total SIC and the observed SIC. Then, analysis increment of total SIC are distributed to each category according to the modeled percentage of the SIC in each category with respect to the total value as method used in Chen et al. (2017). Each category of the SIT is also updated in the same way. Although the EnOI also can facilitate the multivariate scheme as in Zhang et al. (2018) to update each SIT category, it has not been validated whether the static covariance matrix between the aggregated SIT and SIT in each category successfully operates for the multi decades. For this reason, we adopted a relatively simple scheme as aforementioned.

When the SIC increment is positive in the case of both SIC and SIT is 0, the SIC in each category is updated by multiplying the pre-defined ratio for each category. The pre-defined ratios are 60, 20, 10, 5, and 5%, for the category 1–5, respectively. The SIT in each category is updated as a averaged value of the upper and lower thickness limit of each category, which is 0.32, 1.02, 1.93, 3.52, and 6.95 m for the category 1–5, respectively. To update the positive SIT increment in the case of both SIC and SIT is 0, the SIT in each category is updated by multiplying the aforementioned pre-defined ratios, and the SIC in each category is updated by dividing the SIT increment to the aforementioned fixed SIT values.

Experimental Design

The impact of data assimilation of the SIC and SIT is investigated by comparing the quality of the reanalysis in three experiments. “noDA” is a control experiment forced by atmospheric and oceanic forcing fields without any data assimilation. It should be noted that the noDA experiment is not same as long-term control simulation used to obtain the ensemble perturbations in Section Data Assimilation Scheme. The long-term control simulation used to generate ensemble perturbations is the experiment in which sea ice fields are nudged in addition to atmospheric and oceanic forcing, but noDA is driven by forcing only without nudging of sea ice. “Ver1” is an experiment assimilated only SIC, and “Ver2” assimilates both SIC and SIT satellite data. In Ver2, the state vectors for SIC and SIT are constructed and updated separately. The period for the data assimilation in Ver1 is 2000–2019. The data assimilation of the SIT in Ver2 was performed from 2011 to 2019 when SIT observation data is available. Therefore, the reanalysis value before 2011 for Ver2 is identical to Ver1. All data assimilation experiments have a spin-up time of 1-year.

Results

Properties of the Background Error Covariance Matrix

Figure 2 shows the standard deviation (STD) of the SIC perturbation for calculating B matrix at 15th January, April, July, and October. In January, the STD of SIC perturbation is relatively high in the marginal sea ice regions, while it is quite small over the rest of areas including the North Pole (Figure 2A). This denotes that the instability leading the growth of the SIC is relatively high over the marginal sea ice regions including the Atlantic sector of the Barents and Greenland Seas and in the Pacific sector of the Bering Sea and Sea of Okhotsk, and these regions are characterized by strong atmospheric and oceanic variability during boreal winter (Semenov et al., 2015). In April, the overall STD pattern is similar to January, but the amplitude is relatively weak (Figure 2B). The SIC variability in July and October shows a reduced regional difference compared to boreal winter and spring season (Figures 2C,D). These features are evident by the year-to-year variability shown in satellite observations for the same period (not shown). This implies that the background error perturbations would reflect the inter-annual variability of the sea ice in observation and the state-of-the-art numerical models.

FIGURE 2

Figure 2. The year-to-year standard deviation at 15th January (A), April (B), July (C), and October (D) of forecast error perturbations for sea ice concentration (a fraction from 0 to 1) generated by using long-time integration during 1982–2019.

To examine the detailed spatial distribution of the B matrix, Figures 3A,D shows the global covariance patterns of the background error perturbations with respect to two selected points at 15th February, and October, respectively. Note that the localized covariance pattern is shown as in a final form of B matrix for the data assimilation. The covariance pattern with respect to the East Greenland (75°N, 10°W) is tilted in northeast-southwest direction (Figure 3A). This tilted spatial distribution is similarly shown in the SIC and sea surface temperature (SST) regression patterns using monthly anomalies from long-term simulation (i.e., noDA) (Figure 3B). More importantly, the tilted direction of the covariance pattern is consistent with the direction of the climatological sea ice drift, indicating that the East Greenland Current (EGC) is a major contributor to transport sea ice export southward along the continental margin (Figure 3C). This indicates that the sea ice particle mainly moves southwestward, therefore, the covariance with respect to any given point around this region would have relatively high values in southwest direction.

FIGURE 3

Figure 3. The covariance patterns (A,D) in forecast error perturbations of sea ice concentration (a fraction from 0 to 1) generated by using long-time integration, the regression patterns [shaded for sea ice concentration and contour for sea surface temperature (°C); B,E], and mean sea ice drift (cm/s; C,F) in model long-time integration outputs during 1982–2019. The covariance and regression patterns between a selected grid point (A–C: 75N, 10W in Greenland Sea; D–F: 75N, 150W in Beaufort Sea) and global grid points are left and middle column. The upper (lower) panels are results at 15th February (October). All the plots are localized using quasi-Gaussian weight function with localization radius of 800 km.

Similarly, over Beaufort Sea (i.e., 75°N, 150°W), the covariance pattern exhibits zonally-elongated pattern (Figure 3D), and it is also shown in the regression using monthly SST and SIC from long-term simulation (Figure 3E). The Beaufort Gyre (BG), which flows westward, can be responsible for this zonally-elongated pattern by transporting sea ice particle along the BG. This clearly demonstrates that the perturbations for the B matrix reflect the physical properties and balances associated with the sea ice variability.

Figure 4 shows that the time evolution of SIC anomalies over 2000–2019 at the two points used in Figure 3 for the observation and background state of Ver1 and Ver2. The amplitude of the STD of the ensemble perturbations, used to formulate B matrix, are added to, and subtracted from the background state (i.e., before update by assimilation). The observations are generally within the uncertainty range of the model states estimated using stationary ensemble perturbations. It can be also seen that it is also well-represented when zooming in on 2007 (Figures 4B,D). This indicates that the amplitude of B matrix properly measures the amplitude of the model uncertainty to some extent.

FIGURE 4

Figure 4. The temporal evolution over 2000–2019 (A,C) and in 2007 (B,D) of sea ice concentration anomalies (a fraction from 0 to 1) for observation (black) and background states of Ver1 (blue) and Ver2 (red) at a grid point used in Figure 3 (A,B: 75°N, 10°W in East Greenland Current area; C,D: 75°N, 150°W in Beaufort Gyre area). The orange (skyblue) shading marks the uncertainty (standard deviation) of the ensemble perturbations.

The Effect of Data Assimilation

Climatology and Seasonal Cycle

Figure 5 shows that the climatological seasonal cycle of sea ice extent (SIE, the summed area of all Arctic grid points with SIC > 0.15) over 2000–2019 (Figure 5A), and 2011–2019 (Figure 5B). The STD of the analysis errors representing the uncertainty embedded in the final reanalysis product (X_a) is obtained by calculating the diagonal component of the equation $P_{a}^{- 1} = P_{b}^{- 1} + H^{T} R^{- 1} H$ , where P^a, P^b, and R are the covariance matrices of the analysis, background, and observation, respectively, and H is observational operator.

FIGURE 5

Figure 5. The mean seasonal cycle of Arctic sea ice extent index (10⁶km²) during 2000–2019 (A) and during 2011–2019 (B) and temporal evolution over 2000–2019 (C) for SSMIS NT satellite data (black), noDA experiment (gray), Ver1 (blue), and Ver2 (red). The orange shading marks the uncertainty (standard deviation) of the analysis.

Compared to noDA, the seasonal evolution of the SIE in Ver1 and Ver2 becomes quite similar to the satellite observation (i.e., SSMIS NT). For example, in Ver1, the SIE is almost identical to the observation, and exhibited a relatively small difference only from November to February. On the other hand, in noDA, the SIE tends to be systematically underestimated, and overestimated during melting, and freezing season, respectively. The negative SIE bias in noDA during the boreal summer is due to the excessive bottom melting process during previous seasons (not shown). Then, the ice thinning during boreal summer play a role on the increase in the sea ice production during the ice growing season, which may contribute to the excessive SIE during boreal winter (Massonnet et al., 2018). The reduction in the climatological SIE bias in the Ver1 and Ver2 compared to the noDA is particularly robust in August and September, which is a melting season, rather than during January and February, which is the ice freezing season. This improvement in simulating SIE in Ver1 and Ver2 is also clearly shown in monthly time-series (Figure 5C).

It is interesting to note that the SIE in Ver2 is better than that in Ver1, indicating that the assimilation of the SIT observation in Ver2 further reduces the climatological bias in SIC in Ver1 (Fritzner et al., 2019). The SIE in Ver2 exhibits smaller climatological bias during January-February than that in Ver1. The improvement is minimal during 2000–2019, however, it becomes clear during 2011–2019 (Figures 5B,C). On the other hand, SIE in Ver2 is almost identical to that in Ver1 during boreal summer. This might be due to the fact that the SIC is already well-fitted to the observed by assimilating the SIC during boreal summer.

Figure 6 shows the spatial distribution of the climatological SIC of SSMIS NT and bias in noDA and Ver1 in February and August when the improvement in the SIE simulation was robust. The observation shows that the sea ice is distributed not only in polar regions but also in the subpolar regions such as Sea of Okhotsk, Bering Sea, Labrador Sea, and Baffin Bay in February (Figure 6A). The subpolar regions have no sea ice in August (Figure 6D). The climatological bias of SIC in noDA is greater than the observation about 0.3 in most of marginal ice zones, while it is relatively weak over polar region (Figure 6B). The positive bias of SIC in the noDA in most of marginal areas is significantly reduced more than a quarter in Ver1 (Figure 6C). In August, the negative bias of SIC in noDA reached up to −0.4 throughout the central Arctic Ocean with relatively weaker bias amplitude over Canadian Archipelago (CAP) region (Figure 6E). This negative bias is almost disappeared in Ver1 (Figure 6F). This confirms that the SIC data assimilation improves climatological state of the assimilated variable in both freezing and melting seasons.

FIGURE 6

Figure 6. The climatology of sea ice concentration (a fraction from 0 to 1) of SSMIS NT satellite data (A,D) and climatological bias in noDA (B,E) and in Ver1 (C,F) with respect to the satellite observation in February (A–C) and August (D–F) during 2000–2019.

To examine the simulation quality of the climatological SIT, Figures 7A,B show the Northern Hemisphere-averaged monthly SIT climatology, and the spatial distribution of the climatological SIT of CS2SMOS and bias in Ver1 and Ver2 in February. Note that only the calendar months when satellite SIT observation is available are drawn. Over the thin SIT area (SIT of SMOS < 1 m), the SMOS SIT is kept as a nearly constant value (i.e., 0.4 m) from November to March (Figure 7A), possibly due to the cancellation between the SIT thickening of the pre-existing sea ice and the SIT reduction by the new formation of the thin sea ice during freezing season. This constant SIT climatology from November to March is not successfully simulated in noDA, implying that the new sea ice formation is weakly simulated in the model. On the other hand, the SIT climatology is kept as the constant in Ver1. However, the positive SIT bias about 0.2 m is shown throughout all seasons in Ver1. The SIT bias in Ver2 is systematically reduced compared to that in Ver1, while keeping the constant SIT climatology from November to March. This indicates that the simulation quality of the SIT climatology is improved by the direct injection of the SIT observation.

FIGURE 7

Figure 7. The mean seasonal cycle of sea ice thickness (m) averaged over the valid SMOS area (A) and over the valid CryoSat-2 area (B) during 2011–2019 for satellite observation used for assimilation (black), noDA experiment (gray), Ver1 (blue), and Ver2 (red). The orange shading in (A,B) marks the uncertainty (standard deviation) of the analysis. The climatological bias in noDA (C), Ver1 (D), and Ver2 (E) with respect to CS2SMOS in February. The contour lines in (C–E) represent the climatology in each model experiment.

Over the thick SIT area (SIT of CS2 > 1 m), the monthly CS2 climatology is increasing from October to April as it is the freezing season (Figure 7B). The new sea ice formation in those regions is not robust as it is already filled with the sea ice from the beginning of the freezing season. Even though the increasing trend from October to April is shown in noDA, the negative SIT bias is overwhelmed for all seasons. In Ver1, the negative SIT bias is significantly decreased, however, the thickening of SIT during the freezing season is not simulated. For example, the climatological SIT in CS2 is increased from about 1.7 m in October to 2.5 m in April, while that in Ver1 is continuously about 2 m from October to April. On the other hand, in Ver2, the negative SIT bias is reduced by keeping the increasing trend during the freezing season to some extent.

The spatial distribution of the climatological SIT bias in February also demonstrates the positive impact of the SIT data assimilation in Ver2 (Figures 7C–E). The strong negative bias of SIT in noDA is distributed in central Arctic ocean (Figure 7C). The climatological SIT in CS2SMOS exhibited a maximum value over the CAP region (not shown), and this spatial feature is well-simulated in both experiments to some extent (contours in Figures 7D,E). As a result, the SIT bias in both Ver1 and Ver2 is not robust over the CAP region (shadings in Figures 7D,E). In more detail, the amplitude of the SIT bias in Ver2 is <0.5 m over the CAP region, while that is reached up to 1 m in Ver1, supporting the notion that the SIT bias is systematically reduced through the assimilation of the observed SIT data.

The next question is whether the realistic SIT simulation by assimilating SIT observation contribute to improve the simulation quality of the SIC. As it is already pointed out in Figure 5 that the SIE climatology in Ver2 during winter season exhibited smaller bias amplitude compared to that in Ver1, we focused on the change in the SIC and SIT bias by assimilating SIT observation during boreal winter season. Figures 8A,B shows the latitudinally-averaged (65–80°N) absolute value of climatological bias of the SIT, and SIC in Ver1 and Ver2 during DJF season, respectively. Note that our main points are still rigorous with the changes in the latitudinal bands for the averages. The absolute SIT bias in Ver2 is reduced compared to that in Ver1 in most of the Arctic regions; the decrease in the SIT bias is shown around 30–80°E, 120–180°E, 160–90°W, and 30–0°W (Figure 8A). Interestingly, the reduction in the SIC bias by assimilating SIT observation is also evident over 30–80°E, with the relatively stronger changes over the Kara Sea (i.e., 50–80°E) (Figure 8B), even though the reduction of the SIT bias in Ver2 is modest in this region. The large reduction of the SIC bias in this region might be due to the fact that the SIC bias in Ver1 is greatest in those regions, therefore, there would be much room for the improvement in simulating SIC climatology.

FIGURE 8

Figure 8. The meridionally averaged (65–80°N) absolute mean bias (A,B) and mean bias (C,D) between data assimilation experiments (Ver1 represented by blue lines; Ver2 represented by red lines) and satellite observation (CS2SMOS for thickness and SSMIS NT for concentration) for sea ice thickness (m; A,C) and concentration (a fraction from 0 to 1; B,D) during DJF season over 2011–2019. The gray lines in left column represent the difference of absolute mean bias between Ver2 and Ver1.

Over 30–80°E, where the Kara-Barents Sea is located, the sign of the SIT bias in Ver1 is negative, which is corrected to some extent in Ver2 (Figure 8C). The climatological SIC exhibited a strong positive bias in Ver1, and this is significantly reduced in Ver2 (Figure 8D). This indicates that the additional assimilation of the SIT observation leads the increase in the SIT climatology over Kara-Barents Sea, and this contributes to decrease the SIC climatology over this region.

To confirm that the reduction in the bias of the SIC and the SIT in Kara-Barents Sea during boreal winter is led by the increase in the climatological SIT, we conducted an idealized experiment; the individual categories of SIT in initial conditions at 1st February from 2000 to 2019 are increased by 0.05 m, compared to noDA experiment in domain of Kara-Barents Sea (30–80°E, 70–80°N). Note that all other variables are unchanged to examine the sole impact of the SIT change. All-year-averaged monthly (i.e., February) difference in the SIC and SIT between sensitivity experiment and noDA showed that the SIC changes led by the increased SIT exhibited a negative value in this domain (not shown). This confirms that the reduction in the positive SIC bias over Kara-Barents Sea in Ver2 compared to that in Ver1 (i.e., decrease in the SIC climatology in Ver2) is caused by the reduction the negative SIT bias (i.e., increase in the SIT climatology in Ver2) to some extent.

The mechanism representing this inverse relationship between SIC and SIT climatology would be related to thermodynamic or mechanical processes. For instances, the SIT changes can lead the changes in ice conductive flux from ice bottom surface to top surface (Maykut and Untersteiner, 1971). The conductive heat flux of sea ice is inversely proportional to the SIT changes, that is, the increasing in SIT could lead to shrinkage of the conductive flux. Once the conductive heat flux from ice bottom to top surface is decreased, it acts to increase the sea surface temperature because the energy staying in the bottom boundary of the sea ice is rising, which could cause the change in ocean-sea ice heat flux acting on changes in SIC (McPhee, 1992). The mechanical process related to changes in SIT is the ridging which is one of the dominant processes leading the ice deformation. The ridging process can bring out the ice floes to break into smaller blocks that accumulate, resulting in consumption of ice area (Martin, 2007). However, we leave it as a future work as accounting for the detail mechanism is outside the scope of this study.

Interannual Variability

In this sub-section, the interannual variability simulated in the reanalysis products will be evaluated. We first calculated the root-mean-square error (RMSE) of the SIE anomalies between the data assimilation experiments and the reference data (SSMIS NT, or the PIOMAS) with respect to the calendar month (Figure 9). Note that the monthly anomaly is defined by subtracting the climatology at the corresponding calendar month. The result in noDA have largest errors among all experiments in most of seasons, and those in boreal summer is particularly large (gray lines in Figure 9). In Ver1, the RMSE of SIE reduces for most of seasons compared to that in noDA especially during boreal summer and autumn (blue lines in Figure 9). The RMSE in Ver2 is in similar amplitude compared to Ver1, however, RMSE during boreal autumn and winter (October–March) is systematically higher (red lines in Figure 9). This indicates that, even though the assimilation of the SIT observation improves the simulation quality of the climatological SIE (as demonstrated in the previous sub-section), this advantage is not shown in simulating the interannual variability. Nonetheless, we will show that the quality in simulating the dominant SIC variability with spatial distribution is improved in Ver2 in next sub-section.

FIGURE 9

Figure 9. The root-mean-square-errors with respect to the SSMIS NT (solid lines) and PIOMAS (dashed lines) of sea ice extent index (10⁶km²) for noDA (gray), Ver1 (blue and skyblue), and Ver2 (red and pink) during 2000–2019 (A) and 2011–2019 (B).

To investigate the time evolution of the SIE in more detail, we calculated the time-series of the SIE anomaly in February and August in Figure 10. In February, the noDA follows the overall variability well due to influence of the prescribed atmospheric and SST forcings (correlation coefficient between the SIE in noDA and SSMIS NT is 0.77) (Figure 10A). As noDA already exhibited a good performance, the positive impact of data assimilation on the interannual variability of SIE in February is minimal.

FIGURE 10

Figure 10. The time-series of Arctic sea ice extent anomaly (10⁶km²) for SSMIS NT (black), noDA (gray), Ver1 (blue), and Ver2 (red) in February (A) and August (B) during 2000–2019. The “sm_R” in legend shows the correlation coefficient between each experiment and observation data.

However, in August, the quality of the reanalysis with the data assimilation is quite improved (Figure 10B). It is shown that the interannual variability of SIE in noDA is large in Figure 10B. We speculate that atmospheric forcing is the cause of large interannual variability in SIE in “noDA” in August. The downward shortwave radiation of the NCEP reanalysis data is relatively higher than that of other forcing data sets (e.g., JRA-55 and ERA-Interim) especially in summer, which may induce excessive melting of sea ice (Lee et al., 2019). As a result, it was understood that the thinning of sea ice due to the excessive shortwave radiation can lead to increased interannual variability of sea ice (Holland et al., 2008). The anomaly correlation coefficients between the noDA and the SSMIS NT is only 0.51, while the correlation with the SSMIS NT is improved to 0.94 (0.95) for Ver1 (Ver2).

In order to explore the regional features, we calculated the grid point-wise RMSE of the SIC anomalies during 2011–2019 in Figure 11. In February, the RMSEs of SIC in noDA, Ver1, and Ver2 against SSMIS NT tends to be similar from each other, which is consistent with the minimal improvement in simulating SIE in February. However, in August, the improvement by assimilating sea ice variables is clear. The RMSEs in noDA are systematically large over the overall Arctic regions, while that in Ver1 or Ver2 is significantly reduced in most regions (Figures 11D–F).

FIGURE 11

Figure 11. The root-mean-square-errors of sea ice concentration (a fraction from 0 to 1) for noDA, Ver1, and Ver2 with respect to the SSMIS NT observation are calculated in February (A–C) and in August (D–F). The RMSE of noDA (A,D) and the difference of Ver1 (B,E) and Ver2 (C,F) with respect to noDA are calculated during 2011–2019.

Figure 12 shows the RMSEs for the SIT anomalies in February and August in noDA, Ver1, and Ver2 against the CS2SMOS. In February, the noDA show that the RMSEs in CAP regions including east Greenland, Wandel Sea, and east of Beaufort Sea are relatively large than other regions (Figure 12A). In Ver1, the RMSE is slightly large over the polar regions compared to noDA (Figure 12B). In Ver2, the RMSE in CAP region is systematically reduced than that in noDA and Ver1 (Figure 12C). In August, the relatively large RMSEs are mainly shown over the CAP region in noDA, and it is relatively small in the remaining regions (Figure 12D). The RMSE of the SIT anomalies in Ver1 is increased to some extent compared to noDA (Figure 12E). Even though the RMSE of the SIT anomaly in Ver2 is systematically reduced than that in Ver1, however, it is only comparable to noDA (Figure 12F). This modest improvement is probably due to the fact that this is the season when the SIT observation is not injected even in Ver2 due to the lack of SIT observation.

FIGURE 12

Figure 12. The root-mean-square-errors of sea ice thickness (m) for noDA, Ver1, and Ver2 with respect to the CS2SMOS observation and PIOMAS reanalysis are calculated in February (A–C) and in August (D–F). Note that the CS2SMOS is used in February and the PIOMAS is used in August due to the absence of satellite data. The RMSE of noDA (A,D) and the difference of Ver1 (B,E) and Ver2 (C,F) with respect to noDA are calculated during 2011–2019.

Dominant Modes

To analyze the spatial distribution of the dominant interannual variability of the SIC and SIT, we performed the empirical orthogonal function (EOF) analysis for observation and reanalysis. The leading EOF of SIC anomalies for SSMIS NT, PIOMAS, Ver1, and Ver2 in February for 2000–2019 are shown in Figure 13. The leading EOF for SSMIS NT, PIOMAS, Ver1, and Ver2 is well-separated from its second EOF and their explained variance is 23.2, 26.1, 24.4, and 32.5%, respectively. The leading EOF of SIC anomalies for SSMIS NT exhibits dipole patterns in the Greenland-Barents Sea and Baffin Bay-Labrador Seas in Atlantic sector (Figure 13A). The winter Odden sea ice feature, which denotes a large positive tongue of SIC anomaly that advances rapidly northeastward into the Greenland Sea over 8°W−5°E, 73–77°N (Shuchman et al., 1998; Deser et al., 2000), appears in all four products, even though those in the Ver1 and Ver2 exhibited a larger amplitude than the SSMIS NT or PIOMAS. The negative SIC anomalies are shown in Baffin Bay and Labrador Sea in all four products.

FIGURE 13

Figure 13. The spatial distributions of the eigenvector of the first empirical orthogonal function mode of sea ice concentration (a fraction from 0 to 1) for SSMIS NT observation (A), PIOMAS reanalysis (B), Ver1 (C), and Ver2 (D) in February during the period 2000–2019. The pattern correlation represented “PCC” is calculated between SSMIS NT observation and each product. The time-series of the principal component (E) for SSMIS NT (black), PIOMAS (light green), Ver1 (blue), and Ver2 (red) is calculated and the “sm_R” represents the correlation coefficients with SSMIS NT data.

In the Pacific sector, the positive SIC anomalies are shown over the Bering Sea in all four products. In the SSMIS NT and the PIOMAS, this SIC anomalies over the Bering Sea is mostly confined at the north of 60°N, whereas, in the Ver1 and Ver2, it is further extended to the south. The positive SIC anomalies appear in the Sea of Okhotsk in SSMIS NT, PIOMAS, and Ver2, however, it is overwhelmed by the negative values in Ver1. As a result, the SIC anomalies of the dominant EOF in Ver1 exhibited the erroneous dipole pattern. This is also shown in the dominant EOF in noDA (not shown), implying that the dipole pattern can be corrected through SIT data assimilation. In order to examine the similarity with the leading EOF of the SSMIS NT, we calculate the pattern correlation coefficients (PCC) for PIOMAS, Ver1 and Ver2. The dominant EOF in the PIOMAS is most similar to that in the SSMIS NT (PCC = 0.84), followed by Ver2 (0.71) and Ver1 (0.54).

The standardized principal component (PC) time series of the leading EOF for each data is shown in Figure 13E. All PC time series show a decreasing trend. The temporal correlation of the PC time-series from the SSMIS NT is almost similar between noDA and Ver1; the temporal correlation coefficient between the SSMIS NT and noDA is 0.86, where that between SSMIS NT and Ver1 is 0.72. However, assimilating SIT observation in Ver2 further improves the simulation quality of PC time-series (0.89). This means that the indirect effect of SIT data assimilation on SIC simulation contributes to the realistic SIC interannual variations.

The assimilation of the SIT observation plays a critical role in mimicking the observed leading EOF and the associated PC time series in February (Figure 14). The explained variance of the first modes accounts for 28.5, 28.5, 25.5, and 32.4%, for CS2SMOS, PIOMAS, Ver1, and Ver2, respectively. The CS2SMOS and PIOMAS shows that the interannual variability is robust over the CAP region (Figure 14A). In the CS2SMOS, a tongue of positive SIT anomaly extending from west to east in Beaufort Sea between 120°W and 180°E along the latitude of 71–74°N, associated with the Beaufort Gyre, is shown. However, the dominant EOF of the SIT in Ver1 exhibited quite different spatial distribution from the observation (Figure 14C). There is an arbitrary negative, and positive SIT anomalies near the CAP region, and Wandel Sea in Ver1, respectively. As a result, the PCC between the leading EOF in Ver1 and that of CS2SMOS is −0.27, and the correlation of the corresponding PC time-series is only 0.06 (Figure 14E). This clearly indicates the necessity of the SIT data assimilation for the realistic simulation of the SIT anomalies.

FIGURE 14

Figure 14. The first EOF mode of sea ice thickness in February during the period 2011−2019. The CS2SMOS satellite thickness data is used here. The arrangement of the plots is the same as in Figure 13.

On the other hand, in Ver2, the quality in simulating SIT anomalies is systematically improved (Figure 14D). The amplitude of the positive SIT anomalies in CAP region in Ver2 is similar to those in CS2SMOS. In addition, the zonally-elongated positive SIT anomalies associated with Beaufort Gyre is successfully mimicked to some extent only in Ver2. The PCC between the dominant EOF of CS2SMOS and that of Ver2 is 0.80, which exhibited even higher PCC value between CS2SMOS and PIOMAS (i.e., 0.61). The corresponding PC time-series in Ver2 is also well-correlated with the CS2SMOS (i.e., temporal correlation is 0.89) (Figure 14E). It appears that data assimilation of SIT significantly improves the simulation quality of dominant Arctic SIT variability during boreal winter.

The reanalysis quality in simulating SIC anomalies during boreal summer is quite similar between SSMIS NT, PIOMAS, Ver1, and Ver2 (not shown). On the other hand, the dominant EOF of the SIT anomalies in Ver2 exhibited a modest improvement during boreal summer (Figure 15). The leading EOF in PIOMAS shows the Central Arctic Thinning-like pattern (Fučkar et al., 2016) and the PCC between the dominant EOF of PIOMAS and that of Ver2 is 0.82, which exhibited slightly higher PCC between PIOMAS and Ver1 (i.e., 0.74). The corresponding PC time-series in Ver2 also exhibited higher temporal correlation with the PIOMAS (i.e., correlation is 0.87) than that in Ver1 (i.e., correlation is 0.72) (Figure 15D). As there is no observation of SIT in summer, the improvement in simulating SIT anomalies during boreal summer in Ver2 might be originated from the improved reanalysis quality of the SIT anomalies during boreal winter. That is, the realistic simulation of the SIT anomalies during boreal winter season is maintained up to boreal summer season possibly due to the multi-season persistency of the sea ice anomalies (Holland et al., 2011).

FIGURE 15

Figure 15. The spatial distribution of the eigenvector of first EOF mode of sea ice thickness for PIOMAS (A), Ver1 (B), and Ver2 (C) in August during 2011–2019. The pattern correlation represented “PCC” is calculated between PIOMAS and each experiment. The time-series of the principal component (D) for PIOMAS (light green), Ver1 (blue), and Ver2 (red) is calculated and the “pi_R” represents the correlation coefficients with PIOMAS reanalysis data.

Summary and Discussions

In this work, we have exploited the satellite observation of SIT from the CryoSat-2 and SMOS from European Space Agency's (ESA), and SIC from the Special Sensor Microwave Imager/Sounder (SSMIS) aboard Defense Meteorological Satellite (DMSP) to assimilate the SIC and SIT in CICE5 model using ensemble optimal interpolation (EnOI). The SIC, and the SIT satellite observation are assimilated during 2000–2019, and 2011–2019, respectively. In order to use the EnOI, the stationary ensemble perturbations for calculating background error covariance matrix are generated by using long-term model integrations during 1982–2019. We have found that the stationary ensemble perturbations not only represent the seasonal variability of sea ice variable, but also successfully reflect the physical properties and balances related to the sea ice variability.

The quality of the reanalysis produced by sea ice data assimilation are evaluated in terms of climatological field, interannual variability, and dominant interannual variation by comparing assimilated states and the those in the observations and reanalysis. By assimilating the SIC satellite observations, the systematic bias of SIC in control simulation (i.e., noDA) is remarkably corrected in marginal ice zones, and pan-Arctic ocean during boreal winter, and summer season, respectively (Figure 6). The sea ice extent (SIE) anomalies in data assimilation experiments exhibited a similar year-to-year variation to SSMIS NT compared to that in noDA in all calendar months especially during boreal summer (Figures 9, 10). The SIC assimilation have an effect on the correction in climatological bias of SIT, but interannual variations of the SIT are still unrealistic (Figure 7), which requires an additional assimilation of SIT observation.

Not only the climatological bias, but also the RMSEs of SIT anomalies are further reduced particularly during boreal winter by additionally assimilating SIT observations (Figure 12). The empirical orthogonal function (EOF) analysis showed that the dominant interannual variation of SIT anomalies exhibited similar spatial distribution to those in the observation during boreal winter season (Figure 14). The positive SIT anomalies observed in central Arctic ocean, centered over the Canadian Archipelago (CAP), are well-simulated in Ver2 (experiment with SIC and SIT assimilation), where the negative anomalies are distributed in CAP and Beaufort Sea in Ver1 (experiment with SIC assimilation).

Interestingly, the improvement in the SIT states leads to the correction in the SIE bias (Figure 5). The reduction in the negative SIT bias by additionally assimilating SIT observation contributes to reduce the positive SIC bias over Kara-Barents Sea (Figure 8). This inverse relationship between SIC and SIT is also shown in the idealized experiment by increasing the SIT. Even though the assimilation of SIT observation did not result in appreciable differences in interannual variability of SIE anomalies, there has been the improvement in simulating dominant interannual variation for SIC in February (Figure 13). The observed dipole pattern between Greenland-Barents Sea and Baffin Bay-Labrador Sea, and the zonally-elongated positive SIC anomalies, referred as Odden feature, are well-simulated in reanalysis produced by both in Ver1 and Ver2. However, the observed positive anomalies in Sea of Okhotsk are successfully reproduced only in Ver2. The PC time-series of the leading EOF for SIC in Ver2 showed a higher temporal correlation with the observation than that in Ver1.

Through the series of the data assimilation of the sea ice variables, not only the degree of the improvement in simulating the climatology and the interannual variability of the sea ice can be assessed, but also a scientific question can be drawn. The indirect impact of the SIT assimilation on the SIC is one of typical subjects. The simulation of the SIC climatology and the variability tends to be improved once the SIT observation is further assimilated in addition to the SIC observation. This indicates that the non-linear physical processes associated with the corrected SIT help to improve the SIC simulation (Collow et al., 2015). In addition, even though the SIC improvement by assimilating SIT observation is not clear during boreal summer season in this study, several previous studies showed that the knowledge of winter ice thickness can provide some predictive capability for summer ice extent (Kauker et al., 2009; Holland et al., 2011). This means that the reduction the systematic errors in SIT during boreal winter season have a seasonal memory to improve the SIC simulation during boreal summer. To fully understand the mechanism of the improvement led by the data assimilation, the interactions between sea ice variables, and temporal memory of the sea ice variables, which might be dependent on the season, should be investigated further as a future work.

In addition, the improvement of the SIC simulation by assimilating SIT observation implies that the multi-variate assimilation might further improve the reanalysis quality. As it is already demonstrated that the updated SIT leads a positive impact on the SIC through physical processes in the model, the direct updates of the SIC by given SIT increment based on the modeled covariance would lead the improvement in simulating SIC without any time lags.

The reanalysis produced in this study aims to improve the quality of the initial condition of the sea ice component in the Global Seasonal forecast system (GloSea) for the seasonal forecasts (Maclachlan et al., 2015). In particular, as still operational version of the GloSea (i.e., GloSea version 5) for the seasonal forecasts only assimilates the SIC, the improvement in the sea ice initial conditions by assimilating the SIT observation is expected. This expectation is based on the previous studies that a significant improvement in predictive skill of Arctic sea-ice extent and ice edge location for multi-seasons forecasts targeted to September Arctic sea ice by assimilating the CryoSat-2 thickness (Blockley and Andrew Peterson, 2018). However, the temporal coverage of the SIT observation is limited after 2010, and this significantly shortens the hindcasts period compared to the conventional period for the hindcasts for seasonal forecasts (i.e., longer than two decades, e.g., Johnson et al., 2019). Nonetheless, the temporal coverage of a decade (i.e., from 2010 to the present) would allow the preliminary assessment of the seasonal forecast skill and its possible advantages by assimilating the SIT observations.

Data Availability Statement

The original contributions presented in the study are included in the article, further inquiries can be directed to the corresponding author.

Author Contributions

J-GL and Y-GH designed the research and wrote the paper. J-GL performed the experiments and analyzed the data. Both authors discussed the methodology, results and structure of the paper. Both authors contributed to the article and approved the submitted version.

Funding

This study was supported by the Korea Meteorological Administration Research and Development Program under Grant KMI2021-01210.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We thank the NSIDC for providing SSMIS sea ice concentration (http://nsidc.org/data/nsidc-0051 and nsidc-0079), CryoSat-2 sea ice thickness (http://nsidc.org/data/RDEFT4), the University of Hamburg for providing SMOS sea ice thickness (https://icdc.cen.uni-hamburg.de/thredds/catalog/ftpthredds/smos_sea_ice_thickness/catalog.html), and the Alfred Wegener Institute for sharing CS2SMOS sea ice thickness (https://spaces.awi.de/display/CS2SMOS/). NCEP_Reanalysis 2 and OISST V2 data were provided by the National Oceanic and Atmospheric Administration (NOAA) and the PIOMAS reanalysis is provided by Polar Science Center.

References

Aksenov, Y., Popova, E. E., Yool, A., Nurser, A. J. G., Williams, T. D., Bertino, L., et al. (2017). On the future navigability of Arctic sea routes: high-resolution projections of the Arctic Ocean and sea ice. Mar. Policy 75, 300–317. doi: 10.1016/j.marpol.2015.12.027