As shown in Figure 1, the Beaufort Sea region in the Arctic (longitude: 135° to 168°W, latitude: 72° to 79°N) was selected as the study area. This study site has a high latitude, and the cold Beaufort Current ensures the presence of extensive annual sea ice and multi-year sea ice in the area (Wiebke et al., 2020; Timmermans and Toole, 2023), which can fully verify the performance of the models.

The Soil Moisture and Ocean Salinity (SMOS) satellite, launched by the European Space Agency (ESA), carries an L-band two-dimensional synthetic aperture microwave radiometer that is highly sensitive to the emissivity characteristics of sea ice with varying thicknesses. This enables long-term, large-scale, and high-frequency observations of Arctic sea ice thickness (Kaleschke et al., 2012). The SMOS Level 3 Sea Ice Thickness (SMOS L3 SIT) product released by ESA provides both high temporal and high spatial resolution, making it one of the few gridded sea ice thickness datasets that simultaneously offer daily temporal resolution and a spatial resolution of 12.5 km, Table 1. Such characteristics allow the dataset to capture finer-scale temporal variations and spatial patterns of sea ice thickness, which is particularly valuable for the present study.

Moreover, the observational accuracy of the SMOS L3 SIT product has been extensively validated and widely recognized (Kaleschke et al., 2016), and it has been successfully applied in various contexts, including data assimilation in numerical models (Xie et al., 2016; Gupta et al., 2021). ESA has been providing daily sea ice thickness observations from SMOS since the 2010 freeze-up season, and the dataset is publicly accessible. In this study, SMOS L3 SIT data from October 15 to January 15 of each year were collected to support the short-term prediction of Arctic sea ice thickness.

The atmosphere, ocean, and sea ice are interacted and influence each other, making the changes in sea ice thickness not spatially independent. Current research is dedicated to characterizing the local variations of SIT, while the redundancy of observational data in terms of spatial distribution characteristics and temporal trends poses a challenge to the description of the global features of SIT. Spatiotemporal decomposition can decompose the SIT data to identify and separate time and space information, aiding in understanding the intrinsic mechanisms of sea ice thickness and describing its variation characteristics.

As shown in Figure 2a, Empirical Orthogonal Function (EOF) is a spatiotemporal analysis method (Deser et al., 2000) for extracting the main patterns and characteristics from SIT data. It projects the SIT data onto the covariance matrix through covariance matrix calculations and decomposes the time-varying SIT data into time-invariant EOFs and time-dependent PCs through Eigen Value Decomposition (EVD).

Assuming S represents the number of SIT spatial points on a given day, and T represents the number of time series at a certain location, then the SIT data can be represented as X˜S×Tob. Eliminating the influence of long-term states in the SIT variations, the SIT anomaly matrix XS×Tob is calculated using Equation 1:

Then XS×Tob can be expressed as Equation 2:

in which, the element xst represents the SIT data at the t-th time point and the s-th spatial point.

From XS×Tob, the SIT covariance matrix CS×S can be obtained using Equation 3:

Subsequently, through EVD, the SIT spatial eigenvectors EOFS×S and their corresponding time coefficients (λ1,…,λS) are obtained using Equation 4:

in which, ΛS×S is a diagonal matrix with the time coefficients (λ1,…,λS) arranged in descending order, as in Equation 5:

Ultimately, by matrix multiplication, the time components PCS×T corresponding to EOFS×S are obtained using Equation 6:

The SIT observation at a specific moment is viewed as a linear combination of various spatial modes with different weights. The spatiotemporal decomposition module decomposes the SIT data into EOFs and PCs. Consequently, the variation in SIT is transformed into spatial modes that do not change with time and different series components containing spatiotemporal information, which aids in describing its underlying mechanisms.

Characterizing the temporal variation characteristics of SIT is essential for obtaining accurate daily-scale spatiotemporal forecasts of SIT, yet current research lacks consideration of contextual features in SIT temporal variation. Compared to other deep learning models, the Transformer can build associations between any two elements in the sequence through attention mechanisms, characterizing the variation characteristics of SIT sequences without being affected by the length of the dataset (Ashish et al., 2023). As shown in Figure 2b, the Transformer is introduced and used to construct contextual features of SIT PCs in the temporal domain. The PCs derived from SIT data are input into the Transformer, and the predicted PCs are obtained through the encoder and decoder.

Figures 3a, b illustrate the Encoder-Decoder architecture of the Transformer. The Encoder is used to generate the contextual representation of the spatiotemporal variation of SIT for the Decoder. The Encoder consists of multiple identical layers, each containing a multi-head attention mechanism and a feed-forward neural network. Each sublayer is followed by residual connections and layer normalization, which helps alleviate the vanishing gradient problem and stabilize model training. The Decoder predicts the PCs for future moments based on the generated SIT contextual representation and SIT PCs. Similarly, the Decoder is composed of multiple identical layers, with the main sublayers including a masked multi-head attention mechanism, a multi-head attention mechanism, and a feed-forward neural network.

Establishing dependencies between different parts of the SIT PCs in the temporal sequence and capturing relationships between different locations rely on the Multi-Head (MH) attention mechanism. Figure 3c shows the MH attention mechanism used, which includes multiple parallel attention mechanisms to capture different features and relationships within the input SIT data, as in Equation 7:

in which, WH is the weight matrix, and the nth attention mechanism in the cascade, h e a d n = A t t e n t i o n ( Q i c e n , K i c e n , V i c e n ) , can be expressed using Equation 8:

in which, 1dk is the scaling factor, and the elements of the attention mechanism, query Qice, key Kice, and value Vice, are obtained by passing the SIT PCs through three distinct linear transformation layers.

Therefore, when using the Transformer for spatiotemporal prediction of SIT, the SIT PCs with positional encoding are input into the encoder, as in Equation 9:

MH−Attention focuses on the important locations of the SIT PCs input by exploring the complex dependencies between different positions, thereby generating a contextual representation of the sea ice sequence, as in Equation 10:

The feed-forward neural network applies non-linear transformations to Contextice, enhancing the expressive power of Transformer over SIT PCs, as shown in Equation 11:

In the decoder, Masked Multi-Head (MMH) attention is used to gather contextual information from the already generated PCs as in Equation 12:

in which, Targetsice represents the PCs that have been predicted in the sequence to be forecasted, and tar_Contextice represents the contextual representation of the PCs that have been predicted in the sequence to be forecasted. Compared to MH attention, MMH attention ensures that transformer can only access information up to the current moment when processing the contextual information of PCs, and cannot obtain information from future moments, thereby preserving the ability of model to make temporal predictions.

Further utilize MH−Attention and FFN to generate the prediction results for the PCs in next future moment from nolin_Contextice and tar_Contextice:

Iterating Equations 12 and 13, the predicted PCs at each moment in the sequence to be predicted, PC_pred, are passed through the output layer to obtain the final results of the SIT PCs, PC_output, as shown in Equation 14:

Through the aforementioned process, the Transformer is able to utilize the contextual information contained within the SIT PCs for the prediction of SIT.

The PCs obtained from the spatiotemporal prediction module contain information about the future spatiotemporal changes in sea ice thickness, while the spatial distribution patterns of sea ice thickness are embedded in the EOFs output by the spatiotemporal decomposition module. As shown in Figure 2c, the inverse principle of EOFs is used to reconstruct the sea ice thickness prediction results. In actual calculations, multiplying the spatial EOFs by the temporal PCs and adding the baseline 1 T ∑ t = 1 T X ˜ S × T o b , the SIT prediction result X˜S×Tpred can be obtained as in Equation 15:

The SMOS L3 SIT dataset includes data from 14 years, and cross-validation is used to test the predictive performance of models across different years, see Supplementary Material for more details. The spatiotemporal changes of sea ice thickness from 2022–2023 were analyzed in greater detail, where data from 2010–2019 were used as the training set, data from 2020–2021 were used as the validation set, and data from 2022–2023 were used as the test set. The mean squared error (MSE) was used as the loss function and minimized, iterating over 1000 time steps and optimizing EOF-Trans parameters through backpropagation of gradients, as shown in Equation 16:

in which, M is the number of training samples, predi and truei are the predicted and observed values for the i-th sample, respectively.

The performance is evaluated using the root mean square error (RMSE), as in Equation 17, and Pearson correlation coefficient (Corr), as in Equation 18, between the observed and predicted values:

in which N is the number of forecast points in the study area, true¯=1N∑i=1Ntruei and pred¯=1N∑i=1Npredi represent the mean values of the true and predicted values, respectively. At the same time, the Persistence model, Climatology model, TOPAZ model, EOF-Markov (Wang et al., 2023), CNN (Kim et al., 2020), ConvLSTM (He et al., 2022) and U-Net (Andersson et al., 2021) are introduced to validate the performance of EOF-Trans relative to classical models. The hyperparameter settings are provided in the Supplementary Materials.

The number of principal components plays a critical role in the prediction performance of EOF-Trans, so this study first conducts a principal component sensitivity analysis. The explained variance ratios corresponding to principal components characterize the spatiotemporal modes and their contributions to sea ice thickness changes, and the quantity of principal components determines the scale and hierarchy of the variability covered, as shown in Figure 4a. By using varying numbers of PCs to predict sea ice thickness for 2022–2023, this study found that the prediction accuracy follows a U-shaped trend with the number of PCs, as illustrated in Figure 4b. When the number of PCs is 18, they include 95% of the sea ice evolution information, yielding the most accurate prediction with an error of 0.1716 m. When the number of PCs is fewer than 18, the prediction accuracy improves relatively rapidly with an increase in PCs, as a greater number of PCs contain more useful sea ice information, positively contributing to the predictive performance. However, when the number of PCs exceeds 18, the prediction accuracy declines with further increases in PCs, likely due to unpredictable small-scale variability interfering with the predictions. These findings demonstrate that the number of PCs is crucial for the predictive accuracy of the EOF-Trans. Therefore, this study selects 18 PCs for predicting the spatiotemporal variations in sea ice thickness.

As shown in Figures 5a, c, compared to the TOPAZ numerical model, Persistence and Climatology baseline models, the EOF-Trans consistently achieves the lowest root mean square error and the highest correlation coefficient in each year. This indicates that the predictions from EOF-Trans are closest to the actual observations and exhibit the best statistical correlation, fully demonstrating the superior performance in sea ice thickness prediction. As a classical numerical model, TOPAZ also demonstrates strong predictive capability, with its RMSE remaining stable between 0.20–0.25 m in most years and Corr ranging from 0.65 to 0.75. While both Persistence and Climatology provide basic forecasting skill, Climatology performs significantly better than Persistence. Figures 5b, d further illustrate the prediction errors and correlation performance of EOF-Trans under different leadtime in 2022–2023. Benefiting from the ability to capture long-term temporal dependencies in Transformer, EOF-Trans exhibits strong stability. Although TOPAZ performs well in the initial lead days, its correlation declines significantly due to error accumulation. The Persistence model, which assumes that the current ice thickness distribution resembles the previous time step, performs particularly well in the first two lead times. However, its predictive skill deteriorates rapidly when significant changes in sea ice thickness occur. In contrast, Climatology reflects the average trend of ice thickness distribution, resulting in relatively stable but less dynamic predictions. In summary, EOF-Trans delivers the best overall performance, followed by Climatology, with TOPAZ and Persistence ranking successively lower.

Figure 6 provides spatial and temporal examples comparing the differences between EOF-Trans prediction results and actual observations. The predicted results exhibit distinct spatial characteristics. As latitude increases, the SIT predicted by EOF-Trans gradually increases to 1.7 m. Thanks to the descriptive ability of EOFs for the overall characteristics of ice thickness, the results of EOF-Trans are relatively close to actual observations and can accurately predict the spatial distribution of SIT. From a temporal perspective, the prediction errors on the 7th, 14th, and 21st days gradually increase. Although the errors in EOF-Trans predictions increase with time, the spatial distribution of SIT it depicts remains consistent with actual observations. The description of the temporal context of SIT PCs by EOF-Trans ensures that its performance is still reliable even with larger lead days.

Although EOF-Trans outperforms both the baseline models and the numerical model, its relative advantages over other deep learning models require further comparison. Using cross-validation, this study evaluates the interannual variations and lead-time-dependent changes in the RMSE and correlation of EOF-Trans and four additional models, EOF-Markov, ConvLSTM, CNN, and U-Net, as shown in Figure 7. EOF-Trans consistently achieves the lowest RMSE and highest correlation across all years in Figures 7a, b, demonstrating the best overall prediction performance. Despite the increased prediction difficulty in 2010 and 2018, when all models exhibit RMSE values exceeding 0.25 due to strong Beaufort High anomalies (Kwok and Cunningham, 2015; Petty, 2018), EOF-Trans still maintains the lowest error among all models. U-Net generally performs better than ConvLSTM, suggesting that a deeper stacked convolutional structure provides advantages over the combination of convolution and recurrent memory mechanisms for sea ice thickness prediction. However, U-Net lacks explicit temporal context modeling, which may be an important reason why EOF-Trans outperforms U-Net in predictive skill. EOF-Markov and CNN also exhibit certain predictive capability but do not reach the performance level of EOF-Trans.

Further examination of the prediction performance of EOF-Trans and the deep learning models across different leadtime in Figures 7c, d reveals a high degree of consistency with the interannual comparison. The Transformer architecture enables effective associations between information at different time steps, allowing EOF-Trans to maintain relatively high accuracy even at extended leadtimes, thereby demonstrating strong stability in Figure 7c. U-Net, which adopts a direct prediction strategy, effectively mitigates error accumulation with increasing leadtime. In contrast, the memory mechanism in ConvLSTM leads to continuous error accumulation, negatively impacting its predictive performance. This explains why ConvLSTM exhibits higher correlations than EOF-Trans during the first three leadtimes, but subsequently experiences a rapid decline Figure 7d.

To evaluate the predictive reliability of the models, the results of three independent runs of EOF-Trans and the comparative models on the 2022–2023 test dataset are reported in Table 2. The EOF-Trans model attains a mean value of 0.8804 with a standard deviation of 0.0078, indicating a high degree of consistency across repeated experiments. The U-Net and ConvLSTM models exhibit similar levels of variability, however, U-Net demonstrates superior predictive performance relative to ConvLSTM. The performance of EOF-Markov and CNN is comparatively lower.

Furthermore, the spatial distributions of prediction bias for EOF-Trans, ConvLSTM, and U-Net at different leadtime are shown in Figure 8. The prediction biases of all models are primarily concentrated in regions with strong sea ice variability, where they tend to underestimate thin ice and overestimate thick ice. This spatial pattern becomes more pronounced as the leadtime increases. However, EOF-Trans exhibits smaller bias magnitudes and more localized bias regions, indicating superior model stability. In contrast, ConvLSTM and U-Net fail to accurately capture the detailed evolution of sea ice, showing noticeable deviations as early as lead time 7, with the discrepancy further amplified at longer leadtimes. Although ConvLSTM demonstrates prediction skills comparable to EOF-Trans in certain local areas, its insufficient consideration of interregional dependencies leads to suboptimal reconstruction of the overall spatial distribution of ice thickness, indirectly highlighting the effectiveness of EOF-based global feature extraction. U-Net, being a fully convolutional architecture, relies entirely on convolutional operations for spatial feature extraction. While this enables relatively good spatial modeling capability, its limited receptive field results in notable prediction inconsistencies across different local regions, making its overall performance slightly inferior to that of EOF-Trans.

To further examine how EOF-Trans utilizes temporal information, this study conducts a visualization analysis of the multi-head attention weights in the Transformer decoder. Figure 9 presents the attention weight distributions of the three attention heads, based on statistics derived from the 2022–2023 test dataset. Overall, the multi-head attention mechanism enables EOF-Trans to attend to information across different temporal ranges, thereby forming a multi-level representation of short-term dependencies, intermediate lagged relationships, and longer-term temporal structures present in the spatiotemporal variability of SMOS-retrieved sea ice thickness.

Among the three attention heads, Head 1 exhibits a pronounced emphasis on recent inputs, with future predictions assigning the highest weights to the most recent time steps (i.e., t, t-1, and t-2), while contributions from earlier time steps gradually decrease and remain relatively low Figure 9a. This pattern indicates that the model places greater statistical importance on short-term temporal continuity and autocorrelation in the SMOS-retrieved sea ice thickness data, suggesting that this head primarily captures rapid variations at short temporal scales. In contrast, Head 2 shows a more distributed attention pattern, with relatively higher weights assigned to inputs from t-2 to t-4 Figure 9b. This behavior reflects the model’s ability to incorporate lagged temporal dependencies spanning several days, representing accumulated influences embedded in the SMOS-retrieved sea ice thickness data without explicitly resolving specific physical drivers. Head 3 displays the most uniform attention distribution across time steps Figure 9c, indicating a tendency to aggregate information over the entire temporal window of the SMOS-derived sea ice thickness sequence. Such global attention supports the extraction of broader temporal tendencies in the input sequence and contributes to maintaining stable predictions at longer lead times.

To assess the reliability of the model predictions, this study compares the kernel density distributions of the predicted RMSE and the inherent retrieval uncertainty of the SMOS L3 sea ice thickness product over the 2022–2023 test period Figure 10. As shown in Figure 10a, the RMSE distribution exhibits a unimodal pattern, with the density peak corresponding to an RMSE of approximately 0.23 m. More than 90% of the samples fall within the range of 0.10–0.25 m, indicating that the model’s prediction errors are generally low and exhibit limited variability.

In contrast, the uncertainty distribution of the SMOS L3 product Figure 10b presents a distinct bimodal structure. The first density peak occurs within the 0.0–0.5 m range, representing regions with relatively low retrieval uncertainty, while the second peak, located between 1.5 and 2.0 m, corresponds to ice conditions where the retrieval uncertainty is substantially higher. The comparison reveals that the RMSE of EOF-Trans is lower than the primary uncertainty level of the SMOS L3 data, suggesting a potential tendency to fit observational noise or biases. Nevertheless, when considered together with the preceding analyses, these results indicate that EOF-Trans primarily learns and reproduces the spatiotemporal regularities embedded in the SMOS-retrieved sea ice thickness fields. It is worth noting that the proposed framework is trained on satellite-retrieved sea ice thickness fields, and its learned representations therefore reflect the dominant spatiotemporal patterns present in the SMOS product. The extent to which these patterns can be directly associated with the underlying physical processes of sea ice thickness evolution remains uncertain and warrants further investigation.