Seasonal Arctic Sea Ice Prediction Using a Newly Developed Fully Coupled Regional Model With the Assimilation of Satellite Sea Ice Observations

To increase our capability to predict Arctic sea ice and climate, we have developed a coupled atmosphere‐sea ice‐ocean model configured for the pan‐Arctic with sufficient flexibility. The Los Alamos Sea Ice Model is coupled with the Weather Research and Forecasting Model and the Regional Ocean Modeling System in the Coupled Ocean‐Atmosphere‐Wave‐Sediment Transport modeling system. It is well known that dynamic models used to predict Arctic sea ice at short‐term periods strongly depend on model initial conditions. Parallel Data Assimilation Framework is implemented into the new modeling system to assimilate sea ice observations and generate skillful model initialization, which aid in the prediction procedures. The Special Sensor Microwave Imager/Sounder sea ice concentration, the CyroSat‐2, and Soil Moisture and Ocean Salinity sea ice thickness are assimilated with the localized error subspace transform ensemble Kalman filter. We conduct Arctic sea ice prediction for the melting seasons of 2017 and 2018. Predictions with improved initial sea ice conditions show reasonable sea ice evolution and small biases in the minimum sea ice extent, although the ice refreezing is delayed. Our prediction experiments suggest that the use of appropriate uncertainty for the observed sea ice thickness can lead to improved spatial distribution of the initial ice thickness and thus the predicted sea ice distribution. Our new modeling system initialized by the output of the National Centers for Environmental Prediction Climate Forecast System seasonal forecasts with data assimilation can significantly increase the sea ice prediction skills in sea ice extent for the entire Arctic as well as in the Northern Sea Route compared with the predictions by the National Centers for Environmental Prediction Climate Forecast System.


Introduction
Associated with global climate change, the properties of Arctic sea ice have undergone a dramatic change over the past few decades. Arctic sea ice extent/area has been decreasing since the satellite era. The decreasing trend in areal extent is found for every month (e.g., Cavalieri & Parkinson, 2012;Liu et al., 2019). Meanwhile, the thickness of the ice pack has been thinning (e.g., Kwok & Cunningham, 2015;Kwok & Rothrock, 2009;Kwok & Untersteiner, 2011), largely due to the thicker multiyear sea ice being replaced by the thinner first year sea ice (e.g., Maslanik et al., 2007Maslanik et al., , 2011Tschudi et al., 2016).
The rapid changes in the properties of Arctic sea ice have brought challenges to an increasing demand for skillful Arctic sea ice predictions, particularly at a seasonal timescale (e.g., Jung et al., 2013;Stroeve et al., 2014). Figure 1 gives the predicted September sea ice extent based on the July Sea Ice Outlook from 2008 to 2018, a community effort managed by the Sea Ice Prediction Network (https://www.arcus.org/sipn). Here we only show the prediction from a range of dynamical models participating in Sea Ice Prediction Network, including sea ice-ocean coupled models forced by prescribed atmospheric forcings from numerical weather predictions and fully coupled (including atmosphere, sea ice, ocean, and land) regional and global modeling systems. It appears that the median of the predicted September sea ice extent shows substantially differences compared to the observations for most years. Thus, more efforts are needed to improve the prediction capability for Arctic sea ice and the climate system at a seasonal timescale.
Many models have been used for Arctic research. For example, The Arctic Regional Climate Model Intercomparison Project (ARCMIP, Curry & Lynch, 2002) evaluated the performance of atmospheric models. Studies that compared ARCMIP models with Surface Heat Budget of the Arctic Ocean observations (Uttal et al., 2002) showed that there was a large disagreement between ARCMIP models and observations in surface radiation fluxes, which is mainly the result of errors in cloud cover and vertical cloud distribution (e.g., Inoue et al., 2006;Rinke et al., 2006;Tjernstrőm et al., 2008;Wyser et al., 2008). More realistic physical parameterizations are needed to reproduce reasonable results. The Arctic Ocean Model Intercomparison Project (Proshutinsky et al., 2001) assessed the performance of ocean models. Studies have reported that Arctic Ocean Model Intercomparison Project models showed large differences in temperature, salinity, and circulations as compared to observations (e.g., Holloway et al., 2007;Proshutinsky et al., 2011). The Sea Ice Model Intercomparison Project (Notz et al., 2016) has been endorsed as a diagnostic MIP for CMIP6, which will analyze heat, momentum, and mass budgets that govern the evolution of sea ice. Sea Ice Model Intercomparison Project will allow researchers to understand the role of internal variability, external forcing, model tuning, and the formulation of sea ice models for improving the fidelity in simulating and predicting sea ice.
It is well known that the performance of dynamic models to predict Arctic sea ice at short-term periods strongly depends on model initial conditions (e.g., Blanchard-Wrigglesworth et al., 2011;Msadek et al., 2014;Liu et al., 2019;Yang et al., 2016). The accurate sea ice initialization requires not only sea ice concentration but also variables (e.g., sea ice thickness) that influence surface energy fluxes and ocean-atmosphere interactions (Yang et al., 2016). Thus, a data assimilation method is needed to incorporate sea ice observations to produce skillful model initializations, leading to improved Arctic sea ice prediction. Several studies have explored assimilating observed sea ice concentration into coupled sea ice-ocean models. For example, Lisaeter et al. (2003) used an ensemble Kalman filter (KF) to assimilate Special Sensor Microwave/Imager ice concentration, while Lindsay and Zhang (2006) employed a nudging scheme. Stark et al. (2008) used an optimal interpolation method to assimilate the Special Sensor Microwave/Imager ice concentration as well as sea ice drift data. Tietsche et al. (2013) assimilated ice concentration with a simple Newtonian relaxation scheme and updated the sea ice thickness analysis using a proportional dependence between the ice concentration and mean ice thickness in a coupled climate model. These studies suggested that the assimilation of observed ice concentration in the models improves the simulated sea ice cover but minor improvements for the sea ice thickness.
At a seasonal timescale, the initialization of sea ice thickness has been shown to be critical for summer prediction (e.g., . Some studies (e.g., Blanchard-Wrigglesworth et al., 2011;Koenigk & Mikolajewicz, 2009) showed that the persistence of sea ice thickness anomalies is much longer than that of sea ice concentration. Higher predictability of Arctic sea ice thickness (volume) with respect to that of Arctic sea ice cover has been found at longer time scales (e.g., Guemas et al., 2014). Stroeve et al. (2014) suggested that the quality of initial states and the initialization methods that the modeling groups utilize are the key to close the gap between the operational skill and theoretical predictability. Recent studies have shown that assimilating observed sea ice thickness and/or ice concentration can lead to significantly reduced biases in predicted sea ice extent and thickness based on the simulation in the National Centers for Environmental Prediction (NCEP) Climate Forecast System (Chen et al., 2017;Collow et al., 2015). Benefiting from both the CryoSat-2 and Soil Moisture and Ocean Salinity (SMOS) satellite-derived sea ice thickness products, more recent studies have shown the improvements on short-term ice predictions as well as the ice volume estimation by assimilating CyroSat-2 and/or SMOS together with satellite-derived sea ice concentration (Blockley & Peterson, 2018;Xie et al., 2018).
In this study, we have developed a new coupled modeling system configured for the Arctic that can be used for improving predictive capability of Arctic sea ice and the climate system at a seasonal time scale, as well as resolving and understanding important processes and feedbacks in the Arctic. The goal of the new coupled modeling system is sufficient flexibility. Sufficient flexibility means that each model component of the modeling system has different physics options for us to choose (see section 2 for details). For improving the accuracy of initial conditions, we have incorporated a newly developed ensemble KF into the new coupled modeling system. We use this filter to assimilate not only satellite-based sea ice concentration but also recently produced and validated satellite-based sea ice thickness, which provides improved sea ice initial conditions for the model (see section 3 for details).

Model
As mentioned above, we have developed a new coupled modeling system configured for the Arctic with sufficient flexibility. Specifically, the Coupled Ocean-Atmosphere-Wave-Sediment Transport (COAWST) modeling system (Warner et al., 2010 is employed here as the base framework for our regional coupled modeling system. The COAWST modeling system includes three model components: the Weather Research and Forecasting Model (WRF), the Regional Ocean Modeling System (ROMS), and the Simulating Waves Nearshore. The COAWST is designed and used to better understand significant processes affecting our coastlines and how those processes create coastal changes (Warner et al., 2010).
For Arctic sea ice prediction at a seasonal time scale, we couple the Los Alamos sea ice model (CICE) with the WRF model and the ROMS model using the Model Coupling Toolkit in the COAWST modeling system. More details about each model component are described below.

Atmospheric Model Component
The atmospheric model component in the COAWST modeling system is the Advanced Research Weather Research and Forecasting model (ARW, Skamarock et al., 2005). The ARW model is a nonhydrostatic and quasi-compressible model, which uses the Arakawa C-grid in horizontal and sigma-pressure vertical coordinate in vertical with the top of the model at 50 mb. The ARW model includes recent advances in physics parameterizations, that is, P3 microphysics scheme (Morrison & Milbrandt, 2015) and KSAS cumulus scheme (Kwon & Hong, 2017). The WRF model has been widely used in operational weather forecast and atmospheric research (e.g., Bromwich et al., 2009;Done et al., 2004;Hines & Bromwich, 2008;Hines et al., 2011;Liu et al., 2016).

Oceanic Model Component
The COAWST modeling system uses the ROMS as the ocean model component. The ROMS model is a free-surface and terrain-following model, which solves three-dimensional Reynolds-averaged Navier-Stokes equations with the hydrostatic and Boussinesq approximation (Haidvogel et al., 2008;Shchepetkin & McWilliams, 2005). In the horizontal, the ROMS model uses boundary-fitted, orthogonal curvilinear coordinates on a staggered Arakawa C-grid. In the vertical, the equations are discretized over bottom topography with stretching terrain-following coordinates (Song & Haidvogel, 1994). Like the WRF model, the ROMS model provides various options for advection schemes (momentum and tracer), vertical mixing, subgrid-scale parameterizations, and open boundary conditions. The ROMS model has been used in a dives range of applications (e.g., Capet et al., 2008;Fennel et al., 2006;Siegel et al., 2008).

Sea Ice Model Component
The ROMS model provides a simple sea ice module, which combines the elastic-viscous-plastic rheology (Hunke & Dukowicz, 1997) and a simple one-layer ice and snow thermodynamics (Mellor & Kantha, 1989). Here, the Los Alamos sea ice model (CICE) version 5.1.2 is chosen to be coupled with the WRF model and the ROMS model within the COAWST modeling system. The CICE model is designed to be a sea ice component for global coupled climate models. Its dynamic part simulates the movement of sea ice based on the material strength of the ice and forces from the atmosphere, the ocean, and Earth's rotation; its thermodynamic part simulates the growth rate of snow and ice as a result of vertical conductive, radiative, and turbulent fluxes; and its transport part computes the advection of sea ice concentration and volume as well as other state variables. The CICE model also simulates the evolution of the ice thickness distribution (ITD, Thorndike et al., 1975) in time and space since the dynamic and thermodynamic properties of the ice pack are related to how much ice lies in each thickness range. More details can be found in the CICE documentation (Hunke et al., 2015).
The choices of physics in the WRF, ROMS, and CICE models determine the behavior of the simulated atmospheric (i.e., winds, radiation, temperature, humidity, and precipitation), oceanic (i.e., currents, temperature, and salinity), and sea ice (i.e., velocity, ice concentration, ice thickness, and melt pond distribution) variables that have profound influences on simulating sea ice conditions in the Arctic. Thus, a series of sensitivity experiments with different physics options have been performed to determine the "optimal" physics configuration that provides reasonable simulation of Arctic sea ice, serving as the baseline.

Assimilation Procedure
In this study, the Parallel Data Assimilation Framework (PDAF, Nerger & Hiller, 2013), an open-source software for ensemble data assimilation, is implemented into the new coupled modeling system described in section 2. The advantages of PDAF are as follows: (1) PDAF simplifies the implementation of the data assimilation system with numerical models and provides a standardized interface that makes the assimilation routines independent from the numerical models and (2) PDAF provides a variety of optimized data assimilation algorithms, in particular ensemble-based KFs. The KF (Kalman, 1960;Kalman & Bucy, 1961) can be mathematically expressed as follows.
Here x f and x a are the forecast and analysis state vectors. P f and P a are the model forecast and analysis error covariance matrices. K is the Kalman gain. y and R are the observation vector and its error covariance matrix. H is the forward operator, performing the mapping from model space to observation space. A forecast model M is applied to transport the state vector and the model error covariance matrix, with the covariance of forecast model error Q, to the next analysis time step: In contrast to the KF that works with the full distribution of the state explicitly, the ensemble Kalman filter (Evensen, 1994), a flow-dependent data assimilation method, was proposed, which combines the KF with a sequential Monto Carlo method, and uses an ensemble to approximate the covariance of the background error. That is, one never needs to compute a full covariance matrix in model state space in the ensemble Kalman filter approach. In addition, it has an advantage of relative ease of implementation and restoring the ensembles, retaining the benefit of flow-dependent error covariance.
Here, a localized error subspace transform ensemble Kalman filter (LESTKF) is utilized to assimilate satellite-observed sea ice data. Nerger et al. (2012a) showed that the error subspace transform ensemble Kalman filter (ESTKF), as a variant of the singular evolutive interpolated Kalman (Pham, 2001) filter, produces the same ensemble transformation as that of the ensemble transform Kalman filter (ETKF, Bishop et al., 2001). The ESTKF explicitly projects the ensemble onto the error subspace and computes the ensemble transformation directly in the error subspace. The ESTKF also provides minimum transformation of the ensemble members, which performs better than singular evolutive interpolated Kalman, and has lower computational cost. To minimize sampling errors of the estimated error covariance matrix because of small ensembles, especially long-range error covariance matrix, the localization of covariance matrices is applied. The localized ESTKF conducts a regulated observation localization by the weighting of elements of the observation-error covariance matrix with a localization function of variable width (Nerger et al., 2012b). For smaller observation errors, the regulated observation localization prevents the broadening of the effective localization length that can degrade the assimilation performance. Furthermore, its extra computational cost is insignificant compared with the total cost (see details in Nerger et al., 2012b).
Our assimilation procedure includes two steps: initialization and analysis. Using the prediction in this study that starts on 1 July 2017 as an example, for the initialization step, we follow the ensemble-generating method proposed by Pham (2001) to perform 1 month of simulation with the new modeling system (starting on 15 June and on 15 July), which is initialized by the NCEP 9-month operational forecasts starting on 15 June. The daily simulated state vectors of sea ice concentration and thickness are stored sequentially to form the trajectory of sea ice states. The second-order exact sampling is then applied to the singular vectors and values, which are computed by performing the singular value decomposition on the trajectory of sea ice states, to generate eight ensemble perturbations overlaid on the mean state of trajectories. Once we have the initial ensembles that provide an estimate of initial model states and the model uncertainties, the analysis step with the LESTKF is conducted to assimilate satellite-based sea ice concentration and thickness with their observational uncertainties into the initial sea ice conditions (see section 4.1). As described above, the CICE model uses categories to describe the ITD evolution, and thus, a method to distribute the assimilated sea ice concentration and thickness to the ITD space is required. First, the assimilated sea ice concentration and thickness are treated as the mean value of grid cells. Second, we employ the module in the CICE model to distribute the initial mean value to each ice thickness category using a distribution function, P: where A and h are the assimilated sea ice concentration and thickness, H i is the prescribed ice thickness category (0-0.3, 0.3-0.7, 0.7-1.2, 1.2-2.0, and 2.0 m~), N is the number of ice thickness category (here N = 5), and A i and V i are the ice concentration and volume of each ice thickness category.

Initialized and Assimilated Data
The Climate Forecast System version 2 (CFSv2, Saha et al., 2014) 9-month operational forecast is obtained from the NCEP (http://nomads.ncep.noaa.gov/pub/data/nccf/com/cfs/prod/). The CFSv2 forecast provides all variables with a spatial resolution of 1°for initializing all model components of our new modeling system.
We assimilate two satellite-observed sea ice variables (sea ice concentration and thickness) to generate improved estimation of initial sea ice conditions. The near real-time daily Arctic sea ice concentration is obtained from the National Snow and Ice Data Center (NSIDC), which is derived from the Special Sensor Microwave Imager/Sounder (SSMIS), processed by the the National Aeronautics and Space Administration Team algorithm (Maslanik & Stroeve, 1999). The data are on the NSIDC polar stereographic grid with a spatial resolution of 25 km. The uncertainty of near real-time SSMIS sea ice concentration is 10-15% (Tonboe & Nielsen, 2010), and 15% is used as the observational error in sea ice concentration.
One of the most significant biases in the NCEP CFSv2 simulation is sea ice being too thick over most of the Arctic basin with interannual variability that is too weak (Saha et al., 2014). Two satellite-derived sea ice thickness products are employed here. The first product is the monthly sea ice thickness derived from the CryoSat-2 (Laxon et al., 2013), a satellite of the European Space Agency. The CryoSat-2 sea ice thickness (http://data.seaiceportal.de) is on the EASE Grid (Brodzik et al., 2012) with a spatial resolution of 25 km. The uncertainty of CryoSat-2 ice thickness comes from two sources, the freeboard-to-thickness (~0.3 m) and its system uncertainty (~0.6-1.2 m) (Ricker et al., 2014). The CryoSat-2-derived sea ice thickness is more accurate for the central Arctic ice pack. It is also possible for thin ice measurement but with less robustness (Sallila et al., 2019). The second product is the daily sea ice thickness derived from the SMOS brightness temperature measured by the Microwave Imaging Radiometer using Aperture Synthesis using a single-layer emissivity model (Kaleschke et al., 2012;Tian-Kunze et al., 2014; https://icdc.cen.uni-hamburg.de/ thredds/catalog/ftpthredds/smos_sea_ice_thickness/catalog.html). The SMOS sea ice thickness is on the NSIDC polar stereographic grid with a spatial resolution of 12.5 km. The SMOS ice thickness has small uncertainty (~0.7 m) for thin ice (<1 m) but large underestimation for thick ice because of the saturation of SMOS brightness temperature with thickness (Kaleschke et al., 2012;Tian-Kunze et al., 2014). In order to take advantage of the complementary feature of CryoSat-2 and SMOS, we use a simple approach to merge the two ice thickness products. Sea ice with thickness less than 1 m in the CryoSat-2 ice thickness is replaced by SMOS ice thickness. The 1.5 m is used here to represent the observational error in the CryoSat-2/SMOS merged sea ice thickness, which is the largest uncertainty of CryoSat-2.
In this study, we only consider the assimilation of satellite-based sea ice concentration and thickness. The variables of the atmosphere and the ocean are initialized by the CFSv2 operational forecasts (see section 4.1) without assimilating additional observations for two reasons. First, our focus is the impacts of assimilation of key sea ice variables on Arctic sea ice prediction. Second, the CFSv2 data that are used as the initial and boundary conditions for our modeling system has assimilated the variables of the atmosphere and ocean (Saha et al., 2010).

Evaluation Data
In a fully coupled model, sea ice strongly interacts with both the atmosphere above and the ocean below. In this study, we also assess the simulation of the WRF model and the ROMS model in the new coupled predictive model. For the evaluation of the WRF model, the ERA-Interim Re-Analysis (Dee et al., 2011) is chosen. Lindsay et al. (2014) suggested that ERA-Interim has the better performance in the Arctic than other global reanalysis. For the evaluation of the ROMS model, the National Oceanic and Atmospheric Administration (NOAA) Optimum Interpolation Sea Surface Temperature (SST) version 2 (obtained from NOAA/OAR/ESRL PSD, Boulder, Colorado, USA, http://www.esrl.noaa.gov/psd/) with a spatial resolution of 0.25°is utilized.

Model Configurations and Experiment Designs
Currently, seasonal Arctic sea ice prediction typically starts from 1 June, 1 July, and 1 August (see Sea Ice Prediction Network for details, https://www.arcus.org/sipn). In this study, numerical experiments for the 10.1029/2019MS001938

Journal of Advances in Modeling Earth Systems
Pan-Arctic seasonal sea ice prediction have been conducted, starting from 1 July and 1 August to 1 October for the year of 2017 and 2018, respectively. The model domain includes 319 (449) x-(y-) grid points with añ 25-km grid spacing for all model components (WRF, ROMS, and CICE, Figure 2). The WRF model uses 50 vertical levels, the ROMS model uses 40 vertical levels, and the CICE model uses 7 ice layers, 1 snow layer, and 5 categories of sea ice thickness. The coupling frequency across all model components is 30 min. Detailed options of physical parameterizations and model settings for the WRF, ROMS, and CICE models are summarized in Table 1. Initial and boundary conditions for the WRF and ROMS models are from the CFSv2 9-month operational forecasts. Sea ice initial conditions are generated from the assimilation method mentioned in section 3. Ensemble forecasts with eight members are conducted. Table 2 provides Arctic sea ice prediction experiments conducted in this study. These numerical experiments allow us to examine impacts of the assimilation of satellite-based sea ice concentration and thickness on seasonal sea ice prediction, including different approaches to generate initial sea ice conditions and different observational errors of sea ice thickness.

Melting Season Sea Ice Thickness
Sea ice thickness derived from CryoSat-2 and SMOS as mentioned in section 4.1 is not available during the melting season (only available until 15 April). This is primarily due to that sea ice is highly inhomogeneous during the melting season (i.e., a footprint of SMOS can be a mixture of wet snow/sea ice, melt ponds, and open water, Huntemann et al., 2014). Thus, we need an approach to extrapolate the CryoSat-2/SMOS merged ice thickness to the above specific dates (i.e., 1 July and 1 August). Here the Pan-Arctic Ice Ocean Modeling and Assimilation System (PIOMAS) daily sea ice thickness data (Zhang & Rothrock, 2003; http://psc.apl.uw.edu/research/projects/arctic-sea-ice-volume-anomaly/data/) are used for this purpose. Previous studies (e.g., Laxon et al., 2013;Schweiger et al., 2011) have showed that the PIOMAS sea ice thickness is in reasonable agreement with the observations, validated by the observations including submarines, moorings, ICESat (Schutz et al., 2005), and CryoSat-2 ice thickness. We use the PIOMAS daily sea ice thickness during 2000-2016 to calculate the climatology of seasonal evolution of sea ice thickness and then apply the seasonal cycle to the CryoSat-2/SMOS merged ice thickness to approximate ice thickness at the date of the model initialization. We try two approaches to do the extrapolation of sea ice thickness based on the calculation of the ratio of the ice thickness change between the initial date of predictions and the last available date of CryoSat-2 and SMOS data. The first approach is that all grids use the Arctic-averaged PIOMAS ice thickness to calculate seasonal cycle (Figure 3a), and the other approach is that each grid uses its PIOMAS ice thickness time series as the seasonal cycle (Figure 3b). Figure 4 shows the satellite-based sea ice concentration and ensemble mean of initial ice concentration without and with data assimilation for the JUL17 experiment. Before the data assimilation, the initial ice concentration, which is the average of one-month model simulation from 15 June to 15 July 15, exhibits quite smooth spatial distribution in the Arctic Basin with larger concentration in the area covered by the first year sea ice compared to the observations (Figure 4a vs. Figure 4b). After the data assimilation, the initial ice concentration has a spatial pattern in the Arctic Basin that is similar to the observations; that is, some clusters with the relatively low concentration (65-85%) can be found between the Chukchi, Lincoln, and Laptev Seas for JUL17 (Figure 4b vs. Figure 4c). The initial ice concentration at Hudson and Baffin Bays is also reduced

Journal of Advances in Modeling Earth Systems
after the assimilation (Figure 4b vs. Figure 4c). This is also the case for the AUG17 experiment (Figure 4e vs. Figure 4f). Figure 5 shows the satellite-based sea ice thickness and ensemble mean of initial ice thickness without and with assimilation for the JUL17 experiment. As mentioned above, the most significant bias of the NCEP CFSv2 is that sea ice is too thick over much of the Arctic basin with much weaker interannual variability (Saha et al., 2014). It is apparent that the initial ice thickness generated by the CFSv2 forecast is much thicker in the Beaufort Gyre compared with the observations (Figure 5a vs. (Figure 5b). After the assimilation, the thick bias in the Beaufort Gyre is reduced significantly, and the spatial distribution is more comparable with the CryoSat-2/SMOS/PIOMAS blended ice thickness, in which the thickest ice is located in the Canadian Archipelago and northern Greenland Sea (Figure 5a vs. Figure 5c). This is also the case for the AUG17 experiment (Figures 5d-5f)).

Prediction Experiments
The temporal evolution of the ensemble mean and ensemble spread of the predicted Arctic sea ice extent for the JUL17 experiment is shown in Figure 6. The ice extent is calculated as the sum of the area of each grid cells with ice concentration greater than 15%. Along with the total ice extent, we also calculate ice extent for the following subregions: (1) Beaufort and Chukchi Seas (120-180 W, 65-80 N); (2) East Siberian and Laptev Seas (90-180°E, 65-80 N); (3) Barents,Kara,and Greenland Seas (30 W to 90°E,[65][66][67][68][69][70][71][72][73][74][75][76][77][78][79][80]; and (4) Canadian Archipelago and Baffin Bay (30-120 W, 65-80 N). The ensemble mean (thick blue line) has añ 0.4 × 10 6 -km 2 bias at the initial, but the observed ice extent still lies within the ensemble spread. As a result, the predicted ice extent overestimates the observation for about 1 week. For the rest of July, the ensemble mean of the ice extent is in good agreement with the observations. For August, our model slightly overestimates the actual ice extent. For September, the predicted ice extent shows a mixed feature compared to the observations with a delayed ice re-freezing in late September. The observed minimum sea ice extent (4.64 × 10 6 km 2 ) occurred on 13 September 2017, whereas our predicted minimum (4.65× 10 6 km 2 ) occurred on 24 September 2017. At the regional scale, the predicted sea ice extent shows reasonable temporal evolution in the Beaufort-Chukchi Seas, the Barents-Kara-Greenland Seas, and the Baffin Bay-Canadian Archipelago, but large difference is identified in the East Siberian-Laptev Seas, which overestimates the observed ice extent.  Topographic formulation (Flocco & Feltham, 2007) 10.1029/2019MS001938

Journal of Advances in Modeling Earth Systems
As the benchmark for our predictions, we also present the climatology prediction (CLIM, based on the data during 1997-2016) and the damped anomaly persistence prediction (DAMP; Van den Dool, 2006). The damped anomaly persistence prediction is calculated as follows: is the corresponding sea ice extent from the 1997-2016 climatology, SIE anom (0) is the ice

10.1029/2019MS001938
Journal of Advances in Modeling Earth Systems extent anomaly at lag 0, A(t) is the autocorrelation of SIE in the 1997-2016 period, and σ(t) is the standard deviation at lag t. Compared with the CLIM/DAMP predictions, our 2017 model predictions are much better than the one based solely on the climatology and show similar temporal evolution compared with the DAMP predictions before mid-August. After mid-August, the DAMP predictions tend to have faster ice melting rate than our predictions and thus underestimate the seasonal minimum ice extent. Figure 7 shows the spatial distribution of the observed and the difference between the predicted sea ice concentration and the observations for all grid cells that both the predictions and the observations have at least 15% concentration for the JUL17 experiment. The vertical and horizontal lining areas represent difference of the ice edge location. For the significance test, we used the daily time series of each grid cell to conduct the t test for the 1:1 regression line based on the predicted values and the observed values (see Piñeiro et al., 2008 for details). The areas marked with cross represent the predicted values that are in good agreement with the observed values, passing 95% significance level. The predicted ice distribution is very similar to the observations in July (Figures 7a and 7d)). Although the predicted ice distribution in August and September is broadly consistent with the observations, the prediction shows slower retreat of sea ice in the Pacific side of the Arctic Ocean compared with the observations, especially in the Chukchi and East Siberian Seas (Figures 7b, 7c, 7e, and 7f).  Figure 8 shows the ensemble mean of simulated sea ice thickness of the JUL17 experiment. For July (Figure 8a), the ice thickness in the northern Beaufort, Chukchi, and East Siberian Seas is~1-1.8 m thick.

Journal of Advances in Modeling Earth Systems
In the subsequent months of prediction (Figures 8b and 8c), a major thinning of the ice occurs in the northern East Siberian Sea from 1.6 m to 0.6 m, whereas the ice thickness in the northern Beaufort and Chukchi Seas does not decrease significantly, still thicker than 1.2 m. Such spatial evolution of ice thickness contributes to the slow ice retreating along the Pacific side of Arctic Ocean in JUL17, which is also in part associated with the initial sea ice thickness. That is, the sea ice thickness assimilation (Figure 5c) generates thinner ice compared to that of the CFSv2 in the Pacific side of Arctic Ocean, especially in the Chukchi Sea, but still thicker relative to the satellite observations.
In the JUL17 experiment, we employ the largest uncertainty of CryoSat-2 ice thickness (1.5 m) as the observational uncertainty input for PDAF. Here we examine the sensitivity of the spatial distribution of assimilated sea ice thickness to (1) different approaches to extrapolate ice thickness based on the calculated PIOMAS seasonal cycle and (2) different observational uncertainties in ice thickness (see Table 2 for details). Figure 9 shows the resulting initial ice thickness for the first ensemble member of the JUL17 experiment. Regardless of different ways to extrapolate ice thickness using the PIOMAS seasonal cycle (Arctic averaged vs. grid based), the use of the maximum uncertainty (1.5 m) of the CryoSat-2 ice thickness

Journal of Advances in Modeling Earth Systems
much thinner ice compared to those employing 1.5-m uncertainty (Figures 9c and 9d). Thus, we apply the reduced uncertainty of the observed ice thickness and the grid-based PIOMAS seasonal cycle for the AUG17, JUL18, and AUG18 predictions (see Table 2). For example, as demonstrated in Figure 4, the improvement is more clearly visible for the initial ice concentration for the AUG17 experiment after the

10.1029/2019MS001938
Journal of Advances in Modeling Earth Systems data assimilation. In particular, the initial sea ice concentration in the area extending from the Beaufort Sea, through the central Arctic, to north of central Siberia is in very good agreement with the observations (Figure 4e vs. Figure 4f), even though the ensemble mean of sea ice extent of AUG17 ( Figure 6 thick red line) has an initial sea ice extent bias of~0.7 × 10 6 km 2 . The spatial distribution of the initial sea ice thickness resembles the CyroSat-2/SMOS/PIOMAS blended ice thickness (Figure 5e vs. Figure 5f). Because of the bias at the initial step, the predicted sea ice extent for August overestimates the observed ice extent nearly 2 weeks and then tracks the observation closely. The delayed freezing also occurs near the end of September, which is similar to that of JUL17. As for the predicted spatial pattern of sea ice concentration of AUG17, there is reduced ice cover in the Chukchi and East Siberian Seas in August and September, which is better than those of JUL17 (Figure 7g).

Journal of Advances in Modeling Earth Systems
Although assimilating satellite sea ice observations improve the accuracy of initial ice concentration and thickness (spatial pattern) for the 2017 predictions substantially, the initial total ice extent shows a positive bias, especially for AUG17 ( Figure 6, red line). This is partly due to the initial ensemble generating procedure, in which the ensemble perturbations are overlaid on the mean state of sea ice trajectories. The mean state of trajectories is the average of 1-month model simulation using the CFSv2 9-month operational

Journal of Advances in Modeling Earth Systems
forecast as initial and lateral boundary conditions (15 June to 15 July for JUL17 and 15 July to 15 August for AUG17). As mentioned above, sea ice simulated by the CFSv2 is too thick over most of the Arctic basin, which leads to slow ice reduction and higher ice extent during the melting season. Assimilating sea ice observations based on initial ensembles with higher mean state results in the aforementioned bias in the initial total ice extent. To address this issue, we modify the ensemble generating procedure by replacing the mean state of sea ice trajectories calculated from 1-month model integration with the CFSv2 initial ice concentration and thickness for the initial date of 2018 predictions (JUL18 and AUG18). As shown in Figure 10, the initial sea ice extent with modified ensemble generating method has comparable bias for JUL18 (very small compared to the observation), but much smaller bias for AUG18, compared to the original method (Figure 10, stars). Figure 10 shows the temporal evolution of the ensemble mean and ensemble spread of the Arctic sea ice extent of JUL18 and AUG18 as well as the CLIM and DAMP predictions. For JUL18, our model captures the observed evolution but slightly underestimates the observed ice extent throughout the prediction period. The observed seasonal minimum extent for 2018 is 4.63 × 10 6 km 2 on 21 September (Figure 10 black line), whereas our predicted minimum is 4.34 (4.78) × 10 6 km 2 on 26 September (21) for JUL18 (AUG18). AUG18 (Figure 10 red line) has an overestimation at the initial (~0.7 × 10 6 km 2 ), although the bias is reduced as discussed above. As a result, the prediction overestimates the actual ice extent for about 1 week. For the rest of the prediction period, the ensemble mean of the ice extent agrees well with the observation. In contrast to 2017, our 2018 model predictions are much better than both the climatology and DAMP predictions. The DAMP predictions of 2018 clearly overestimate (underestimate) the observed sea ice extent for the July (August) prediction and has a larger bias in seasonal minimum of ice extent than our predictions. With the modified ensemble generating procedure, the overestimation of regional ice extent in the East Siberian-Laptev Seas is reduced for JUL18 compared with JUL17 but there is still an overestimation for AUG18.
The modification of the generating procedure for the initial sea ice of 2018 predictions (JUL18 and AUG18) allows us to directly compare our predictions with the CFSv2 predictions. Figure 11 shows the temporal evolution of the predicted sea ice extent of CFSv2 in conjunction with our predictions. Our predictions without data assimilation are also shown for reference (Figure 11 dashed lines). For JUL18, both predictions give reasonable date of the summer minimum around 24 September (Figure 11a). Our summer extent minimum is slightly smaller than the NSIDC observation. However, a large overestimation of 3.83 × 10 6 km 2 (a relative overestimation of 83%) is found for the CFSv2 prediction due to slow ice melting. We also calculate the root-mean-square error (RMSE) of sea ice extent, which is a measure of prediction errors. As shown in Figure 11, both predictions have comparable and small RMSE at the initial as well as the first week of integration. After that, our prediction shows a similar level of RMSE as the first week of integration. However, the RMSE of CFSv2 increases quickly in July, reaching about 3 × 10 6 km 2 in early August, and then gradually increases to about 3.7 × 10 6 km 2 in late September. This is also the case for AUG18.

Journal of Advances in Modeling Earth Systems
The Arctic shipping routes have become increasingly navigable in recent years because of the rapid retreat of sea ice (e.g., Melia et al., 2016;Smith & Stephenson, 2013). The navigability of the ice-covered Arctic Ocean is mainly determined by the real-time ice conditions, including the ice concentration and thickness, and the engineering capabilities of a particular vessel class (Smith & Stephenson, 2013). The temporal evolution of the predicted sea ice extent in the Northern Sea Route (NSR, here the NSR region is defined as the area of 0-180, 65 N-80 N) by our modeling system and CFSv2 is shown in Figure 12. For JUL18, our predicted NSR sea ice extent agrees well the observations and the RMSE stays less thañ 0.2 × 10 6 km 2 during the entire prediction period, whereas the CFSv2 prediction still shows a large overestimation in the NSR region (Figures 12a and 12b). However, for AUG18, our model overestimates the NSR sea ice extent (~0.4 × 10 6 km 2 ), although our prediction still has better skill in the NSR regions than CFSv2 prediction (Figures 12c  and 12d). Although the limitation of the CFSv2 sea ice prediction that is largely due to its assimilation system does not implement sea ice thickness assimilation because of limited availability of real-time thickness observations (Saha et al., 2010), the CFSv2 has irreplaceable value since it provides real-time operational forecasts that can be used as the initial and boundary conditions for regional forecast modeling systems.

Initialization Sensitivity Experiments
As demonstrated above, assimilating satellite-based sea ice concentration and thickness can lead to better estimation of initial ice mass fields and increased predictive skill. The assimilation procedures described in section 3 reduce the biases inherited from CFSv2. The bias correction can be also achieved by the direct insertion of observations. We perform additional numerical experiments (Table 2) using our modeling system including (1) without any assimilation of satellite sea ice data (Exp-NoAss), (2) direct insertion of the SSMIS sea ice concentration (Exp-A), and (3) direct insertion of SSMIS sea ice concentration and CryoSat-2/SMOS/PIOMAS blended sea ice thickness (Exp-B). All missing data in the observed/blended data are replaced with the original CFSv2 data. To maintain the physical consistency between sea ice concentration and thickness, data adjustments were conducted, that is, the ice concentration and thickness in a grid cell are set to 0 if either the ice concentration or ice thickness is 0. Figure 13 shows time series of sea ice

10.1029/2019MS001938
Journal of Advances in Modeling Earth Systems extent from additional numerical experiments along with the JUL18 prediction. The prediction without data assimilation (blue line) significantly overestimates sea ice cover than the observations, although it is less excessive than that of the CFSv2 prediction ( Figure 11). Directly inserting SSMIS sea ice concentration (yellow line) produces adjusted ice extent closer to the observation compared to that of Exp-NoAss, but

Journal of Advances in Modeling Earth Systems
the improvement gradually diminishes after a few weeks' integration, still leading to a significant overestimation of the observed ice extent. This implies that only inserting observed sea ice concentration cannot eliminate the biases in the ice thickness inherited from the CFSv2. By contrast, by direct inserting SSMIS ice concentration and CryoSat-2/SMOS/PIOMAS blended sea ice thickness, the Exp-B (red line)

Journal of Advances in Modeling Earth Systems
shows smaller initial ice extent than those of the Exp-A and the observations and underestimates the observed ice extent throughout the prediction period. Note that its initial underestimation is related to the adjustment that needs to maintain the model physical consistency. Overall, our prediction that assimilates satellite sea ice concentration and thickness and takes account of their covariance and uncertainty leads to improved Arctic sea ice prediction.
As demonstrated by Collow et al. (2015), using PIOMAS ice thickness to initialize NCEP CFS can lead to increased interannual predictive skill. The JUL18_P and the AUG18_P experiments (Table 2) are designed to investigate the difference between assimilating Cyrosat-2 and SMOS ice thickness and assimilating PIOMAS ice thickness in our modeling system. Figure 14 shows the temporal evolution and the ensemble mean and ensemble spread of the Arctic sea ice extent and the RMSE of JUL18_P and AUG18_P in conjunction with JUL18 and AUG18. As shown in Figure 14, the experiments only assimilating PIOMAS data produce larger sea ice extents compared to both the experiments assimilating the CrySat-2/SMOS/PIOMAS blended ice thickness and the observations. Comparing the RMSE of these experiments suggests that all experiments have comparable and small RMSE at the initial as well as during a few weeks' integration. After that, the experiments only assimilating PIOMAS ice thickness have larger RMSE than those assimilating the CrySat-2 and SMOS ice thickness, particularly for the experiment initialized on 1 August (Figure 14b). This suggests that assimilating satellite-based sea ice thickness might be beneficial for Arctic sea ice prediction.

Conclusions and Discussions
In this study, a new regional coupled model that couples the Los Alamos Sea Ice Model (CICE) with the WRF model and the ROMS model under the framework of the COAWST system has been developed for Arctic sea ice prediction. Moreover, we implement the LESTKF in the new coupled model to assimilate SSMIS near-real-time sea ice concentration and CryoSat-2/SMOS merged ice thickness to generate improved model initial conditions. Two forecasting experiments starting from 1 July and 1 August to the end of September are performed for 2017 and 2018, respectively. We examine the impacts of different approaches to extrapolate ice thickness to the starting date of the prediction, different observation errors in ice thickness, and different procedures to generate the ensemble initial state.
The results of 2017 predictions show that the new modeling system configured for the Arctic combined with the assimilation of satellite sea ice data can predict reasonable evolution of the total ice extent compared with the observations, though the predictions exhibit a delayed ice-refreezing timing. Spatially, the JUL17 shows slower ice retreating along the Pacific side of the Arctic Ocean, especially at the Beaufort and Chukchi Seas, which is partly attributed to the initial ice thickness distribution. Sensitivity experiments suggest that the use of appropriate uncertainty in the observed ice thickness can lead to improved spatial distribution of assimilated initial ice thickness. It is demonstrated that our new modeling system initialized by the output of CFSv2 seasonal forecasts with our data assimilation procedures can significantly reduce the sea ice prediction errors compared with the prediction conducted by the CFSv2.
Other important sea ice variables that influence surface energy fluxes and ocean-atmosphere interaction may also need for better initial sea ice conditions, that is, snow depth and melt pond fraction that changes surface albedo property, radiation transfer into the ice layers, and sea ice motion that influences the ice dynamics. The initial states of the atmospheric model and the ocean model also play a crucial role in sea ice predictions. Recent studies showed that assimilating additional atmospheric observations (i.e., radiosonde) have significant influences on predicting atmospheric variation and thus lead to improved weather and sea ice prediction (Inoue et al., 2013(Inoue et al., , 2015. Bushuk et al. (2017) showed that the initialization of the ocean subsurface provides significant skill for regional winter sea ice extent at a seasonal timescale. Further experiments will be performed to understand the impacts of data assimilation of other variables with our new modeling system.
In the fully coupled predictive model, sea ice strongly interacts with the atmosphere above and the ocean below. Thus, biases in the atmosphere and the ocean may potentially degrade the predictive capability for sea ice. Here we conduct a preliminary evaluation of the simulation of the atmosphere and ocean, focusing on near-surface air temperature and SST, which play important roles in determining the melting of sea ice during the prediction period. Figure 15 shows the spatial distribution of the ERA-interim 2-m air temperature (T2) and predicted anomalies (ensemble mean minuses ERA-interim) for July, August, and September in 2017. For JUL17, compared to the ERA-interim, the predicted temperature in July has small cold (warm) biases over the Arctic Ocean and northern Pacific and Atlantic (Canada) (Figures 15a and 15d). Moderate-to-large cold biases (~4-6°) are found over central and western Siberia. For August (Figure 15e), the cold biases over the northern Pacific and Atlantic are increased relatively, and large cold biases (~6-10°) covers almost entire Siberia and eastern Europe. In contrast to JUL17, AUG17 shows a similar bias pattern with relatively enhanced magnitude for cold biases (Figure 15g). However, the cold biases over the land surface become smaller in September (Figure 15h). Although the significant test for near-air surface temperature suggests that our predictions show a good agreement with the observations over large part of land surface, efforts are needed to reduce the large cold biases over the land surface, particularly predictions with longer lead times. Recently, Kuo et al. (2018) suggested that surface emissivity is important for reducing high-latitude temperature biases in the Community Earth System Model. Figure 16 shows the spatial distribution of the NOAA optimum interpolation sea surface temperature and predicted anomalies (ensemble mean minuses NOAA) for July, August, and September in 2017. In general, both JUL17 and AUG17 show cold biases in the North Pacific and Atlantic. For JUL17, the cold bias is increased from July to August compared to optimum interpolation sea surface temperature, especially in the large portion of the North Pacific, Greenland Sea, Labrador Sea, and Sea of Okhotsk in August (~3-5°, Figures 16d-16f). Such bias pattern maintains in September. In contrast to JUL17, AUG17 shows smaller cold biases over the ocean, but the cold biases still enlarge relatively over time, especially in the Greenland Sea and Labrador Sea (Figures 16g and 16h). Zhu et al. (2015) discussed reducing SST bias in CFSv2 prediction with the ocean initialization scheme. They suggested that using a simple SST-nudging scheme is an effective way to enhance predictive skill on SST for most of ocean basins. The significant test suggests that sea ice edge zones have good consistency between the predicted values and the observed values.

Journal of Advances in Modeling Earth Systems
We also find that the cold biases in both predicted T2 and SST in our new modeling system are partly due to the underestimation of the downward shortwave radiation at the surface compared to the ERA-interim (not shown). This underestimation of downward shortwave radiation is associated with the simulated cloudiness by the Grell-Freitas cumulus scheme and Morrison two-moment microphysics, which tend to produce too much cloud that greatly reduces the downward shortwave radiative flux. Further investigation of simulated radiation for our modeling system will be conducted in future research.
Although this new modeling system has shown great potential for seasonal Arctic sea ice prediction, more experiments with this modeling system for seasonal Arctic sea ice and climate prediction will be conducted to examine the robustness of the modeling system and reduce the biases in three components of the model. For instance, physical processes which are missing or poorly represented in the new modeling system need to be improved, and a more sophisticated data assimilation involving more state variables requires further investigation.