Daily Scale Prediction of Arctic Sea Ice Concentration Based on Recurrent Neural Network Models

1. Introduction

Global warming has accelerated the rate of melting of the Arctic sea ice [1]. During 1970-2010, Arctic sea ice area decreased by an average of 4% per decade [2]. The rate of sea ice area decline increased dramatically into the 21st century [3,4,5]. The extent of the Arctic sea ice reached its smallest value in the recorded history of satellite data in 2012, at about 3.34 million ${\mathrm{km}}^2$, while the second lowest value occurred in 2020, at about 3.74 million ${\mathrm{km}}^2$[6]. The reduced extent of the sea ice presents new opportunities for a number of industries, including Arctic shipping, tourism, fisheries, and oil and gas exploration [7]. Predicting seasonal and daily scale changes in Arctic sea ice is of great practical significance for the safe operation of Arctic shipping lanes and the development and utilization of Arctic resources.

Physical interactions between the atmosphere, ocean, and sea ice are the basis for predicting sea ice. Several studies have already been conducted to predict Arctic sea ice on different spatial and temporal scales and to explore the predictability in different seasons. Blanchart et al. [8] used the outputs of the Community Climate System Model, version 3 (CCSM3) to assess the mechanisms of sea ice persistence, and comparisons with actual observations demonstrated that the model can be used for seasonal to annual predictions of sea ice. Krikken et al. [9] used 15 climate models of CMIP5 to analyze the natural variability in Arctic sea ice from an energy balance perspective and found a strong correlation between the energy balance and the reappearance of sea ice anomalies from the sea ice melting season to the growing season. Guemas et al. [10] reviewed potential sources of Arctic sea ice predictability from months to years, including the persistence and advection of sea ice anomalies, interactions with the oceanic atmosphere, and changes in the radiative forcing. Mohammadi et al. [11] determined the potential predictability of Arctic winter sea ice using a sea ice–ocean coupling model, noting high predictability of sea ice concentration (SIC) and sea ice edge position over a 10-day period. Cruz et al.[12] investigated sea-ice predictability from seasonal to annual scales using a variety of climate models, emphasizing the importance of the reoccurring effects of sea ice anomalies, and observed that anomalies in SIC in the Barents Sea are highly negatively correlated with local sea surface temperature anomalies. Onarheim et al.[13] showed that ocean heat transport (OHT) variability plays an important role in winter sea ice variability in the Barents Sea, and that the use of the OHT can lead to predictions two years in advance. These studies provide knowledge on the predictability of sea ice associated with a wide range of physical processes, while emphasizing the importance of selecting predictors that are relevant to the target location and time scale.

The main approaches to Arctic sea ice prediction are numerical simulations, statistical predictions and deep learning. Numerical simulation methods are based on physical links between temperature changes, humidity transport, wind field models, cloud cover and ocean heat fluxes. Major climate simulation centers around the world have released some atmospheric and oceanic simulation data. Their adoption of climate models relies on real-time inputs of observational conditions in the data assimilation process. At the same time, due to the limitations of a single climate model and the large differences in the results of different models, it is necessary to determine the assigned weights of each climate model based on its contribution to the simulation of the current climate mean, and then take a weighted average to improve the accuracy of sea ice prediction. In practice, this treatment does not eliminate the effect of model bias on sea ice prediction. In addition, many processes in the dynamic model of sea ice need parameterization, and the current model lacks modeling of rheology, ice thickness distribution, wave–ice interaction, landing ice, melting water and size distribution of floating ice[14]. Statistical methods are used to predict the state of sea ice according to historical data. With the thinning of sea ice, the average state of sea ice has changed significantly. In addition, most statistical models are linear models, which cannot learn the nonlinear relationship between variables in the Arctic climate system. Because the nonlinear feedback mechanism plays an important role in the coupling system of Arctic atmosphere, ocean and sea ice, it is necessary to predict Arctic sea ice using a nonlinear model. Deep learning technology has a strong nonlinear learning ability. Chi et al. and Choi et al. input the monthly average SIC data of the National Snow and Ice Data Center (NSIDC) into the multilayer perceptron (MLP) and the long and short-term memory (LSTM) models to predict the monthly average of SIC, and found that the results are better than the traditional autoregressive (AR) model[15,16]. Kim et al.[17] used the integrated data of a regional climate model (RCM) as input variables, and used the deep neural network (DNN) method to deal with the nonlinear relationship between SIC and climate variables. As a result, they predicted the SIC in the Kara Sea and Barents Sea in the next 10-20 years. Fritzner et al.[18] compared the prediction accuracy of a high-resolution dynamic assimilation model, K-NN model and FCN model for the next 7 days, and pointed out that the FCN model can provide similar prediction results to the dynamic assimilation model. Kim et al.[19] input eight predictors into the CNN model to predict the monthly average SIC in the next month, and the results were better than that of the RF model. Andersson et al.[20] put forward a model of a probabilistic sea ice prediction system. The model used climate simulation and observation data as the input data to predict the monthly average SIC in the next six months. The results showed that the IceNet model has a high accuracy in predicting the sea ice range, and it is superior to the SEAS5 dynamic model in predicting extreme sea ice events in summer. Liu et al.[21] used the convolutional LSTM (ConvLSTM) model to predict SIC in the Barents Sea over the next six weeks. They added ERA-Interim reanalysis data to the training dataset, and used the covariance between different variables and the spatiotemporal correlation to complete the prediction of regional SIC. The results were better than the linear regression model. At present, the deep learning method is mainly used to predict the sub-seasonal scale of regional sea ice. Because the daily short-term forecast of SIC is very important for maritime shipping decision-making, it is urgent to attain the accurate daily short-term forecast of SIC.

The ConvLSTM and predictive recurrent neural network (PredRNN) models exhibit the capacity to capture spatiotemporal correlations among diverse input parameters, enabling them to theoretically predict spatiotemporal sequence data. This study introduces these models into the realm of the high-precision daily scale short-term prediction of Arctic sea ice. Initially, we compare the predictive performance of the ConvLSTM and PredRNN models when only SIC is utilized as the input. Subsequently, we enhance the input data by incorporating meteorological parameters that influence both SIC and the sea boundary, leading to the formation of ConvLSTM-multi and PredRNN-multi models. Through an investigation of the spatiotemporal correlations between SIC and meteorological parameters, we observe a substantial enhancement in the model's predictive capability for the sea ice edge region. Furthermore, this paper conducts a quantitative analysis to discern the model's sensitivity to the input meteorological parameters, pinpointing the key meteorological variables that affect the prediction accuracy of SIC.

2. Data and methods

2.1. Data

We used SIC data from the NSIDC and reanalysis data from ERA5, provided by the European Centre for Medium-Range Weather Forecasts (ECMWF), as the training set, spanning the period from 1988 to 2021. The NSIDC SIC dataset is derived from observations made by the scanning multichannel microwave radiometer (SMMR) carried by the Nimbus-7 satellite, the special sensor microwave imager (SSM/I) sensors carried by the National Defense Meteorological Satellite Program (DMSP)-F8, -F11 and -F13 satellites and the special sensor microwave imager/sounder (SSMIS) sensors carried by DMSP-F17. Among the NSIDC datasets, we utilized the NASA Bootstrap algorithm-derived dataset for our research. The product is provided in the form of a daily average, with polar stereo projection (45°W, 70°N), spatial resolution of 25km×25km and grid number of 448×304. The values within this dataset range from 0, signifying the absence of ice in the grid, to 100, indicating complete sea ice coverage. The primary sources of error in these data arise from thin ice (ranging from 30% to 50%) and surface melting (ranging from 10% to 30%) [22]. The Arctic Ocean covered by this dataset predominantly encompasses the Sea of Okholtsk, Bering Sea, Chukehi Sea, Beaufort Sea, Canadian Archipelago, Hudson Bay, Baffin Bay, Greenland Sea, Norwegian Sea, Barents Sea, Kara Sea, Laptev Sea, East Siberian Sea and Central Arctic Sea. The distribution of sea areas is illustrated in Figure 1. The NSIDC takes SIC=15% as the sea ice boundary threshold, and the region with SIC≥15% is considered as the sea ice area.

The ERA5 dataset is a global meteorological dataset in which the ECMWF combines meteorological model data with observation data from all over the world, with a spatial resolution of $0.25°\times 0.25°$ [23]. Table 1 details the specific parameters within the ERA5 dataset utilized as inputs for our model, include sea surface temperature (SST), 2m temperature (T2M), skin temperature (SKT), surface solar radiation downwards (SSRD), mean sea level pressure (MSL), 10m u-component of wind (U10) and 10m v-component of wind (V10). The hourly value of ERA5 data is converted into a daily average, resampled into a grid consistent with SIC data, and the value is normalized to [0,1].

SIC, ERA5 15 days in advance, land mask, cosine and sine values of dates are fused together to form a series of 10 consecutive days as training samples. The dataset is divided into training dataset (1988-2018), verification dataset (2019) and test dataset (2020-2021). The Arctic sea ice area in September 2020 was the second lowest level in the recorded history of satellite data (greater only than in September 2012)[24], which serves as a fitting evaluation point for assessing the model's predictive performance under extreme conditions. All three datasets were partitioned into 10-day sequences to achieve a random input of 10 days of historical data to predict the SIC for the next 10 days. To analyze to what extent and how the atmospheric conditions and oceanic variables affected the accuracy of the model predictions, the models were trained using the mode of inputting both single and multiple predictors (Table 1).

In addition, the sixth phase of the Coupled Model Intercomparison Project (CMIP6) "selected models" nominated by the Sea-Ice Model Intercomparison Project (SIMIP) community[25] were used. SIMIP aims to compare and evaluate the performance of different sea-ice models, collect different simulation results and establish a standardized database, identify model strengths and weaknesses and uncertainties, improve the accuracy of sea ice models, and lay the foundation for future sea-ice prediction studies. The “selected models” in Table A 1 provide the best estimates of the future evolution of Arctic sea ice, and the daily average SIC provided by three CO2 emission scenarios, SSP126, SSP245, and SSP585, were chosen for comparison with the deep learning model.

2.2. Model

SIC data should be predicted using the spatiotemporal sequence prediction model. The spatiotemporal sequence is a dynamic system in which historical observations with arbitrary length J evolve over time, and observations at each moment can be represented on an M×N grid. Thus, the observation at any time can be represented by a tensor $\mathcal{X}\in {ℝ}^{J\times M\times N}$ , where J denotes the domain of the observed features. If we record observations periodically, we will obtain a sequence of tensors, while observations over period of T are denoted as ${\mathcal{X}}_{in}=\left\{{\mathcal{X}}_1,\dots ,{\mathcal{X}}_T\right\}$. The model is designed to predict the sequence ${\widehat{\mathcal{X}}}_{\mathrm{out}}=\left\{{\widehat{\mathcal{X}}}_{T+1},\dots ,{\widehat{\mathcal{X}}}_{T+K}\right\}$ for the next K time steps, given ${\mathcal{X}}_{in}$. For the training pairs ${\left\{\left({\mathcal{X}}_{\mathrm{in}}^n,{\mathcal{X}}_{\mathrm{out}}^n\right)\right\}}_n$ formed by all SIC data, a set of parameters ${\theta }^{\prime }$ is found by using random gradient descent to ensure that the log-likelihood of the generated target sequence ${\mathcal{X}}_{out}$ is the maximum when the input data ${\mathcal{X}}_{in}$ are provided:

$$\theta \text{'}=\underset{\theta }{\arg \max }\sum _{\left({\mathcal{X}}_{\mathit{in}}^n,{\mathcal{X}}_{\mathit{out}}^n\right)}\log P\left({\mathcal{X}}_{\mathit{out}}^n\mid {\mathcal{X}}_{\mathit{in}}^n;\theta \right)$$

(1)

In this study, we present the ConvLSTM model [26] and PredRNN model [27], which are commonly employed for spatiotemporal sequence data prediction, to realize the prediction of SIC. The ConvLSTM model replaces matrix multiplication with convolution operations within the LSTM gating structure unit. This alteration allows it to simultaneously capture both temporal and spatial features in the data, effectively transforming the traditional encoding–decoding structure into an encoding–forecasting structure. The gating formula for the ConvLSTM unit is:

$$\begin{array}{c}{g}_t=\tanh \left(W_{xg}\ast {\mathcal{X}}_t+W_{hg}\ast {\mathcal{H}}_{t-1}+b_g\right)\\ i_t=\sigma \left(W_{xi}\ast {\mathcal{X}}_t+W_{hi}\ast {\mathcal{H}}_{t-1}+W_{ci}\odot {\mathcal{C}}_{t-1}+b_i\right)\\ f_t=\sigma \left(W_{xf}\ast {\mathcal{X}}_t+W_{hf}\ast {\mathcal{H}}_{t-1}+W_{cf}\odot {\mathcal{C}}_{t-1}+b_f\right)\\ {\mathcal{C}}_t=f_t\odot {\mathcal{C}}_{t-1}+i_t\odot g_t\\ o_t=\sigma \left(W_{xo}\ast {\mathcal{X}}_t+W_{ho}\ast {\mathcal{H}}_{t-1}+W_{co}\odot {\mathcal{C}}_t+b_o\right)\\ {\mathcal{H}}_t=o_t\odot \tanh \left({\mathcal{C}}_t\right)\end{array}$$

(2)

where $i_t$ is the input gate; $f_t$ is the forgetting gate; ${\complement }_t$ is the unit storage state; $o_t$ is the output gate; ${\mathcal{H}}_t$ is the hidden state; $\mathrm{W}$ is the weight matrix, with the subscript describing the correspondence between the weight matrix and the state of each gate; $\mathsf{\chi }$ is the input; $\mathrm{b}$ is the bias; ∗ is the convolution operation; $\odot$ is the Hadamard product; $\mathsf{\sigma }$ is the sigmoid activation function; $\tan \mathrm{h}$ is the hyperbolic tangent activation function; and the subscript $\mathrm{t}$ denotes the time step. Figure 2 illustrates the ConvLSTM network framework, which consists of a three-layer stack of ConvLSTM units employed in the encoding–forecasting framework. In this framework, the encoder incorporates a down-sampling operation, while the forecaster utilizes an up-sampling operation before each layer of input at each time step. The training data input follows the traditional sequence-to-sequence approach [28]. Historical observations are input into the encoder during the training phase, and the state layers generated at each layer of the encoder (the shaded region in Figure 2) are transmitted to the forecaster. The loss function is determined by comparing predicted values to actual values, and model parameters are adjusted through backpropagation until Equation 1 is satisfied.

The PredRNN exhibits three distinctions from the ConvLSTM framework: (1) It introduces a spatiotemporal memory-state circulation method, depicted by the blue arrows in Figure 3, which enhances the lower layer's capability to learn top-layer features from the previous time step; (2) The unit responsible for the spatiotemporal memory stream employs a dual-stream memory transition mechanism involving ${\complement }_t$and ${\mathsf{{\rm M}}}_t$. This results in memory states that cannot be decoupled spontaneously. To address this, a convolutional layer is added to the ${\complement }_t$and ${\mathsf{{\rm M}}}_t$ increments at each time step, and the spatial distances between them are extended using a novel decoupling loss function. This approach trains different memory states to focus on long-term and short-term spatiotemporal changes. (3) The model training approach incorporates the reverse scheduled sampling (RSS) [29] method for data input. This technique forces the model to randomly conceal true observations in the encoder to learn long-term dynamics, with the probability of concealing true observations decreasing with the number of iterations. This ensures that the model has the same likelihood of inputting true observations during both the training and prediction phases. In (1), the network is enabled to learn the complex nonlinear variations in short-term motions. However, the state layer transfer path stretching across time brings the problem of gradient vanishing, making it challenging to capture long-term dependencies. Therefore, a dual-stream memory transition mechanism is introduced to achieve a short-term recursion depth and long-term consistency by combining the original memory unit ${\complement }_t$ and the new memory unit ${\mathsf{{\rm M}}}_t$ to form a unit named ST-LSTM, which is calculated as follows:

$$\begin{array}{l}{g}_t=tanh(W_{xg}\ast {\chi }_t+W_{hg}\ast {\mathcal{H}}_{t-1}^l+b_g)\\ i_t=\sigma (W_{xi}\ast {\chi }_t+W_{hi}\ast {\mathcal{H}}_{t-1}^l+b_i)\\ f_t=\sigma (W_{xf}\ast {\chi }_t+W_{hf}\ast {\mathcal{H}}_{t-1}^l+b_f)\\ C_t^l=f_t\odot C_{t-1}^l+i_t\odot g_t\\ g_t^{\text{'}}=tanh(W_{xg}^{\text{'}}\ast {\chi }_t+W_{mg}\ast M_t^{l-1}+b_g^{\text{'}})\\ i_t^{\text{'}}=\sigma (W_{xi}^{\text{'}}\ast {\chi }_t+W_{mi}\ast M_t^{l-1}+b_i^{\text{'}})\\ f_t^{\text{'}}=\sigma (W_{xf}^{\text{'}}\ast {\chi }_t+W_{mf}\ast M_t^{l-1}+b_t^{\text{'}})\\ M_t^l=f_t^{\text{'}}\odot M_t^{l-1}+i_t^{\text{'}}\odot g_t^{\text{'}}\\ o_t=\sigma (W_{xo}\ast {\chi }_t+W_{ho}\ast {\mathcal{H}}_{t-1}^l+W_{co}\ast C_t^l+W_{mo}\ast M_t^l+b_o)\\ H_t^l=o_t\odot tanh(W_{1\times 1}\ast [C_t^l,M_t^l])\end{array}$$

(3)

The PredRNN network framework comprises four layers of interconnected ST-LSTM units. As illustrated in Figure 3, the input data undergo reverse scheduled sampling (RSS), progressively increasing the likelihood of true observations being incorporated in the encoder, while inversely decreasing in the forecaster. The state unit ${\mathcal{H}}_t^l$ and memory unit $M_t^l$ circulate along the orange arrows in the diagram. Subsequently, the memory unit $M_t^l$ flows directly along the blue arrows, transitioning from the uppermost layer of the preceding time step to the lower layer of the subsequent time step, thereby establishing the circulation of short-term memory states.

2.3. Evaluation metrics

Mean absolute error (MAE) (Equation 4), root mean square error (RMSE) (Equation 5), normalized root mean square error (nRMSE) (Equation 6), anomalous correlation coefficient (ACC) (Equation 7), Nash–Sutcliffe efficiency (NSE) (Equation 8), and structure similarity index measure (SSIM) (Equation 9) are used to assess the model's performance in predicting SIC. MAE and RMSE are used to measure the absolute difference between the model's predictions and observed values, with lower values indicating superior predictive capability. nRMSE is computed by dividing RMSE by the standard deviation of observed values, thereby offering a more precise depiction of the residuals in the sea-ice edge region within the model's predictions. It is expressed as a percentage, with lower values signifying reduced residuals. ACC serves as an indicator of the fidelity of predicted anomalies and the degree to which predicted values align with the actual data, producing values within the range of +1 to -1. A value closer to +1 denotes greater consistency between predicted and observed values. NSE is employed to gauge the precision of the model's output values within the range of $-\infty$ to 1. Values approaching 1 signify a more accurate model. SSIM quantifies the structural resemblance between predicted and observed values, with values approaching 1 indicating a higher structural similarity.

$$MAE=\frac{\left|X_{obs,i}-X_{model,i}\right|}{n}$$

(4)

$$RMSE=\sqrt{\frac{{\displaystyle {\sum }_{i=1}^n{(X_{obs\text{},i}-X_{model,i})}^2}}{n}}$$

(5)

$$nRMSE=\sqrt{\frac{{\displaystyle {\sum }_{i=1}^n{(X_{obs\text{},i}-X_{model,i})}^2}}{{\displaystyle {\sum }_{i=1}^n{(X_{obs\text{},i}-\overline{X_{obs}})}^2}}}$$

(6)

$$ACC=\frac{{\displaystyle {\sum }_{i=1}^n(X_{model,i}-\overline{X_{model}})(X_{obs,i}-\overline{X_{obs}})}}{\sqrt{{{\displaystyle {\sum }_{i=1}^n(X_{model,i}-\overline{X_{model}})}}^2{{\displaystyle {\sum }_{i=1}^n(X_{obs,i}-\overline{X_{obs}})}}^2}}$$

(7)

$$NSE=1-\frac{{\displaystyle {\sum }_{i=1}^n{(X_{obs,i}-X_{model,i})}^2}}{{\displaystyle {\sum }_{i=1}^n{(X_{obs,i}-\overline{X_{obs}})}^2}}$$

(8)

$$SSIM\left(obs,model\right)=\frac{\left(2{\mu }_{obs}{\mu }_{model}+C_1\right)\left(2{\sigma }_{obsmodel}+C_2\right)}{\left({\mu }_{obs}^2+{\mu }_{model}^2+C_1\right)\left({\sigma }_{obs}^2+{\sigma }_{model}^2+C_2\right)}$$

(9)

Model prediction errors primarily manifest at the sea ice edge, an area of significance for Arctic shipping and navigation. Therefore, it is imperative to validate the prediction accuracy of the sea-ice edge position. Melsom et al. [30] proposed employing three metrics, namely, the mean ice edge displacement (${\mathrm{D}}_{\mathrm{AVG}}^{\mathrm{IE}}$), the integrated ice edge error (IIEE) average displacement ($D_{\mathrm{AVG}}^{\mathrm{IIEE}}$), and the IIEE bias (${\mathsf{\Delta }}^{\mathrm{IIEE}}$), to assess the accuracy of the model's sea-ice edge prediction. ${\mathrm{D}}_{\mathrm{AVG}}^{\mathrm{IE}}$ measures the shortest Euclidean distance between the actual observed sea-ice edge grid cell points and the model's sea-ice edge grid cell points (Equation 10). Here, $N_o$, $N_m$ denote the number of grids for observed and predicted values, while $d_o^{\mathrm{n}}$, $d_m^{\mathrm{n}}$ represent the distance displacements between the observed and predicted values corresponding to the nth grid of the edge. $D_{\mathrm{AVG}}^{\mathrm{IIEE}}$ defines the integral displacement between the observed and predicted ice edges (Equations 11-12). Error estimation by integration reduces the effect of small-sized localized ice features (e.g., polygonal openings) on the total displacement[31]. ${\mathsf{\Delta }}^{\mathrm{IIEE}}$ quantifies the disparity between observed and predicted ice (Equation 13), with a positive deviation indicating that the predicted ice exceeds the observed value, and vice versa.

In these equations, $c_m,c_o$ and $c_e$ represent the predicted values at grid points, the observed values at grid points, and the constant values defining the sea ice boundaries, respectively. $L_O$, $L_M$are the observed and predicted ice edge lengths. $A^{IIEE}$ encompasses the total number of grid cells in the over-predicted and under-predicted ice regions, with ${\alpha }^{IIEE}$ representing the difference between them. ${\gamma }_{AVG}$ determines the robustness of sea-ice edge error measurement results, with larger values indicating greater sensitivity to the formulation of the sea-ice edge displacement error (Equation 14).

$$D_{\mathit{AVG}}^{\mathit{IE}}=\frac{1}{2}\left[\frac{1}{N_O}{\displaystyle \sum }_{n=1}^{N_O}{d}_o^n+\frac{1}{N_M}{\displaystyle \sum }_{n=1}^{N_M}{d}_m^n\right]$$

(10)

$$\begin{array}{l}{A}^{\mathit{IIEE}}=A^{+}+A^{-},A^{+}={\displaystyle \sum _A{a}^{+}},a^{+}=\left\{\begin{array}{l}a,c_m>c_{\mathrm{e}}\wedge c_o<c_{\mathrm{e}}\\ 0\end{array}\right.\\ {\alpha }^{\mathit{IIEE}}=A^{+}-A^{-},A^{-}={\displaystyle \sum _A{a}^{-}},a^{-}=\left\{\begin{array}{l}a,c_o>c_{\mathrm{e}}\wedge c_m<c_{\mathrm{e}}\\ 0\end{array}\right.\end{array}$$

(11)

$$D_{\mathit{AVG}}^{\mathit{IIEE}}=\frac{2}{L_O+L_M}{A}^{\mathit{IIEE}}$$

(12)

$${\mathsf{\Delta }}^{\mathit{IIEE}}=\frac{2}{L_O+L_M}{\alpha }^{\mathit{IIEE}}$$

(13)

$$r_{\mathit{AVG}}=\frac{D_{\mathit{AVG}}^{\mathit{IE}}}{D_{\mathit{AVG}}^{\mathit{IIEE}}}$$

(14)

3. Results and discussion

3.1. Daily scale predictions of sea ice concentration

Figure 4 depicts the distribution of average prediction accuracy metrics (MAE, RMSE, nRMSE, ACC, NSE, and SSIM) values for the ConvLSTM, PredRNN, ConvLSTM-multi, and PredRNN-multi models for the years 2020 and 2021. As shown in Figures 4a and 4b, the daily average prediction accuracy of SIC in 2020 and 2021 follows the order of ConvLSTM-multi, PredRNN-multi, ConvLSTM, and PredRNN, from high to low, for MAE and RMSE. The metrics nRMSE, ACC, and NSE exhibit sensitivity to the prediction accuracy of the sea ice edge region. As evidenced by Figure 4c, d, and e, the prediction accuracy of the multi-predictor models ConvLSTM-multi and PredRNN-multi in the thin ice region significantly outperforms that of the single predictor models. PredRNN-multi demonstrates superior prediction accuracy in the initial six days, but it is subsequently surpassed by ConvLSTM-multi. Figure 4f illustrates that the ConvLSTM model outperforms the PredRNN model in terms of the SSIM between predicted values and observed values. Notably, the ConvLSTM-multi model excels, capturing the shape changes in Arctic sea ice coverage for 2020 and 2021. Figure A 2 shows the SIC distribution of observed and predicted values from September 6 to September 15, 2020, for ConvLSTM, PredRNN, ConvLSTM-multi, and PredRNN-multi models. It is evident that the sea ice edge's shape in the ConvLSTM-multi model aligns most closely with the observed values.

Overall, the prediction accuracy of the four models significantly surpasses the average prediction accuracy of the CMIP6 model in three prospective climate scenarios (Table 2). On the 10th day, the prediction accuracy of the four models reaches its weakest values (MAE: 8.45%, RMSE: 15.96%, nRMSE: 48%, ACC: 0.88, NSE: 0.76), which is superior to the ensemble prediction accuracy of the best model (SSP126 in 2020: MAE: 19.67%, RMSE: 29.13%, nRMSE: 69%, ACC: 0.76, NSE: 0.53).

Figure 5 and Figure 6 display comparative histograms illustrating the distribution of SIC intervals, observed from NSIDC data, and those predicted by four distinct models during the melting season spanning from June to September in the years 2020 and 2021. Since NSIDC assigns all regions with less than 15% SIC to 0, the number of grids with 0 in the predicted SIC values of the four models is obviously less than that of NSIDC SIC. The comparison results show that the number of grids predicted by ConvLSTM and ConvLSTM-multi models is high in 40%~50% and low in 50%~80%. The predicted values of PredRNN and PredRNN-multi models exhibit a remarkable level of agreement with NSIDC SIC in 30%~90%; The predicted values of the four models are low in 90%~95% and high in 95%~100%. The difference between the predicted grid number and the observed grid number of the PredRNN-multi model in 80%~100% is smaller than that of the PredRNN model. The interval distribution of SIC predicted by the PredRNN-multi model is most consistent with the distribution of NSIDC SIC in the melting season.

The RMSE of the daily scale SIC predicted values of the four models in 2020–2021 changes with the prediction time, as shown in Figure 7. The prediction error is concentrated within the one-year ice region, with minimal error observed in the multi-year ice region. Specifically, RMSE values initially rise in the Kara Sea and Barents Sea, followed by a gradual increase in the Laptev Sea and East Siberian Sea. The ConvLSTM-multi model demonstrates superior performance, displaying the slowest increase in RMSE with respect to the number of prediction days, particularly in the Kara Sea and Barents Sea. The inflow of warm salt water from the Atlantic Ocean[32,33] and positive solar radiation accumulated in summer[34] caused the sea ice in Kara Sea and Barents Sea to shrink, and the SIC to change rapidly, so the model prediction error of these two sea areas was large. There is annual sea ice drift from the Laptev Sea and East Siberian Sea to Fram Strait[35,36], which causes the sea ice along the line to change rapidly, and correspondingly increases the prediction error of the model along the line.

The ACC spatial distribution of the daily scale SIC predicted values of the four models in 2020 -2021 changes with the prediction time, as shown in Figure 8. The ACC of the one-year ice region is high, which shows that the anomalies predicted by the model are in good agreement with those observed in practice, and the model demonstrates its efficacy in accurately capturing the dynamic transition from sea ice melting to regrowth within the one-year ice region. SIC in the perennial ice area remains stable throughout the year, and the prediction error is small, so the ACC value is low. With the increase in forecast days, the consistency between the forecast anomaly and the observation anomaly in a year's ice area decreases. The ConvLSTM-multi model exhibits the highest performance, and the ACC of its prediction results decreases the slowest with the number of prediction days.

Figure A1 shows the RMSE and ACC of the "selected models" in CMIP6 for three future CO2 emission scenarios (SSP126, SSP245, and SSP585), with the RMSE of the predicted values in the SSP126 scenario being the best, and the RMSE of the predicted values for the one-year ice ranging from 20% to 50%, while the RMSE of the predicted values of the four models in the 10th day of the one-year ice ranges from 12.5% to 25%, indicating that the deep learning models can better predict the process of sea ice from melting to regrowth in the one-year ice. The ACC of the predicted values in the SSP126 scenario is optimal and is mainly distributed between 0.7 and 0.8 in the one-year ice, whereas the ACC of the predicted values of the four models on day 10 in the one-year ice region is mainly distributed between 0.9 and 1, which indicates that the deep learning models are more capable of predicting the anomalies that need to be captured.

3.2. Sea ice edge prediction accuracy

Figure 9 shows how the IIEE values of the four models change with the prediction time, and the IIEE values of the ConvLSTM-multi model increase the slowest with the prediction time. By contrasting the predictive outcomes of the ConvLSTM model with those of the ConvLSTM model supplemented with reanalysis data in the training dataset, it becomes evident that the inclusion of reanalysis data effectively enhances predictive accuracy for the sea ice edge region and mitigates the decline in accuracy as the prediction time lengthens. Similarly, the incorporation of reanalysis data in the PredRNN model also bolsters predictive accuracy for the sea ice edge region. In the extreme year (2020), the PredRNN-multi model yields higher IIEE values in predictions compared to the ConvLSTM model, while in the standard year (2021), it delivers lower values. This implies that the PredRNN-multi model falls short of matching the predictive capabilities of the ConvLSTM-multi model in thin ice region forecasting. Figure A3 provides a visual representation of the spatial distribution of IIEE for the four SIC prediction models for the period September 6-15, 2020, using SIC=15% as the threshold for distinguishing over-predicted ice regions ($A^{+}$) and under-predicted ice regions ($A^{-}$). The visual analysis corroborates that the ConvLSTM-multi model consistently demonstrates superior performance by exhibiting the least pronounced increase in IIEE values as the prediction time extends.

Figure 10 shows the distribution of ${\mathrm{D}}_{\mathrm{AVG}}^{\mathrm{IE}}$,$D_{\mathrm{AVG}}^{\mathrm{IIEE}}$ and ${\mathsf{\Delta }}^{\mathrm{IIEE}}$ across four model predictions during the winter–spring (freezing) and summer–autumn (melting) periods. ${\mathrm{D}}_{\mathrm{AVG}}^{\mathrm{IE}}$ is the point-to-point displacement between the predicted and observed ice edge, representing the upper limit of the sea-ice edge position displacements. $D_{\mathrm{AVG}}^{\mathrm{IIEE}}$ represents the integral area displacement between the predicted and observed ice edge, delineating the lower limit of the sea-ice edge position displacement. ${\mathsf{\Delta }}^{\mathrm{IIEE}}$ quantifies the deviation in the predicted total ice content from the observed total ice content, and ${\gamma }_{AVG}$ reflects the sensitivity of the actual error in the predicted ice-edge position concerning the applied displacement error measure. In the comparison of results from Figure 10 a, c, and e, it is evident that the deviation between the sea ice edge predicted by the four models during winter–spring and the observed values is minimal, exhibiting a gradual increase with the extension of the prediction time. Conversely, when contrasted with the outcomes from Figure 10b, d, and f, it becomes apparent that the deviation in the sea ice edge predicted by the four models during summer–autumn increases approximately three times faster than during winter–spring, and the trajectory of each index displays marked divergence. This discrepancy is primarily attributed to the more substantial changes in the sea ice edge area during summer–autumn, rendering predictions inherently more challenging. Analysis of the results from Figure 10 f reveals that the ConvLSTM model demonstrates a greater ability to anticipate the reduction in total sea ice content during the melting season compared to the PredRNN model. However, it exhibits inadequate learning capacity during the early stages of sea ice freezing (October–November). Referring to the sea ice edge displacement results presented in Table 3, the ConvLSTM-multi model exhibits superior performance during summer–autumn, while the PredRNN-multi model excels during spring–winter. Moreover, the mean values of ${\mathrm{D}}_{\mathrm{AVG}}^{\mathrm{IE}}$, $D_{\mathrm{AVG}}^{\mathrm{IIEE}}$, and $\left|{\mathsf{\Delta }}^{\mathrm{IIEE}}\right|$ for the four models during both winter–spring and summer–autumn in 2021 surpass those of the CMIP6 models within the "selected models". The accuracy of the four models in predicting the sea ice edge's location is notably superior during winter–spring when compared to summer–autumn, as well as in contrast to the CMIP6 model.

3.3. Parameter Sensitivity Analysis

Examining the sensitivity of model prediction results to input parameters contributes to enhancing the selection of model input variables and improving prediction accuracy and efficiency. To assess the influence of specific parameters, the target parameters and SIC are retained, while all other parameters are replaced with noise. The RMSE is then for the model's predicted values. Comparative analysis of the RMSE values generated by different parameter types offers insights into the significance of each parameter's impact on the model's predictions. Smaller RMSE values indicate a stronger influence of the parameter on the model's predictive accuracy. Figure 11 and Figure 12 depict the RMSE distributions of prediction results for various input parameters in the ConvLSTM-multi and PredRNN-multi models. SST exhibits the most pronounced influence on model accuracy, followed by date. Notably, the PredRNN-multi model displays heightened sensitivity to SST in comparison to the ConvLSTM-multi model, while the ConvLSTM-multi model exhibits greater sensitivity to date than the PredRNN-multi model. Additionally, the ConvLSTM-multi model demonstrates a slight but noticeable sensitivity to parameters like T2M, SKT, SSRD, and MSL, with slightly more sensitivity to these parameters than to the U10/V10. In contrast, the PredRNN-multi model shows consistent sensitivity to T2M, SKT, SSRD, MSL, and U10/V10.

Comparison of RMSE distributions for the first-day SIC predictions reveals insights into the influence of input parameters on model predictions. This assessment is conducted under three scenarios: when no noise is introduced into the two models, when reanalysis data are treated as noise, and when SIC is considered as noise, Figure 13 and Figure 14 illustrate the outcomes. When meteorological data are treated as noise, both the ConvLSTM-multi and PredRNN-multi models exhibit challenges in accurately predicting the extent of thin ice near the sea ice edge. Notably, the PredRNN-multi model is more affected, aligning with the prior observation that this model displays heightened sensitivity to SST. On the other hand, the response of these models differs when input SIC are noisy. In the case of the ConvLSTM-multi model, the ice-free region's characteristics exert a more significant impact than those of the sea ice region. This outcome is attributed to the model's capacity to learn from the SIC, subsequently isolating the ice-free regions. In contrast, the PredRNN-multi model is primarily influenced by the sea ice region, with minimal impact observed in the ice-free region. This phenomenon arises from the model's ability to glean predictive insights from the distribution of sea ice within the SIC.

3.4. Prediction ability of the model under extreme conditions

Since the commencement of SIC data recording by the NSIDC, the Arctic's sea ice area in September reached historic lows in 2012 and 2020. To ensure the continuity and adequacy of the model's training data, the data from 2020 and 2021 are reserved for testing to evaluate the model's predictive performance in both extreme and normal years. In Figure 4a and 4b, the prediction accuracy of the four models for 2020 is behind that for 2021, with the ConvLSTM-multi model exhibiting superior predictive capabilities for the abrupt changes in September 2020's sea ice area. Figure 5 and Figure 6reveal a decrease in the grid count for high-value SIC from June to September 2020, with values approaching 0%. Notably, low-value SIC in 2020 was significantly less than in 2021. Consequently, metrics such as normalized nRMSE, ACC, NSE, and SSIM, which are sensitive to sea-ice edge prediction accuracy, performed better in 2020 compared to 2021 (Figure 4c-f). ComparingFigure 5d and Figure 6d, the PredRNN-multi model tends to overestimate SIC values from June to September 2020, primarily in the range of 95% to 100%, while its underestimating of values tends to be clustered around 0%. This indicates that the PredRNN-multi model falls short in predicting the sharp decline in the sea ice area in 2020, which also elucidates the significant difference in IIEE values between its 2020 and 2021 predictions (Figure 9). Conversely, the ConvLSTM-multi model's IIEE values for SIC predictions in 2020 and 2021 closely align with each other and exhibit less sensitivity to extreme years. Additionally, in Figure 11, Figure 12, Figure 13 and Figure 14, the trends in prediction accuracy corresponding to various input parameters in relation to prediction duration and spatial distribution remain consistent across different years, suggesting that the influence of input parameters on the model remains largely independent of yearly variations.

4. Conclusions

This study explores the integration of multiple predictors into the ConvLSTM and PredRNN models, thereby creating ConvLSTM-multi and PredRNN-multi models. These models are designed for the 10-day prediction of SIC across the entire Arctic Ocean. The findings underscore that the incorporation of meteorological reanalysis data within the ConvLSTM and PredRNN models leads to a significant enhancement in the daily scale prediction accuracy of SIC. Notably, the predictive accuracy of the model is most influenced by SST, followed by the date of SIC. The ConvLSTM-multi model demonstrates the lowest MAE, RMSE and IIEE, with the least increase with the extension in the prediction time. Moreover, the ConvLSTM-multi model exhibits commendable accuracy, particularly during extreme years, albeit with a slightly inferior performance in predicting the distribution of SIC values compared to the PredRNN-multi model. The evaluation of sea-ice edge displacement indicates that the ConvLSTM-multi model excels during summer and autumn, while the PredRNN-multi model performs optimally during spring and winter. Consequently, this comparative analysis highlights that the ConvLSTM-multi model is better suited for predicting SIC over the next 10 days.

Despite the progress made in this study, it is essential to acknowledge several limitations inherent to the model. The predictive capacity for SIC hinges on both the quantity of available SIC datasets and the diversity of input parameters. However, the quantity of SIC datasets and the spatial resolution offered by NSIDC are constrained, necessitating a higher volume of datasets to further refine the model's parameters. Additionally, the model exhibits limited predictive capability for summer sea ice due to its exclusive reliance on reanalysis data. Crucial parameters such as sea ice thickness and ice melting pools, which wield a substantial influence on sea ice dynamics, have not been incorporated into the model. This omission can be attributed to the limited availability and accessibility of these parameters. Lastly, the study falls short in comprehensively analyzing the impact of input parameters on prediction outcomes and discerning the intricate nonlinear relationships between them. Future research endeavors will focus on intensifying the analysis of physical parameters influencing sea ice variations and enhancing the model's capacity to assimilate input data, to enhance the prediction ability of the model.