CRITER 1.0: a coarse reconstruction with iterative refinement network for sparse spatio-temporal satellite data

Satellite observations of sea surface temperature (SST) are essential for accurate weather forecasting and climate modeling. However, these data often suffer from incomplete coverage due to cloud obstruction and limited satellite swath width, which requires development of dense reconstruction algorithms. The current state of the art struggles to accurately recover high-frequency variability, particularly in SST gradients in ocean fronts, eddies, and filaments, which are crucial for downstream processing and predictive tasks. To address this challenge, we propose a novel two-stage method CRITER (Coarse Reconstruction with ITerative Refinement Network), which consists of two stages. First, it reconstructs low-frequency SST components utilizing a Vision Transformer-based model, leveraging global spatio-temporal correlations in the available observations. Second, a UNet type of network iteratively refines the estimate by recovering high-frequency details. Extensive analysis on datasets from the Mediterranean, Adriatic, and Atlantic seas demonstrates CRITER's superior performance over the current state of the art. Specifically, CRITER achieves up to 44 % lower reconstruction errors of the missing values and over 80 % lower reconstruction errors of the observed values compared to the state of the art.

1 Introduction

Infrared satellite sea surface temperature (SST) data are critical for ocean modeling, climate monitoring, fisheries management, and marine ecology (O'Carroll et al., 2019). On the one hand, the SST is a key boundary condition for atmospheric models extending from classical numerical weather prediction (Senatore et al., 2020; Chelton, 2005) to extreme storms (Ricchi et al., 2023) and climate variability (Garcia-Soto et al., 2021). In the ocean realm, continuous description of SST is vital for analyses of mesoscale (Bishop et al., 2017) and submesoscale baroclinic processes like fronts and eddies but also for implementations of atmosphere–ocean couplings through turbulent heat fluxes (Strajnar et al., 2019; Ličer et al., 2016). Furthermore, vertical temperature profiles are a critical driver of heat, carbon, and nutrient exchange between the surface and the deep ocean and thus for a wide plethora of biogeochemical processes (Mogen et al., 2022) in the ocean surface boundary layer which depend on the temperatures above the pycnocline. Last but not least, SST is a key parameter for the detection, mapping, and analysis of marine heatwaves (Hobday et al., 2016), and reconstructed satellite fields are imperative for determining the regional extent and intensity of such extreme events (Pastor and Khodayar, 2023; Darmaraki et al., 2019), which can have enormous impacts on aquaculture, fisheries, and other aspects of economy (Gómez-Gras et al., 2021; Garrabou et al., 2022).

Such downstream applications therefore often require complete, dense SST fields but cloud cover and sparse satellite coverage invariably lead to gappy and sparse data in both space and time. Reconstruction of gaps in the observations is therefore essential for a continuous description of ocean temperature fields and for many daily operational processes. These can be categorized into two groups: (i) extensions of the optimal interpolation (OI) scheme (Taburet et al., 2019; Ubelmann et al., 2021), and (ii) data-driven approaches. The latter includes methods based on empirical orthogonal functions (EOFs), such as DINEOF (Alvera-Azcárate et al., 2005), and, more recently, end-to-end deep learning techniques. Notable deep learning methods include DINCAE1 (Barth et al., 2020), dADRSR (Buongiorno Nardelli et al., 2022; Fanelli et al., 2024), TS-RBFNN (Young et al., 2024), DINCAE2 (Barth et al., 2022), 4DVarNet (Fablet et al., 2021), 4DVarNet-SSH (Beauchamp et al., 2023), the SSH reconstruction method by Martin et al. (2023), NeurOST (Martin et al., 2024), and MAESSTRO (Goh et al., 2024).

Traditional methods like DINEOF (Alvera-Azcárate et al., 2005) have been widely adopted, iteratively filling in missing data using truncated EOF decomposition. While effective for large-scale patterns, DINEOF struggles with fine-scale features, mostly because of their transient nature. Deep learning approaches have since emerged, surpassing traditional methods' performance. DINCAE1 (Barth et al., 2020) introduced a UNet-based (Ronneberger et al., 2015) model with probabilistic output, while 4DVarNet (Fablet et al., 2021) proposed an energy-based formulation for interpolation, achieving comparable SST reconstruction performance to a convolutional autoencoder architecturally similar to DINCAE1. Recently, Young et al. (2024) proposed a physically informed neural network that reconstructs daily SSTs in both cloudy and cloud-free areas, outperforming DINEOF. Beyond gap filling, super-resolution techniques have been developed to enhance SST resolution: Lloyd et al. (2021) designed a network that fuses optical and thermal satellite imagery, and, more recently, Fanelli et al. (2024) applied a convolutional super-resolution network (originally proposed by Buongiorno Nardelli et al., 2022) to super-resolve small low-resolution SST tiles obtained through optimal interpolation, improving fine-scale feature reconstruction.

DINCAE2 (Barth et al., 2022), the current state of the art and successor to DINCAE1, extended the original implementation with an additional refinement UNet. It operates on temporally consecutive partial SST observations, gradually improving central SST field reconstruction. However, its finite receptive field limits long-range spatio-temporal dependency exploitation, resulting in oversmoothed reconstructions lacking high-frequency details. Recently, MAESSTRO (Goh et al., 2024) addressed some limitations by adapting the Masked Autoencoder (MAE) (He et al., 2022) framework for SST reconstruction. It employs a Vision Transformer (ViT) (Dosovitskiy et al., 2021) architecture to capture global spatial dependencies. However, its single-time-step approach neglects temporal correlations, potentially compromising reconstruction quality for large, contiguous cloud occlusions. Furthermore, MAESSTRO's random patch masking strategy during training and evaluation may inadequately represent real cloud patterns, potentially yielding optimistic error estimates.

To address these limitations, we propose a two-stage Coarse Reconstruction with ITerative Refinement network (CRITER). A transformer-based module first leverages long-range spatio-temporal dependencies to estimate a low-frequency reconstruction. Subsequently, an iterative refinement module enhances high-frequency content. Unlike previous methods, which attempt full signal reconstruction in each block, CRITER decomposes the problem into a sequence of networks, each reducing the residual error of its predecessor, thus optimizing network capacity for local error reduction.

The paper is structured as follows. Section 2 contains descriptions of employed datasets together with preprocessing steps executed prior to the training. Section 3 describes the CRITER architecture, focusing on coarse reconstruction step in Sect. 3.1, its iterative refinement in Sect. 3.2, and residual estimation network in the “Residual estimation network (REN)” section. The training strategy is described in Sect. 3.3, and results are listed in Sect. 4, including an in-depth ablation study investigating the role of individual architectural components (Sect. 4.5).

2 Input data: sea surface temperature

2.1 Evaluation datasets

For our study we utilize Level 3 (L3) sea surface temperature (SST) satellite observation products. L3 level of product refers to the satellite product where spatially sparse and irregular point observations of the ocean surface are gridded into a fixed grid across space and/or time. Such products may combine multiple satellite overpasses or even multiple sensors for the same observed quantity.

Specifically we consider the following three datasets corresponding to three different geographic regions:

1. Central Mediterranean: The SST_MED_SST_L3S_NRT_OBSERVATIONS_010_ 012_a (E.U. Copernicus Marine Service Information, 2023a) dataset contains daily near-real-time (NRT) SST measurements over the Mediterranean sea from 1 January 2008 to 31 December 2021. The dataset is provided on a remapped grid with a spatial resolution of 0.0625° × 0.0625°.

2. Adriatic: The SST_MED_PHY_L3S_MY_010 _042 (Pisano et al., 2016; Casey et al., 2010) dataset contains daily multi-year reprocessed (MY) SST measurements over the Adriatic Sea from 25 August 1981 to 31 December 2022.

   The dataset is provided on a remapped grid with a spatial resolution of 0.05° × 0.05°.

3. Atlantic: The SST_ATL_PHY_L3S_MY_010 _038 (E.U. Copernicus Marine Service Information, 2023b) dataset contains daily multi-year reprocessed (MY) SST measurements from 1 January 1982–1 January 2022. The dataset is provided on a remapped grid with a spatial resolution of 0.05° × 0.05°.

These regions were chosen due to their oceanographic variety. The Adriatic is an elongated semi-enclosed basin with correspondingly poor satellite coverage, and the Central Mediterranean exhibits a wide variety of oceanographic regimes (from regions of freshwater influence in the northern Adriatic to a much deeper Ionian where Levantine and Adriatic water masses communicate), while the Atlantic region is essentially an open ocean region, very different from the Adriatic. These regions should demonstrate generalization abilities of CRITER under a variety of oceanographic conditions. The geographic areas of the three datasets are shown in Fig. 1. It is worth noting that two different satellite products are used in this study, a near-real-time (NRT) and a multi-year (MY) reprocessed dataset. This was done to show that like DINCAE2, CRITER also generalizes well across various datasets of SST. Furthermore, multi-year reprocessed datasets come at a higher resolution and span significantly longer periods of time, which gives access to a larger train and, more importantly, test set.

Figure 1 The map shows the spatial extent of the Central Mediterranean, Adriatic, and Atlantic datasets, highlighting the distinct geographic areas covered by each dataset.

2.2 Input data preprocessing

2.2.1 Filtering out days with excessive cloud coverage

The satellite products corresponding to Level 3 (L3) SST are provided on a fixed grid but are spatially sparse over a subset of spatial locations (mainly due to clouds and land pixels). For training and evaluation of the method in this work, additional missing values need to be simulated to test network performance on values which are hidden to the network but are otherwise known. If the original SST observation field already contains a large number of missing measurements, it becomes difficult to effectively simulate additional missing data. Consequently, to ensure that the dataset is suitable for training and evaluating models, observations that are too sparse need to be filtered out. In the preprocessing stage, we first construct sequences of 3 temporally consecutive days of observed SST fields, as proposed by Barth et al. (2020). The observation sequences are then filtered. Specifically, any 3 d observation sequence is discarded according to the following rule: if the cloud coverage, defined as the fraction of pixels that are missing in the central observation field, relative to the total number of pixels belonging to the sea, is greater than or equal to a certain threshold, the corresponding observation sequence is discarded. The appropriate threshold is selected by considering the total number of samples in each dataset. Specifically, we use a threshold of 100 % for the Mediterranean dataset, resulting in a total of 5114 samples. For the Adriatic dataset, we apply a threshold of 60 %, which yields 7800 samples. Finally, we use a threshold of 75 % for the Atlantic dataset, resulting in 3454 samples.

2.2.2 Train, validation, and test datasets

The filtered satellite SST observations are chronologically split into three subsets: the train set, which comprises the first 90 % of the samples; the validation set, which comprises the next 5 % of the samples; and the test set, which consists of the last 5 % of the samples. The models are trained on the train set, the hyper-parameters are tuned on the validation set, and the performance is assessed on the test set. This approach ensures evaluation on future, unseen data with no temporal overlap between training and test phases.

3 CRITER – Coarse Reconstruction with ITerative Refinement network

Given a sequence of spatially sparse sea surface temperature observations ${\mathbf{X}}_m=[{\mathit{x}}_{t-{\mathrm{\Delta }}_t},\mathrm{\dots },{\mathit{x}}_t,\mathrm{\dots },{\mathit{x}}_{t+{\mathrm{\Delta }}_t}]$, where ${\mathit{x}}_t\in {\mathbb{R}}^{\mathrm{1}\times W\times H}$ is the potentially sparse observation field of width W and height H at time step t and $[t-{\mathrm{\Delta }}_t,t+{\mathrm{\Delta }}_t]$ defines the observed time interval, the task is to estimate the dense reconstruction $\stackrel{\mathrm{̃}}{\mathit{x}}$ at time step t and the uncertainty specified by the variance σ². Following Barth et al. (2022), we set the temporal horizon to Δ_t=1 d, thus in reconstruction of x_t, the days before and after day t are considered.

The proposed Coarse Reconstruction with ITerative Refinement network (CRITER) is a two-stage method composed of a Coarse Reconstruction Module (CRM), described in Sect. 3.1, and an Iterative Refinement Module (IRM), described in Sect. 3.2. An overview of the architecture is provided in Fig. 2.

Figure 2 Given observations for 3 consecutive days $[{\mathit{x}}_{t-\mathrm{1}},{\mathit{x}}_t,{\mathit{x}}_{t+\mathrm{1}}]$ and a binary mask M_t indicating missing pixels, CRITER densely reconstructs x_t in two phases. First, the CRM estimates a coarse reconstruction ${\widehat{\mathit{x}}}_t$, which the IRM then iteratively refines to produce the final reconstruction $\stackrel{\mathrm{̃}}{\mathit{x}}$ and uncertainty σ². CRM tokenizes the input into tokens requiring reconstruction T_r and contextual tokens T_c. These contextual tokens are encoded by a ViT-based encoder into ${\mathbf{T}}_{\mathrm{c}}^{\mathcal{E}}$, combined with T_r, and decoded by a ViT-based decoder into decoded tokens ${\mathbf{T}}_{\mathrm{t}}^{\mathcal{D}}$, which are finally mapped to ${\widehat{\mathit{x}}}_t$. In the IRM, dashed lines indicate the iterative refinement process. At each iteration i, the current reconstruction estimate ${\stackrel{\mathrm{̃}}{\mathit{x}}}^{\left(i\right)}$ and uncertainty estimate σ²⁽ⁱ⁾ are refined by adding the predicted residuals: reconstruction residual ${\mathit{\delta }}_x^{\left(i\right)}$ and uncertainty residual ${\mathit{\delta }}_{{\mathit{\sigma }}^{\mathrm{2}}}^{\left(i\right)}$. The index in REN⁽ⁱ⁾ indicates the change in network parameters in each iteration.

3.1 Coarse Reconstruction Module (CRM)

The Coarse Reconstruction Module (CRM, Fig. 2) follows the ViT encoder–decoder architecture (Dosovitskiy et al., 2021), similar to spatio-temporal MAE (Feichtenhofer et al., 2022). The input observation fields ${\mathbf{X}}_m=[{\mathit{x}}_{t-\mathrm{1}},{\mathit{x}}_t,{\mathit{x}}_{t+\mathrm{1}}]\in {\mathbb{R}}^{\mathrm{3}\times \mathrm{1}\times W\times H}$ are first fed to a tokenization process. To encode information about the yearly temperature cycle, each observation field x_t is concatenated channel-wise with a day-of-the-year auxiliary tensor ${\mathit{a}}_{\mathrm{t}}=\left[\sin \right(d_t\frac{\mathrm{2}\mathit{\pi }}{\mathrm{365}}),\cos (d_t\frac{\mathrm{2}\mathit{\pi }}{\mathrm{365}}\left)\right]\in {\mathbb{R}}^{\mathrm{2}\times W\times H}$, where the two channels contain constants, and d_t is the numerical day of year index (between 1 and 365). The resulting fields are split into non-overlapping $\mathrm{3}\times \mathrm{8}\times \mathrm{8}$ patches, which are then flattened and linearly projected into tokens of shape 1×D_t, where D_t is the dimension of tokens used in ViT blocks, thus creating the list of tokens $\mathbf{T}=\mathit{\{}{\mathbf{T}}_{\mathrm{r}},{\mathbf{T}}_{\mathrm{c}}\mathit{\}}$. Tokens T_r correspond to patches in x_t with at least one unobserved pixel and thus have to be reconstructed. Tokens T_c are the remaining tokens, and they are used as a context for reconstruction. To encode the extent of missing values in a token, all tokens in x_t are summed with their corresponding mask tokens. These are obtained by splitting the binary mask indicating missing pixels ${\mathbf{M}}_{\mathrm{t}}\in \mathit{\{}\mathrm{0},\mathrm{1}{\mathit{\}}}^{W\times H}$ into 8×8 non-overlapping patches, which are then flattened and projected into mask tokens of shape 1×D_t. To maintain the necessary spatio-temporal location of each token, all tokens in T are summed with a spatio-temporal positional embedding as in Feichtenhofer et al. (2022).

After obtaining tokens T, the context tokens T_c are encoded by a ViT (Dosovitskiy et al., 2021) encoder ℰ_ViT into ${\mathbf{T}}_{\mathrm{c}}^{\mathcal{E}}$ (Fig. 2). Then, the list of tokens T_r requiring reconstruction is concatenated with the list of the encoded tokens ${\mathbf{T}}_{\mathrm{c}}^{\mathcal{E}}$. The set of all tokens is again summed with the spatio-temporal positional embedding and passed through a ViT decoder 𝒟_ViT, producing the decoded tokens T^𝒟. The decoded tokens not corresponding to the central observation x_t are removed from T^𝒟, resulting in ${\mathbf{T}}_{\mathrm{t}}^{\mathcal{D}}$. Tokens in ${\mathbf{T}}_{\mathrm{t}}^{\mathcal{D}}$ are then linearly projected into 1×8² vectors and reshaped into $\mathrm{1}\times \mathrm{8}\times \mathrm{8}$ patches. Finally, the patches are reassembled into a grid to form the coarse reconstruction ${\widehat{\mathit{x}}}_t$. All pixel values corresponding to land areas are set to zero using the land mask ${\mathbf{M}}_{\mathrm{l}}\in \mathit{\{}\mathrm{0},\mathrm{1}{\mathit{\}}}^{W\times H}$ that accompanies the data.

3.2 Iterative refinement module (IRM)

To improve the reconstruction accuracy, the coarse reconstruction ${\widehat{\mathit{x}}}_t$ is refined by an iterative refinement module (IRM, Fig. 2) through a sequence of residual improvements, producing the final reconstruction $\stackrel{\mathrm{̃}}{\mathit{x}}$ and the corresponding uncertainty characterized by the variance σ². Per pixel j, we model the reconstructed SST as a Gaussian distribution parameterized by predicted mean ${\stackrel{\mathrm{̃}}{\mathit{x}}}_{\left(j\right)}$ and standard deviation σ_(j), following Barth et al. (2020). Note that σ² emerges from training the model to minimize Eq. (4), which penalizes over- and underestimation of the error variance σ².

Let ${\stackrel{\mathrm{̃}}{\mathit{x}}}^{\left(i\right)}$ and σ²⁽ⁱ⁾ be the reconstruction of the observation field x_t and its estimated uncertainty at ith refinement iteration. An iteration of IRM proceeds as follows. The input observation fields ${\mathbf{X}}_m=[{\mathit{x}}_{t-\mathrm{1}},{\mathit{x}}_t,{\mathit{x}}_{t+\mathrm{1}}]\in {\mathbb{R}}^{\mathrm{3}\times W\times H}$ and the refined estimates $[{\stackrel{\mathrm{̃}}{\mathit{x}}}^{\left(i\right)},{\mathit{\sigma }}^{\mathrm{2}\left(i\right)}]\in {\mathbb{R}}^{\mathrm{2}\times W\times H}$ from the previous iteration are concatenated channel-wise and passed to a residual estimation network REN⁽ⁱ⁾ (detailed in Sect. “Residual estimation network (REN)”) alongside the tokens ${\mathbf{T}}_{\mathrm{t}}^{\mathcal{D}}$ produced by CRM, to produce a two-channel output ${\mathbf{Y}}^{\left(i\right)}=[{\mathbf{Y}}_{\mathrm{1}}^{\left(i\right)},{\mathbf{Y}}_{\mathrm{2}}^{\left(i\right)}]\in {\mathbb{R}}^{\mathrm{2}\times W\times H}$. Following the formulation of Barth et al. (2020), ${\mathbf{Y}}^{\left(i\right)}=[{\mathbf{Y}}_{\mathrm{1}}^{\left(i\right)},{\mathbf{Y}}_{\mathrm{2}}^{\left(i\right)}]$ are decoded into reconstruction ${\mathit{\delta }}_x^{\left(i\right)}$ and uncertainty ${\mathit{\delta }}_{{\mathit{\sigma }}^{\mathrm{2}}}^{\left(i\right)}$ residuals:

$$\begin{array}{}\text{(1)}& {\displaystyle }{\mathit{\delta }}_{{\mathit{\sigma }}^{\mathrm{2}}}^{\left(i\right)}={\displaystyle \frac{\mathrm{1}}{\text{max}\left(\text{exp}\right(\text{min}({\mathbf{Y}}_{\mathrm{1}}^{\left(i\right)},{\mathit{\theta }}_{\mathrm{1}})),{\mathit{\theta }}_{\mathrm{2}})}},\text{(2)}& {\displaystyle }{\mathit{\delta }}_{\mathit{x}}^{\left(i\right)}={\mathbf{Y}}_{\mathrm{2}}^{\left(i\right)}\odot {\mathit{\delta }}_{{\mathit{\sigma }}^{\mathrm{2}}}^{\left(i\right)},\end{array}$$

where ⊙ denotes element-wise tensor multiplication (the Hadamard product), while θ₁ and θ₂, ${\mathit{\theta }}_{\mathrm{1}}>{\mathit{\theta }}_{\mathrm{2}}>\mathrm{0}$ are hyperparameters ensuring training stability. The reconstruction and uncertainty estimates at iteration i=0 are initialized with the coarse reconstruction ${\stackrel{\mathrm{̃}}{\mathit{x}}}^{\left(\mathrm{0}\right)}={\widehat{\mathit{x}}}_t$ and a zero σ²⁽⁰⁾=0. The reconstruction and uncertainty estimated at the (i+1)th refinement iteration are thus ${\stackrel{\mathrm{̃}}{\mathit{x}}}^{(i+\mathrm{1})}={\stackrel{\mathrm{̃}}{\mathit{x}}}^{\left(i\right)}+{\mathit{\delta }}_x^{\left(i\right)}$ and ${\mathit{\sigma }}^{\mathrm{2}(i+\mathrm{1})}={\mathit{\sigma }}^{\mathrm{2}\left(i\right)}+{\mathit{\delta }}_{{\mathit{\sigma }}^{\mathrm{2}}}^{\left(i\right)}$, respectively. IRM runs for N_IRM iterations, with each REN⁽ⁱ⁾ having its own set of trained parameters, allowing each to specialize to its respective residual estimation, finally producing the refined reconstruction $\stackrel{\mathrm{̃}}{\mathit{x}}$ and uncertainty σ².

Residual estimation network (REN)

The residual estimation network REN⁽ⁱ⁾ is a UNet-type architecture (Ronneberger et al., 2015). The encoder ℰ_REN takes the reconstruction and uncertainty estimates $[{\stackrel{\mathrm{̃}}{\mathit{x}}}^{\left(i\right)},{\mathit{\sigma }}^{\mathrm{2}\left(i\right)}]$ as well as the observation fields ${\mathbf{X}}_m=[{\mathit{x}}_{t-\mathrm{1}},{\mathit{x}}_t,{\mathit{x}}_{t+\mathrm{1}}]$ as input and produces the latent features z⁽ⁱ⁾ with an 8-fold reduction in spatial resolution compared to the input. The latent features are then enriched with spatio-temporally aggregated features ${\mathbf{T}}_{\mathrm{t}}^{\mathcal{D}}$ from CRM. Specifically, the tokens ${\mathbf{T}}_{\mathrm{t}}^{\mathcal{D}}$ (see Fig. 2) are spatially reshaped and bilinearly upsampled to match the dimensions of z⁽ⁱ⁾. The two tensors are concatenated and fused by the feature fusion module (FFM) (Yu et al., 2018), yielding the enriched bottleneck features ${\stackrel{\mathrm{̃}}{\mathit{z}}}^{\left(i\right)}$.

The resulting features are then input in the decoder 𝒟_REN and decoded to the same dimensions as the input $[{\stackrel{\mathrm{̃}}{\mathit{x}}}^{\left(i\right)},{\mathit{\sigma }}^{\mathrm{2}\left(i\right)}]$ via convolutional and upsampling blocks, while incorporating intermediate encoder features at multiple scales through UNet skip connections. The resulting decoded features are transformed with two 1×1 convolutional layers to produce a two-channel output ${\mathbf{Y}}^{\left(i\right)}\in {\mathbb{R}}^{\mathrm{2}\times W\times H}$.

3.3 Training strategy

CRITER is trained in two stages to train both CRM and IRM. First, supervised learning with automatically generated targets is used to train CRM. In this setup, part of the input signal is deleted, and the network is trained to reconstruct the entire input signal. The input training samples are created by sampling triples of consecutive observations $[{\mathit{x}}_{t-\mathrm{1}},{\mathit{x}}_t,{\mathit{x}}_{t+\mathrm{1}}]$ and deleting parts of the central observation x_t, resulting in $[{\mathit{x}}_{t-\mathrm{1}},{\mathit{x}}_t\odot {\mathbf{M}}_{\mathrm{m}},{\mathit{x}}_{t+\mathrm{1}}]$, where ${\mathbf{M}}_{\mathrm{m}}\in \mathit{\{}\mathrm{0},\mathrm{1}{\mathit{\}}}^{W\times H}$ is a generated binary mask with 0 corresponding to missing values. Following Barth et al. (2022), the masks M_m are generated by copying clouds from a random day not included in the triplet to maintain mask simulation realism. CRM is trained to minimize the following reconstruction error:

$$\begin{array}{}\text{(3)}& {\mathcal{L}}_{\mathrm{CRM}}={\displaystyle \frac{\mathrm{1}}{|{\mathbf{M}}_{\mathrm{t}}\odot {\mathbf{M}}_{\mathrm{l}}|}}{\sum }_{i=\mathrm{1}}^N\left[({\mathit{x}}_{t\left(i\right)}-{\widehat{\mathit{x}}}_{t\left(i\right)}{)}^{\mathrm{2}}{\mathbf{M}}_{\mathrm{t}\left(i\right)}{\mathbf{M}}_{\mathrm{l}\left(i\right)}\right],\end{array}$$

where ${\widehat{\mathit{x}}}_t$ is the coarse reconstruction generated by CRM, and mask M_t has zeros at locations where ground truth measurements within the observation field x_t are missing, while M_l has zeros at spatial locations belonging to land, and $|{\mathbf{M}}_{\mathrm{t}}\odot {\mathbf{M}}_{\mathrm{l}}|$ denotes the number of ground truth measurements. The summation goes over the N pixels in each of x_t, ${\widehat{\mathit{x}}}_t$, M_t, and M_l. The operator (⋅)_(i) indexes the ith element of a matrix. The consecutive observations used as the model input and the masks M_t, M_l, and M_m used in the training process are visualized in Fig. 3.

In the second stage, the parameters of CRM are fixed and only the parameters of IRM are trained. The training samples are generated as in CRM training, but since IRM produces the mean and variance of the reconstruction, the following negative log-likelihood loss is minimized as in DINCAE (Barth et al., 2020, 2022):

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{\mathcal{L}}_{\mathrm{IRM}}& ={\displaystyle \frac{\mathrm{1}}{|{\mathbf{M}}_{\mathrm{t}}\odot {\mathbf{M}}_{\mathrm{l}}|}}{\sum }_{i=\mathrm{1}}^N\\ & \left[{\displaystyle \frac{({\mathit{x}}_{t\left(i\right)}-{\stackrel{\mathrm{̃}}{\mathit{x}}}_{\left(i\right)}{)}^{\mathrm{2}}}{{\mathit{\sigma }}_{\left(i\right)}^{\mathrm{2}}}}+\log \left({\mathit{\sigma }}_{\left(i\right)}^{\mathrm{2}}\right)\right]{\mathbf{M}}_{\mathrm{t}\left(i\right)}{\mathbf{M}}_{\mathrm{l}\left(i\right)},\end{array}\end{array}$$

where $\stackrel{\mathrm{̃}}{\mathit{x}}$ and σ² are the reconstruction and variance estimated after the last iteration in IRM, and the summation goes over the N pixels in each of $\stackrel{\mathrm{̃}}{\mathit{x}}$, σ², and x_t. This loss thus trains the model to assign higher variance to areas with greater than expected reconstruction error. We validate the variance prediction quality in Sect. 4.4, by demonstrating its correlation with empirical errors.

Figure 3 Top row: a sequence of three consecutive observation fields ${\mathit{x}}_{t-\mathrm{1}},{\mathit{x}}_t,{\mathit{x}}_{t+\mathrm{1}}$, and the central observation x_t⊙M_m, with additional missing values deleted by the sampled mask M_m. Bottom row: the land mask M_l with zeros at land locations; the missing data mask M_t with zeros at locations with missing measurements in x_t; and M_m, which is a randomly sampled M_t from an observation field not included in the input.

3.4 Implementation details

CRM (Sect. 3.1) consists of 12 encoder and decoder transformer blocks, with 3 multi-head attention (MHA) heads, a token dimension of D_t=192, and a patch size of $\mathrm{3}\times \mathrm{8}\times \mathrm{8}$, where 3 denotes the number of channels, while 8×8 represents the width and height, respectively. IRM (Sect. 3.2) consists of a CNN-based encoder with 3 double conv blocks, each followed by a 2×2 max pooling operation. The double conv block is composed of two 3×3 convolutional layers, each followed by a batch normalization layer and a ReLU activation function. The number of convolutional kernels in each block is 32,64, and 128, respectively. This is followed by another double conv block, with 256 kernels, at the bottleneck of the network, a Feature Fusion Module (FFM), and a decoder with 3 transpose convolution layers, each followed by a concatenation based skip connection and a double conv block. The number of kernels in each block is 128,64, and 32, respectively. IRM utilizes N_IRM=3 refinement iterations – this value is selected based on the results of the ablation study in Sect. 4.5.4. Hyperparameters θ₁ and θ₂ are set as ${\stackrel{\mathrm{̃}}{\mathit{\theta }}}_{\mathrm{1}}=\text{ln}\left(N_{\text{IRM}}\right)+{\mathit{\theta }}_{\mathrm{1}}$ and ${\stackrel{\mathrm{̃}}{\mathit{\theta }}}_{\mathrm{2}}=N_{\text{IRM}}{\mathit{\theta }}_{\mathrm{2}}$ to ensure that the variance σ² is bounded between $\mathrm{1}/\text{exp}\left({\mathit{\theta }}_{\mathrm{1}}\right)$ and $\mathrm{1}/{\mathit{\theta }}_{\mathrm{2}}$ for an arbitrary number of refinement iterations N_IRM≥1.

4 Results

4.1 Implementation details

CRITER is implemented using the PyTorch library (Paszke et al., 2017) and trained on an NVIDIA Tesla V100 GPU. The CRM block is trained with a batch size of 8 using the AdamW optimizer with a learning rate $\mathit{\alpha }=\mathrm{3}\times {\mathrm{10}}^{-\mathrm{4}}$, β₁=0.9, and β₂=0.95 for 60 epochs (warm-up period) and then with a cosine decay scheduler (Loshchilov and Hutter, 2016) with step size 30 for another 140 epochs. In the next phase the IRM block is trained using the pre-trained CRM with fixed parameters. We train IRM using the Adam optimizer, with $\mathit{\alpha }=\mathrm{3}\times {\mathrm{10}}^{-\mathrm{4}}$, β₁=0.9, and β₂=0.999 for 300 epochs, using a step learning rate scheduler with step size 50 and multiplicative factor γ=0.5.

4.2 Performance measures

The performance of CRITER is assessed on an independent test set. Reconstruction quality is computed in terms of root-mean-squared error (RMSE) between the ground truth x_t and the reconstruction $\stackrel{\mathrm{̃}}{\mathit{x}}$. In particular, the overall reconstruction error RMSE_all is defined as

$$\begin{array}{}\text{(5)}& {\text{RMSE}}_{\text{all}}=\sqrt{{\displaystyle \frac{{\sum }_{i=\mathrm{1}}^N\left[({\mathit{x}}_{t\left(i\right)}-{\stackrel{\mathrm{̃}}{\mathit{x}}}_{\left(i\right)}{)}^{\mathrm{2}}{\mathbf{M}}_{\mathrm{t}\left(i\right)}{\mathbf{M}}_{\mathrm{l}\left(i\right)}\right]}{|{\mathbf{M}}_{\mathrm{t}}\odot {\mathbf{M}}_{\mathrm{l}}|}}}.\end{array}$$

For additional insights we compute the RMSE separately for (i) deleted regions, corresponding to observations artificially removed by simulated clouds in the L3 SST product and thus withheld during the training, and (ii) visible regions, corresponding to remaining observations post-deletion in x_t.

The reconstruction error of deleted regions is defined as

$$\begin{array}{}\text{(6)}& {\mathrm{\text{RMSE}}}_{\text{mis}}=\sqrt{{\displaystyle \frac{{\sum }_{i=\mathrm{1}}^N\left[({\mathit{x}}_{t\left(i\right)}-{\stackrel{\mathrm{̃}}{\mathit{x}}}_{\left(i\right)}{)}^{\mathrm{2}}{\mathbf{M}}_{\mathrm{t}\left(i\right)}{\mathbf{M}}_{\mathrm{l}\left(i\right)}(\mathbf{1}-{\mathbf{M}}_{\mathrm{m}\left(i\right)})\right]}{|{\mathbf{M}}_{\mathrm{t}}\odot {\mathbf{M}}_{\mathrm{l}}\odot (\mathbf{1}-{\mathbf{M}}_{\mathrm{m}}\left)\right|}}},\end{array}$$

M_m is the mask of deleted regions, and $|{\mathbf{M}}_{\mathrm{t}}\odot {\mathbf{M}}_{\mathrm{l}}\odot (\mathbf{1}-{\mathbf{M}}_{\mathrm{m}}\left)\right|$ denotes the number of deleted ground truth measurements. The reconstruction error of visible regions is defined as

$$\begin{array}{}\text{(7)}& {\text{RMSE}}_{\text{vis}}=\sqrt{{\displaystyle \frac{{\sum }_{i=\mathrm{1}}^N\left[({\mathit{x}}_{t\left(i\right)}-{\stackrel{\mathrm{̃}}{\mathit{x}}}_{\left(i\right)}{)}^{\mathrm{2}}{\mathbf{M}}_{\mathrm{t}\left(i\right)}{\mathbf{M}}_{\mathrm{l}\left(i\right)}{\mathbf{M}}_{\mathrm{m}\left(i\right)}\right]}{|{\mathbf{M}}_{\mathrm{t}}\odot {\mathbf{M}}_{\mathrm{l}}\odot {\mathbf{M}}_{\mathrm{m}}|}}},\end{array}$$

where $|{\mathbf{M}}_{\mathrm{t}}\odot {\mathbf{M}}_{\mathrm{l}}\odot {\mathbf{M}}_{\mathrm{m}}|$ is the number of visible ground truth measurements. To enhance the metric stability, we sample 10 distinct cloud masks for each test SST field, simulating realistic observational variability. We thus evaluate the performance on 2560, 3900, and 1720 masked SST fields for the respective regions, ensuring robust statistical validation.

4.3 Comparison with state of the art

We compare CRITER with DINCAE2 (Barth et al., 2022), a well-known and highly competitive SST reconstruction method, serving as a widely recognized benchmark in recent studies (Barth et al., 2024), and with the recently presented MAESSTRO (Goh et al., 2024) on the three datasets from Sect. 2.1. We reimplemented both DINCAE2 (originally in Julia) following Barth et al. (2022) and MAESSTRO (public implementation unavailable) following Goh et al. (2024) in Pytorch. To ensure a fair evaluation, both methods were trained using the same dataset splits, with tuned hyperparameters, and employed the same loss function computed over identical regions to CRITER. For MAESSTRO, architectural modifications were necessary to ensure comparability. Please refer to Appendix D for the implementation details of baseline models.

Results in Table 1 demonstrate CRITER's consistent superior performance across all datasets. Compared to the current state of the art, DINCAE2, CRITER achieves error reductions in deleted and visible regions of 20 % and 89 % for the Mediterranean, 44 % and 80 % for the Adriatic, and 1 % and 88 % for the Atlantic dataset, respectively. MAESSTRO's significantly lower performance is attributed to its single time step reconstruction approach. This hypothesis is confirmed by our ablation study, detailed in Sect. 4.5, which examines the importance of modeling spatio-temporal data dependencies.

The relative improvements of CRITER compared to the related methods vary across the datasets. This can be attributed to the differing amounts of information available for reconstruction, which is inversely proportional with the extent of missing values. Our analysis of missing values (Appendix A) reveals that the datasets can be ranked by the average amount of information available in each observation triplet, from highest to lowest: Mediterranean, Adriatic, and Atlantic. Notably, the Adriatic dataset shows the greatest decrease in reconstruction error, suggesting that CRITER achieves optimal improvement when the available information is moderate. In contrast, the Atlantic dataset, with the lowest amount of available information, likely requires additional data to be effectively reconstructed. To address this, we propose increasing the temporal horizon Δ_t and incorporating supplementary or proxy variables, such as chlorophyll a and surface winds. We leave the exploration of this approach to future work.

Table 1 Comparison of CRITER, DINCAE2, and MAESSTRO. We report the overall reconstruction error (RMSE_all), as well as the error over deleted (RMSE_mis) and observed regions (RMSE_vis), where the two numbers in parentheses correspond to the 10 % and 90 % percentiles of the error. Bold font indicates the best result for each metric and dataset.

4.3.1 Qualitative comparison

For further insights we visualize the CRITER and DINCAE2 reconstructions in Figs. 4 and 5. We showcase examples from the Mediterranean and the Adriatic test set, respectively, highlighting the masked SST (x_t⊙M_m), target SST (x_t), full reconstruction ($\stackrel{\mathrm{̃}}{\mathit{x}}$), standard deviation (σ), and RMSE computed over the entire target (RMSE_all). Notice that CRITER preserves fine details in cloud-free regions, ensuring minimal distortion of the original input data. In contrast, obscured (deleted) regions require the model to infer missing SST values using spatio-temporal context from adjacent days/pixels. These reconstructed regions exhibit reduced sharpness as a result of the inherent uncertainty caused by sparse observations. However, CRITER demonstrates an excellent ability to reconstruct high-frequency components of the target SST under deleted regions compared to DINCAE2. Additionally, CRITER proves robust to clouds of arbitrary shape, whether small and scattered (Fig. 4, first and last comparison) or large and contiguous (Fig. 4, second and third comparisons). Similar observations can be drawn from the comparisons on the Adriatic dataset presented in Fig. 5. On the Atlantic test set, both models face challenges in reconstructing high-frequency components under deleted regions, as illustrated in Fig. 6. However, we observe that CRITER is able to preserve the SST measurements over visible regions, whereas DINCAE2 introduces significant smoothing. Additional comparison figures are shown in Appendix B (Figs. B1, B2, and B3).

Figure 4 Comparison of sea surface temperature (SST) reconstructions generated by CRITER and DINCAE2 on the Mediterranean dataset. The columns display (1) the original SST field with simulated missing values, (2) the original SST field, (3, 4) full reconstruction of the SST field and the associated standard deviation, and (5) the absolute error map, highlighting the differences between the original and reconstructed fields. All panel values are in °C. Note that color scales for σ and RMSE_all are truncated at the 90th percentile of the data to improve visibility.

Figure 5 Same as Fig. 4 but for the Adriatic domain.

Figure 6 Same as Fig. 4 but for the Atlantic domain.

4.3.2 Spatial spectral analysis

We conduct spatial spectral analysis by comparing the power spectral density (PSD) of ground truth observations against reconstructions from CRITER and DINCAE2, focusing on the Ionian Sea region due to its significant SST variability.

First, we identify observation fields with maximum number of known measurements within the ROI (Region Of Interest) and compute their PSDs over the ROI. Following Fanelli et al. (2024), we compute PSD using a fast Fourier transform (FFT) with a Blackman–Harris window. We then sample 30 cloud masks with distinct coverage over the ROI, with the fraction of missing values ranging from 50 % to 98 %. For each mask, we simulate missing data in the observation fields, reconstruct them using both methods, and compute PSD over the reconstructed ROI.

Figure 7 shows an observation sequence with few available measurements. Both methods maintain PSD values near the target at low wavenumbers, indicating comparable low-frequency reconstruction. For wavenumbers k≥4 cycles per degree, however, CRITER's PSD remains closer to the target than DINCAE2's, demonstrating its superior ability to resolve high-frequency components. Figure 8 depicts a case with more measurements, where both methods generally align closer to the target. Nevertheless, CRITER still outperforms DINCAE2 at high wavenumbers ($k\ge \mathrm{5}\frac{\text{cycles}}{\mathrm{deg}}$). Additional results are provided in Appendix C.

Figure 7 Visualization of reconstruction performance. Row 1 shows the full fields (left to right: masked SST, target SST, CRITER reconstruction, and DINCAE2 reconstruction) with the Region of Interest (ROI) marked by a dashed black rectangle. Row 2 displays the corresponding ROI fields: target SST, CRITER reconstruction, and DINCAE2 reconstruction. Row 3 presents gradient magnitudes within the ROI for target, CRITER, and DINCAE2 outputs. Row 4 compares power spectral densities: target ROI (black), CRITER mean ± SD (orange band), and DINCAE2 mean ± SD (blue band), with solid orange and dotted blue lines showing CRITER's and DINCAE2's PSDs for the selected example.

Figure 8 Same as Fig. 7 but for another sample.

4.3.3 Comparison under different cloud coverage levels

The qualitative results presented in Sect. 4.3.1 suggest that CRITER is robust to clouds of various size. To test this, we compare the reconstruction error of CRITER and DINCAE2 on images with different coverage levels. The cloud coverage is given by the fraction of pixels that are missing or deleted relative to the total number of pixels belonging to the sea. Specifically, we categorize clouds into three distinct groups based on their coverage: low coverage (0 %, 60 %], moderate coverage (60 %, 75 %], and high coverage (75 %, 100 %). We then compute the reconstruction error within each group to assess the performance of both models under varying cloud conditions.

On the Mediterranean test set, the cloud coverage ranged from a minimum of 8.7 % to a maximum of 99 %. CRITER outperformed DINCAE2 across all cloud coverage groups, achieving significant reductions in reconstruction error over deleted regions. Specifically, the error was reduced by 21 % in the low-coverage group, 18 % in the moderate-coverage group, and 16 % in the high-coverage group. Similarly, on the Adriatic test set, the cloud coverage ranged from a minimum of 3.4 % to a maximum of 93 %. Here, CRITER substantially reduced the reconstruction error over deleted regions by 38 % in the low-coverage group, 49 % in the moderate-coverage group, and 54 % in the high-coverage group. Finally, on the Atlantic test set, the cloud coverage ranged from a minimum of 37 % to a maximum of 97 %. CRITER achieved a 4 % decrease in the low coverage group and around a 1.3 % decrease in moderate- and high-coverage groups.

Figure 9 Reconstruction error comparison between CRITER and DINCAE2 across different cloud coverage groups (low, moderate, and high) on the Mediterranean, Adriatic, and Atlantic test sets. The three rows correspond to the RMSE computed over (1) all ground truth measurements, (2) missing measurements, and (3) observed measurements. The error bars indicate the 10 % percentile, mean, and 90 % percentile of the error, respectively.

4.4 Uncertainty estimation and bias analysis

CRITER and DINCAE2 estimate both the reconstruction of missing values and the associated uncertainty (i.e., the standard deviation) for each pixel. To assess the reliability of the estimated standard deviation, we employ the scaled error metric

$$\begin{array}{}\text{(8)}& {\mathit{ϵ}}_{\left(i\right)}={\displaystyle \frac{{\mathit{x}}_{\left(i\right)}-{\stackrel{\mathrm{̃}}{\mathit{x}}}_{\left(i\right)}}{{\mathit{\sigma }}_{\left(i\right)}}},\end{array}$$

as proposed by Barth et al. (2020). This metric quantifies the difference between the ground truth observation x_(i) and the reconstruction ${\stackrel{\mathrm{̃}}{\mathit{x}}}_{\left(i\right)}$, normalized by the estimated standard deviation σ_(i), where i is the pixel index. We calculate the mean, μ_ϵ, and standard deviation, σ_ϵ, of the scaled error over the entire test set. Furthermore, we compute the bias, defined as the (non-normalized) mean difference between the ground truth observations and reconstructions. An ideal reconstruction method would thus have the bias equal to zero (i.e., predicted values are not globally under or over estimated). and standard deviation of the scaled error σ_ϵ equal to 1 (i.e., per-pixel disparities match the predicted uncertainties). Standard deviation of the scaled error σ_ϵ<1 indicates that the predicted standard deviation σ is overestimated, while σ_ϵ>1 indicates that σ is underestimated.

Figure 10 displays the histogram of the scaled error metric ϵ_(i) for each test set, along with the corresponding Gaussian distribution, characterized by the estimated mean μ_ϵ and standard deviation σ_ϵ. The mean (μ_ϵ), standard deviation (σ_ϵ), and the bias for each dataset are provided in Table 2. Notably, CRITER moderately underestimates the standard deviation, with standard deviation of the scaled error σ_ϵ values of 1.116, 1.082, and 1.156 on the Mediterranean, Adriatic, and Atlantic datasets, respectively, ranging from 8 % to 16 %. In contrast, on average, DINCAE2 significantly overestimates the standard deviation, with σ_ϵ values of $\mathrm{0.334},\mathrm{0.996},$ and 0.801 across the three datasets. The over-estimation thus ranges from as little as 0.4% to substantial over-estimates of 66 %. CRITER consistently exhibits a very low bias (of the order of 0.01 °C or lower) over all datasets. Furthermore, CRITER exhibits a significantly smaller bias on the Mediterranean and Adriatic datasets than DINCAE2, whereas DINCAE2 achieves a smaller bias on the Atlantic dataset. Note that, on the Adriatic dataset, DINCAE2 exhibits 18× larger bias than CRITER.

Table 2 Comparison of CRITER and DINCAE2 on each test set, showing the mean of the scaled error (μ_ϵ), standard deviation of the scaled error (σ_ϵ) – both unitless and bias in °C. Bold font indicates the best result for each metric and dataset.

Figure 10 Histograms of the scaled error ϵ_(i) for the Mediterranean, Adriatic, and Atlantic datasets, overlaid with the corresponding Gaussian distributions, which are characterized by the estimated mean (μ_ϵ) and standard deviation (σ_ϵ). Additionally, an ideal model is shown in black.

4.5 Ablation study

We analyze the proposed CRITER architecture by ablating or replacing individual parts. All model variants are trained for a total of 500 epochs (CRM and IRM are trained for 200 and 300 epochs, respectively) using the hyper-parameters described in the “Implementation details” section (Sect. 4.1). The variants are evaluated on the Mediterranean dataset.

4.5.1 Importance of the coarse reconstruction stage

We evaluate CRITER${}_{\stackrel{\mathrm{‾}}{\text{CRM}}}$, a variant without CRM. In this configuration, IRM is initialized with an uninformative prior ${\stackrel{\mathrm{̃}}{\mathit{x}}}^{\left(\mathrm{0}\right)}=\mathbf{0}$ and operates without FFM. Consequently, the first Residual Estimation Network (REN) assumes responsibility for the low-frequency reconstruction, a task previously performed by the transformer-based CRM. Table 3 demonstrates that incorporating CRM reduces the error over deleted and visible regions by 24 % and 87 %, respectively, validating the use of a transformer-based model for estimating the low frequency components.

Table 3 Performance of CRITER and CRITER${}_{\stackrel{\mathrm{‾}}{\text{CRM}}}$, which does not utilize CRM. Bold font indicates the best result for each metric.

4.5.2 Architectural design of CRM

Vision Transformer-based backbone

CRM (Sect. 3.1) utilizes a Vision Transformer-based architecture to compute a coarse reconstruction. The main argument for the transformer-based design is to allow direct information flow from all observed measurements into all corresponding tokens and the final coarse reconstruction. To evaluate the transformer design choice, we replace it by a convolutional counterpart, which maintains the same spatial reduction as the original CRM.

The convolutional variant, denoted with CRITER_CNN, utilizes a CNN-based CRM which accepts the same input as the original CRM and reduces the spatial resolution 8× by three double conv blocks, each followed by a 2×2 max pooling operation. The number of convolutional kernels in each block is 64,128, and 256, respectively. This is followed by a bottleneck layer, consisting of a double conv block with 512 convolutional kernels. The output latent features are upsampled to the original spatial resolution by applying 3 transpose convolution layers, each followed by a double conv block. The number of convolutional kernels in each block is 256,128, and 64, respectively. Finally, a single 1×1 convolutional layer computes the coarse reconstruction, which is passed to the IRM along with the latent features.

Results in Table 4 show that using CRITER_CNN leads to a substantial increase in reconstruction error over deleted and visible regions. This verifies the importance of the transformer-based design of CRM and suggests that global information flow plays an important role in obtaining good latent features and the coarse reconstruction.

Table 4 Performance of CRITER variants with different backbones. CRITER_CNN utilizes a CNN-based CRM, while CRITER utilizes the proposed ViT-based CRM. Bold font indicates the best result for each metric.

Modeling spatio-temporal data dependencies

We next inspect the importance of using spatio-temporal information in the coarse reconstruction. For this reason we remove the IRM from CRITER, leading to only using our proposed spatio-temporal masked-auto-encoder-based CRM architecture for reconstruction. We compare the reconstruction capabilities of CRM with the recent MAESSTRO (Goh et al., 2024), which also employs a Vision Transformer (ViT) (Dosovitskiy et al., 2021) and is based on a masked autoencoder (He et al., 2022). In fact, the major difference is that CRM utilizes three temporally consecutive SST fields to reconstruct the central field, while MAESSTRO uses only the central field.

Table 5 demonstrates that CRITER${}_{\stackrel{\mathrm{‾}}{\text{IRM}}}$ reduces reconstruction error by 44 % and 56 % over deleted and visible regions, respectively, compared to MAESSTRO. These results confirm the CRM modeling capability of spatio-temporal data dependencies, which considerably improves reconstruction performance.

Table 5 Performance of MAESSTRO and CRM. Bold font indicates the best result for each metric.

4.5.3 Importance of the refinement stage

To investigate the importance of refinement, we compare CRITER with two variants. The first variant, CRITER${}_{\stackrel{\mathrm{‾}}{\text{IRM}}}$, does not utilize refinement and takes the output of CRM as the final reconstruction. The second variant CRITER${}_{\stackrel{\mathrm{‾}}{\text{res}}}$ modifies IRM to estimate the full reconstruction at each iteration (in contrast to the proposed IRM that estimates a sequence of residuals).

Table 6 shows that utilizing refinement consistently leads to improved reconstruction. In particular the proposed IRM reduces CRM reconstruction error by 24 % and 91 % over deleted and visible regions, respectively. Furthermore, the results confirm that our proposed approach of consecutive residual estimation leads to lower errors than when the full signal is reconstructed at each refinement step. We hypothesize two reasons for this result. First, the residual estimation approach better exploits the individual REN networks in IRM, allowing each network to dedicate the full capacity for correction of the errors from the previous REN, thus gradually focusing on the high-frequency content reconstruction. Secondly, since the final reconstruction is obtained by summing the residuals, this enables a better gradient flow directly to each REN, thus enabling better training.

Table 6 Performance of CRITER variants using different refinement approaches. CRITER${}_{\stackrel{\mathrm{‾}}{\text{IRM}}}$ does not utilize refinement, and CRITER${}_{\stackrel{\mathrm{‾}}{\text{res}}}$ modifies IRM to estimate the full reconstruction at each iteration, while CRITER utilizes the proposed IRM. Bold font indicates the best result for each metric.

4.5.4 Influence of refinement iteration steps

We next investigate the impact of varying the number of refinement steps in IRM (Sect. 3.2) on the reconstruction quality. Figure 11 shows results of CRITER retrained with different number of steps in IRM. The lowest reconstruction error is reached at three refinement steps. In particular the RMSE_all is reduced by 8 % compared to using a single refinement step. Using more refinement steps does not improve performance but leads to increased error. This is likely due to the parameter increase, since each refinement step introduces a new REN network, which makes training less efficient on the limited dataset size. We defer explorations of more resilient IRM architectures to future work.

Figure 11 Performance of CRITER variants with increasing number of refinement iterations.

4.5.5 Importance of the CRM latent features

In IRM (Sect. 3.2), the latent features computed by CRM are fused with the bottleneck features to improve injection of global coarse information in the refinement steps. To evaluate the importance of this, we retrained CRITER without the coarse latent features fusion in IRMs – this variant is denoted as CRITER${}_{\stackrel{\mathrm{‾}}{\text{fus}}}$.

Results in Table 7 show that the reconstruction error over deleted and visible regions of CRITER with feature fusion reduces by 1.9 % and 19 %, respectively, compared to the feature fusion free counterpart.

Table 7 Comparison of CRITER, which utilizes latent features computed by CRM, with CRITER${}_{\stackrel{\mathrm{‾}}{\text{fus}}}$, which does not. Bold font indicates the best result for each metric.

4.5.6 Importance of time auxiliary features

CRM (Sect. 3.1) takes as input a sequence of consecutive observation fields, that are concatenated with auxiliary features, particularly the cosine and sine of the day of the year that encode the yearly cycle of SST. The auxiliary features offer additional information which CRM can incorporate when computing the latent features and generating the coarse reconstruction. To evaluate the importance of this, we train a CRM variant which does not leverage auxiliary features, denoted by CRITER${}_{\stackrel{\mathrm{‾}}{\text{aux}}}$. Results in Table 8 show that augmenting the input with auxiliary features leads to a 1.5 % and 26 % decrease in reconstruction error over deleted and visible regions, respectively.

Table 8 Comparison of CRITER, which utilizes auxiliary features (cosine and sine of the day of the year), with CRITER${}_{\stackrel{\mathrm{‾}}{\text{aux}}}$, which does not. Bold font indicates the best result for each metric.

5 Conclusions

This study introduced CRITER, a novel two-stage model for reconstructing sea surface temperature (SST) from sparse satellite observations. High performance of the CRITER method stems from a Coarse Reconstruction Module (CRM) utilizing a Vision Transformer (ViT) architecture for initial reconstruction, followed by an Iterative Refinement Module (IRM) to refine the reconstruction with a focus on high-frequency information. The global receptive field of the ViT enables modeling of long-range dependencies in the data, while iterative refinement allows each network to focus its full capacity on modeling high-frequency corrections. This combination leads to significant enhancements in overall performance. The introduction of CRM's ViT global attention mechanism proved crucial for effective long-range dependency modeling, addressing limitations of convolutional architectures.

Our results show that CRITER surpasses the state-of-the-art DINCAE2 model by a significant margin across three diverse SST datasets: Mediterranean, Adriatic, and Atlantic. Notably, CRITER achieves substantial reductions in reconstruction error, with improvements of up to 89 % in observed regions and up to 44 % in missing regions.

The iterative refinement process of IRM, focusing on residual estimation, further enhanced reconstruction accuracy by efficiently utilizing model capacity for high-frequency variability in the SST observations. Ablation studies confirmed the importance of CRM's transformer-based design, the effectiveness of iterative residual estimation in IRM, and the utility of incorporating auxiliary features such as the day-of-year encoding.

Overall, CRITER sets a new benchmark for SST reconstruction, providing a robust framework that leverages the strengths of both transformer and convolutional architectures to deliver superior performance. Future work will explore extending CRITER's applicability by incorporating additional environmental proxy variables (like chlorophyll a, which often serves as a complementary variable to SST in ocean state estimates) and increasing the temporal horizon for even more accurate sparse data reconstructions.

Appendix A: Analysis of missing values in evaluation datasets

We analyze the extent of missing values in each dataset described in Sect. 2.1. To quantify the number of missing data, we define the cloud coverage A_t of an observation x_t as

$$\begin{array}{}\text{(A1)}& A_t={\displaystyle \frac{|\mathbf{1}-{\mathbf{M}}_{\mathrm{t}}|}{\left|{\mathbf{M}}_{\mathrm{l}}\right|}},\end{array}$$

where ${\mathbf{M}}_{\mathrm{t}}\in \mathit{\{}\mathrm{0},\mathrm{1}{\mathit{\}}}^{W\times H}$ is the missing data mask corresponding to observation x_t, and ${\mathbf{M}}_{\mathrm{l}}\in \mathit{\{}\mathrm{0},\mathrm{1}{\mathit{\}}}^{W\times H}$ is the land mask. Cloud coverage is computed as the fraction of pixels that are missing relative to the number of pixels belonging to sea areas. We then calculate the mean cloud coverage over Δ_t=3 consecutive observation fields as $A=\frac{\mathrm{1}}{\mathrm{3}}(A_{t-\mathrm{1}}+A_t+A_{t+\mathrm{1}})$. Note that the proportion of available information in the entire observation triplet is thus given by 1−A. Figure A1 presents a histogram of the mean cloud coverage A for all three filtered datasets. The results show that the datasets can be ranked by the average amount of available information in each observation triplet, from highest to lowest: Mediterranean, Adriatic, and Atlantic.

Figure A1 Histogram (100 bins) of cloud coverage for all three (filtered) datasets.

Appendix B: Additional qualitative analysis figures

This section presents additional reconstructions generated by CRITER and DINCAE2 (Barth et al., 2022). For a detailed discussion of the qualitative comparison, refer to Sect. 4.3.1.

Figure B1 Same as Fig. 4 on different samples.

Figure B2 Same as Fig. 4 but for the Adriatic domain.

Figure B3 Same as Fig. 4 but for the Atlantic domain.

Appendix C: Extended spatial spectral analysis

This section presents supplementary power spectral density (PSD) comparisons. Figure C1 shows a challenging case with sparse measurements, where CRITER's PSD remains closer to the target (on average) for wavenumbers $k\ge \mathrm{4}\frac{\text{cycles}}{\mathrm{deg}}$. Figure C2 depicts a high-measurement scenario featuring a failure case for CRITER: minor noise amplification beyond $k\ge \mathrm{5}\frac{\text{cycles}}{\mathrm{deg}}$. A similar issue occurs with DINCAE2 but in a different wavenumber band: Fig. C3 shows significant noise amplification within $k\in [\mathrm{2},\mathrm{4}]\frac{\text{cycles}}{\mathrm{deg}}$. For a detailed discussion of the comparison, refer to Sect. 4.3.2.

Figure C1 Same as Fig. 7 but for a different sample.

Figure C2 Same as Fig. 7 but for a different sample.

Figure C3 Same as Fig. 7 but for a different sample.

Appendix D: Implementation details of baseline models

MAESSTRO (Goh et al., 2024) is trained using mean squared error loss, as described in Sect. 3.3, for consistency in our comparison. The model processes only the current time step SST field without auxiliary features. We modified MAESSTRO's original random patch masking to inject sampled real cloud masks, enhancing real-world applicability. An SST patch is masked if its corresponding cloud mask patch contains any zero values. MAESSTRO employs a ViT-Tiny backbone with 12 encoder and decoder layers, 3 multi-head attention (MHA) heads, a token dimension of D_t=192, layer-norm epsilon of $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{12}}$, and patch size of 8×8. MAESSTRO is trained with a batch size of 8 using the AdamW optimizer with a learning rate $\mathit{\alpha }=\mathrm{3}\times {\mathrm{10}}^{-\mathrm{4}}$, β₁=0.9, and β₂=0.95 for 100 epochs (warm-up period) and then with a cosine decay scheduler (Loshchilov and Hutter, 2016) with a step size of 50 for another 300 epochs.

DINCAE2 (Barth et al., 2022) is trained using the negative log-likelihood loss, as described in Sect. 3.3, to maintain consistency in our comparison. The model utilizes a sequence of three temporally consecutive SST fields, along with day-of-the-year auxiliary features, to reconstruct the central SST field. Hyperparameters of the re-implemented DINCAE2 differ slightly between the datasets. On the Mediterranean and Atlantic, DINCAE2 is trained using the Adam optimizer, with an initial learning rate of $\mathit{\alpha }=\mathrm{4}\times {\mathrm{10}}^{-\mathrm{3}}$, β₁=0.90, and β₂=0.999 and a batch size of 8 for a total of 1000 epochs, using a step learning rate scheduler with a step size of 100 epochs and a multiplicative factor of γ=0.5. On the Adriatic we use an initial learning rate of $\mathit{\alpha }=\mathrm{7}\times {\mathrm{10}}^{-\mathrm{3}}$ and a step size of 150; all other hyperparameters remain unchanged.