The observations that underpin the system presented here include along-track satellite altimetry from the Radar Altimeter Database System (RADS, Scharroo et al., 2013). The along-track data is corrected for various geophysical and instrumental effects (Scharroo, 2018), including the effects of atmospheric pressure. Additionally, the sea-level observations are converted to SLA by removing an appropriate mean sea surface height. The experiments presented here span 2023, when altimeter data is available from six altimeter missions, including Sentinel-3A, -3B, and -6A, SARAL/Altika, Jason-3, and Cryosat-2. The altimeter data analysed in this study have a sampling frequency of 1 Hz, which translates to a spatial resolution of roughly 7 km along-track. The repeat cycle for Sentinel-3A and -3B is 27 days; for Sentinel-6A and Jason-3 is 10 days; and for SARAL/AltiKa is 35 days; while CryoSat-2 has no exact repeat cycle.

The gridded SLA fields used as an archive for our system are from many different sources. The archive of states from which analogs are identified is sometimes called a catalogue (e.g., Platzer and Chapron, 2024) or a library (e.g., Ding et al., 2018). This includes many different models, namely version 3 and version 4 of the Ocean Forecasting Australia Model (OFAM, Oke et al., 2013), version 3 of the LASG/IAP Climate System Ocean Model (LICOM3, Ding et al., 2022), the 0.1^◦-resolution version of the Australian Community Climate and Earth System Simulator Ocean Model (ACCESS-OM2-01, Kiss et al., 2020), and the Ocean general circulation mode For Earth Simulations (OFES, Masumoto et al., 2004). Relevant details of the specific model runs are summarised in Table 1. All of the models used here are based on either version 3 or 4 of the Modular Ocean Model (MOM, Griffies, 2009). Model runs are forced with atmospheric fields from JRA55 (Kobayashi et al., 2015), ERA-Interim (Dee and Uppala, 2009), ERA-20C (Poli et al., 2016), NCEP (Kalnay and Kanamitsu, 1996), or a 17-member ensemble-mean CMIP5-RCP8.5 (e.g., Kharin et al., 2013). The horizontal resolution for all models is 0.1^◦, although the data from C100 (using LICOM3.0) are only freely available at 0.25^◦-resolution. Importantly, all model runs provide daily-averaged SLA fields. Some models report sea-level that includes the impact of atmospheric pressure. For these elements of the archive, an inverse barometer correction has been applied.

The model runs included in Table 1 are a mix of standard model integrations, with realistic interannual forcing (IAF) and some are model simulations with repeat year forcing (RYF). Of the RYF runs, Y79R is a 123-year model run with repeat-year forcing from 1979 (Zhang et al., 2016); AOMR is a 230-year model run with repeat-year forcing from 1990/1991 (Stewart et al., 2020); and OFSC is an 8-year model run with climatological forcing fields from an NCEP reanalysis (Masumoto et al., 2004). The runs with ACCESS-OM2-01 (AOM1, AOM2, AOM3, and AOM4) are IAF runs, but include four cycles, each spanning 61 years, with forcing from JRA-55 (Kiss et al., 2020).

In addition to data from model runs, for some experiments, we also use data from ocean reanalysis (RA) and ocean analysis (OA) products. The RA data sets include two versions of the Bluelink ReANalysis (BRAN) - version 2020 (BR20; Chamberlain et al., 2021b, a) and version 2023 (BR23; which includes incremental improvements compared to BRAN2020). The OA products are both global gridded estimates of SLA, produced using optimal interpolation. These include IMOS Ocean Current (OCUR; http://oceancurrent.imos.org.au) and Ssalto/DUACS (SSDU), produced under Copernicus Marine Environment Monitoring Service (CMEMS; https://marine.copernicus.eu).OceanCurrent is a 0.2^◦-resolution product, and Ssalto/DUACS is a 0.125^◦-resolution product. The RA and OA products all span 2023 - which is the test period for this study. When part of the archive, gridded fields after September 2022 are excluded from the archive to ensure a “fair” assessment of the analog forecasts.

For this study, we assemble a global archive of daily SLA fields. This permits efficient application to any region of the world. However, to assess the system, we produce forecasts of SLA for 12 regional domains around Australia (Figure 1). Each domain spans 5x5^◦ – approximately 500x500 km at 30^◦S – and includes a range of dynamic regimes, including western and eastern boundary currents, eddy-rich fields, and wind-driven regions in both tropical and mid-latitude zones.

The most computationally intensive step in our system is reading the archived data into memory. This is made more efficient by “cutting out” regional subsets of each model in advance for each domain.

As noted in the introduction, most NWP centres have abandoned analog forecasting, citing the difficulty finding “good analogs”. Here, we use data from many different sources, hoping to expand the catalogue of states from which analogs can be identified. The full aggregated archive of states is assembled from 16 different model runs or analysis products (Table 1) from over 1000 years of data. To evaluate how data from additional sources contribute independent information to the archive, we perform a singular value decomposition (SVD) on various subsets of the archived fields. We analyse data from Region 1, which we assume to be representative of other regions. We perform the SVD for each model run or analysis product individually; and we also calculate the SVD for aggregated archives, starting with a single model run (m=1) and progressively including data from additional model runs and analysis products, up to m=16.

Since the square of each singular value of a data matrix is proportional to the variance explained by its corresponding mode, the singular value spectrum shows how the effective dimension of the state space changes with the addition of each new source of SLA data. To prevent artificial inflation of the variance due to increasing sample size, we normalise each data matrix by 1/n, where n is the number of daily fields in the archive. This normalisation ensures that increases in singular values truly reflect new patterns of variability, not simply more data. When a new dataset contributes no independent information, its addition doesn’t change the singular values; and the data from the new source merely projects onto existing modes. By contrast, an increase in singular values indicates an increase in the effective dimension of the subspace spanned by the archive.

Results are presented in Figure 2, showing the spectra of singular values for: (a) all OFAM runs, (b) both ocean analyses, (c) all ACCESS-OM2–01 runs, and (d) all model runs and analyses listed in Table 1. In Figure 2a, although six OFAM runs are included, they fall into three distinct groups with independent variability: the ERA-Interim–forced run (EISP), the JRA-55–forced runs (JRSP, Y79R), and the reanalyses (BR20, BR23). We find that CM85 adds no new dimensions on top of EISP, JRSP, and Y79R. We also see that both reanalyses (BR20 and BR23) project onto equivalent modes, together increasing the overall dimension of the archived dataset.

The two observation-based ocean analyses (OCUR and SSDU; Figure 2b) contribute distinct modes, likely due to differences in gridding methods and spatial resolution. All of the ACCESS-OM2–01 runs (Figure 2c) project onto the same subspace; and data from the OFES model introduces new, independent variability (Figure 2d), despite the shorter duration of those runs. We find that the C100 dataset, generated using LICOM3.0, also contributes unique modes, expanding the effective dimension of the archive.

Based on analysis of singular values (Figure 2), we conclude that multiple runs from the same model configuration generally span the same subspace. Consequently, adding more runs from the same configuration does not increase the effective dimensionality of the archive. To further explore the contribution of information to the aggregated dataset, we also examine the associated singular vectors.

We find that the right singular vectors, representing the dominant spatial modes, are effectively equivalent across all runs from the same model configuration (though the exact order of modes and the details are not precisely the same). This confirms that these runs share the same underlying spatial structures - consistent with the analysis of singular values. By contrast, the left singular vectors, representing the temporal variability, are uncorrelated across different runs - even for the repeat cycles of AOM* that span the same time periods and forced with the same fluxes. This indicates that while each run projects onto the same set of spatial modes, the combination and weighting of those modes varies in time for different runs.

The implication of this result is important: although the spatial subspace is unchanged, each new run contributes a unique combination of modes. As a result, the full set of archived fields becomes richer in its diversity of states. For analog forecasting, this means that adding more runs - even from the same model - broadens the catalogue of historical states that can be searched for analogs, increasing the likelihood of identifying a close match to the observed state.

This finding is consistent with a well-known property of high-dimensional spaces - that independently drawn random vectors from a high-dimensional state are likely to be nearly orthogonal (Christiansen, 2018). It is also supported by our sensitivity experiments that show that forecast skill consistently improves as more runs are included in the archive.

To benchmark the skill of our analog forecasting system, we compare its predictions against archived forecasts from OceanMAPS version 4.0i (OMv4.0i; Brassington et al., 2023), accessed as they were available at the time of each forecast (last retrieved 14 February 2025). OMv4.0i is a traditional ocean forecast system that uses a hybrid-EnKF data assimilation system to initialise a near-global ocean general circulation model (OFAM). OMv4.0i uses a 48-member dynamic ensemble and a 144-member stationary ensemble, similar to Counillon et al. (2009). This approach to data assimilation combines the Ensemble Kalman Filter (Evensen, 2003) and Ensemble Optimal Interpolation (Oke et al., 2002b; Evensen, 2003; Oke et al., 2010) for estimating the system’s background error covariance. OMv4.0i produces forecasts that are of comparable performance to other operational systems run around the world (e.g., Aijaz et al., 2023) and are significantly better than version 3 of the same system that preceded it (Brassington et al., 2023).

To independently assess each forecast produced here, we compare forecast SLA fields to a verifying analyses from IMOS OceanCurrent (http://oceancurrent.imos.org.au, last accessed on 14 February 2025). SLA analyses from OceanCurrent are available on a 0.2^◦-resolution grid, with analyses every day. OceanCurrent analyses are produced by merging along-track SLA from all available satellite altimeter missions together with SLA observations from coastal tide gauges around Australia using optimal interpolation.

To evaluate our system, we conduct 25 independent 15-day forecasts, initialised every 15 days throughout 2023. These forecasts span 12 regional domains around Australia (Figure 1), encompassing a range of dynamic regimes. Most results presented here are from Experiment 1 (Exp. 1), which uses sea level anomaly (SLA) fields from the free-running models listed in Table 1 - EISP, JRSP, Y79R, C100, AOM*, and OFS* - comprising a total of 902 years of daily SLA data.

Figures 6, 7 summarise the system performance across all domains. Figure 6 compares MAD between SLA in the verifying analysis and forecasts from three sources, including the kNN-based equal-weighted ensemble mean (kNN EqW); OceanMAPS version 4.0i (OMv4.0i); and the first nearest neighbour (NN1). Figure 7 shows an equivalent analysis, but using ACC to evaluate forecasts.

For Figures 6, 7, we show the mean metrics computed from 25 forecasts. The shaded bands represent the standard error of the mean across those forecasts, and the associated p-values are also reported in each panel. To identify the “best” forecast performance for each region, we compare the metrics (averaged over forecast days 1-6 - limited to 6-day forecasts because that is available from OMv4.0i) between three forecast types: kNN EqW, OMv4.0i, and NN1.

For each pairwise comparison, the null hypothesis is that the two forecast methods have equal mean performance. We evaluate this hypothesis using a standard two-sample t-test (e.g., Browne, 2010) applied to the distribution of metrics across the 25 forecasts. A p-value (e.g., Vidgen and Yasseri, 2016) below 0.05 indicates that the observed difference in mean metrics would occur by chance less than 5% of the time under the null hypothesis. When the p-value is below this threshold, we regard the difference between forecast methods as statistically significant.

For regions where one forecast type significantly outperforms the benchmark – or where the benchmark significantly outperforms the kNN-based forecasts – we denote the region with a coloured box in Figures 6, 7.

Using MAD as a metric, kNN EqW outperforms OMv4.0i in 4 of 12 domains, NN1 outperforms OMv4.0i in 1 of 12 domains, and OMv4.0i outperforms the kNN-based forecasts for 3 of 12 domains (Figure 6). For 3 of 12 domains, the mean MAD are not significantly different.

Using ACC, kNN EqW outperforms OMv4.0i in 6 of 12 domains; NN1 outperforms OMv4.0i in 1 of 12 domains, and OMv4.0i outperforms the kNN-based forecasts in 1 of 12 domains (Figure 7). For 4 of 12 domains, the mean ACC are not significantly different.

Based on these results, we conclude that the kNN/analog system outperforms traditional operational forecasts in 40-60% of cases, performs equally well (no statistical difference) in about 30% of cases, and is outperformed in about 10-25% of cases.

A distinct feature of OMv4.0i forecasts is a sharp increase in MAD and drop in ACC between days 0 and 1 (Figures 6, 7). We attribute this to initialisation issues that are caused by the dynamically inconsistent initial conditions, as noted in the introduction. By contrast, analog forecasts avoid such degradation since each forecast is from a free-running model, and is therefore dynamically consistent.

It is interesting to note the relative number of NNs that are identified from each source in the assembled archive (Table 2). In total, we identify 3600 NNs for this study, comprised of an ensemble of 12 for 25 forecasts over 12 domains. We report results from two configurations: Exp. 1, where only dynamically consistent SLA fields are used (from EISP, JRSP, Y79R, CM85, C100, AOM*, and OFS*); and Exp. 17, using all data in the archive (the same as Exp. 1, plus BR20, BR23, OCUR, and SSDU). For Exp. 1, we find that there is a disproportionally high number of NNs from AOM* IAF runs – accounting for 38.1% of NNs in Exp. 1, having contributed 30.3% of the archive. For Exp.17, when we include ocean reanalyses and observation-based analyses in the archive, this weighting towards the AOM* IAF runs reduces, and the observation-based analyses contribute a disproportionally large percentage of NNs. Together, OCUR and SSDU account for 11.5% of NNs, having only comprised 5.9% of the archive; and fields from the BRAN experiments account for 7.4% of NNs, having contributed only 4.4% to the archive.

To examine the mean characteristics of the ensemble forecasts, we present mean rank histograms (Hamill, 2001) for each region – averaged over forecast days 0–15 for 25 forecasts – in Figure 8. The histograms reveal clear and spatially systematic differences in ensemble performance across the 12 regions. Several regions exhibit a pronounced U-shape (e.g., regions 3, 4, 8, and 9), indicating that the ensemble spread is generally under-dispersive and that observations frequently fall outside the ensemble range (e.g., Hamill, 2001). Region 2 displays a dome-shaped distribution, suggesting an over-dispersive ensemble. Other regions – particularly regions 6, 7, and 10 – show flatter, more uniform histograms that imply a more appropriate spread. Regions with strong skewness towards the lowest or highest ranks (e.g., regions 3, 4, 8, and 12) indicate localised biases in the ensemble mean.

Overall, the rank histograms indicate that the analog-based ensembles capture important aspects of forecast uncertainty, although the spread is insufficient in a number of regions and occasionally exhibits systematic skewness. These deficiencies likely reflect limitations in the diversity or representativeness of the analog archive, rather than issues associated with ensemble data assimilation. The poorest-performing regions (3, 4, 8, 9, and 12) show clear U-shaped histograms, consistent with under-dispersion. In three of these regions (3, 4, and 9), the analog ensemble usually underperforms relative to OMv4.0i (Figures 6, 7). This suggests that kNN performance could be improved by expanding or diversifying the archive so that the selected states span a broader range of independent dynamical conditions.

To further illustrate forecast performance and practical application, we present a number of detailed examples (Figures 9–13). Each example shows the verifying analyses, equal-weighted NN ensemble means, and forecasts from 12 individual NNs at t+0, t+5, t+10, and t+15. Results described in this section are from Exp. 1, using only SLA fields from free-running models. The examples discussed below include a mix of cases when the analog forecast system performs well, and some where it performs poorly. To further demonstrate the performance, we also offer results from 35 other examples (2 or 3 for each region) in the Supplementary Material.

Interestingly, for most cases, the initial conditions of each NN (at t+0) are often quite different from the verifying analyses at t+0. This is not surprising, because each NN is identified based on their similarity to observed evolution over a 10-day period, not necessarily their initial state (at t+0).

A key question in ensemble data assimilation studies relates to the optimal ensemble size. For this study, this is the same as seeking the optimal number of NNs for each forecast. The efficiency of the kNN approach readily permits performance of a comprehensive set of sensitivity tests.

For each forecast, we calculate the MAD between the ensemble mean of the first 1, 2···100 NNs and the verifying analyses (from OceanCurrent). Figure 14 shows the average MAD as a function of ensemble size, averaged over all 300 forecasts (25 forecasts across 12 domains). While no single optimal ensemble size is evident, averaging 10–25 NNs seems to produces the most accurate ensemble-mean forecast. This result guided our choice of 12 NNs for results presented in this study.

Compared to traditional ocean forecast systems, there are relatively few “tuneable parameters” in an analog forecast system. However, results can be sensitive to the data included in the archive and the method used to identify NNs. We summarise a series of sensitivity experiments in Table 3. This includes experiments using different sources of SLA in the archive, using a different metric for identifying NNs, using a different observation window for identification of the NNs, and using a different ensemble size. We present these mean statistics to demonstrate the sensitivity, but we have not applied the same rigorous statistical tests that we applied the results in Figures 6, 7. Here, we aim merely to quantify the sensitivity, rather than definitively identify the “optimal” configuration.

Several general conclusions can be drawn from these sensitivity experiments. The observation window used for identifying NNs influences the forecast skill. Comparing experiments with an observation window of 4d (Exp. 3), 7d (Exp. 5), 10d (Exp. 1), and 15d (Exp. 4), suggests that 15d is too long, and 4d is too short. Further, we conclude that a 7d window produces the best forecasts. Inclusion of SLA data from additional, or fewer, sources in the archive also influences the forecast skill. We also find that there is some improvement when a larger ensemble is used – though a more comprehensive assessment is warranted to strongly support this claim. Compare results with 12 members (Exp. 17), 6 members (Exp. 21), and 18 members (Exp. 22). For most cases, the results are intuitive, with inclusion of more data improving the forecast skill, and using less data degrading the forecast skill.

Overall, we find some sensitivities to different choices of configuration of the analog forecast system. Different domains seem to warrant different configurations. Noting the computational efficiency of this system, it’s possible that every individual domain could be tuned for the best performance. However, we also note that the range of ACC for experiments presented in Table 3 is small, so the improvement in forecast skill gained by tuning is not great. In fact, most of the differences reported here are not statistically significant, again highlighting the insensitivity of the kNN/analog approach to these choices.

Although kNN or analog forecasting is not yet widely used in oceanography, this study demonstrates its significant potential for mesoscale ocean forecasting. Compared to traditional ocean forecast systems, analog forecasting offers several advantages.

First, analog forecasting can produce dynamically consistent forecasts, avoiding the rapid degradation often seen with sequential data assimilation. This may be particularly beneficial for scenarios that are sensitive to dynamical imbalance, including forecasts of the mesoscale ocean circulation (e.g., Pilo et al., 2018) and biogeochemical (BGC) forecasting (e.g., Raghukumar et al., 2015). Analog forecasting can completely mitigate these issues by ensuring dynamic consistency in the forecast fields.

Second, analog forecasting allows targeted predictions for a subset of ocean variables (e.g., SLA, SST, or surface phytoplankton) without requiring a full ocean state to be predicted for every forecast. In this study, we chose to forecast only SLA, making our system computational efficiency, but still relevant to a number of end users. By contrast, traditional numerical models have to compute the entire ocean state for each forecast, even if the intent of the system is to forecast a subset of variables (e.g., SLA in one region).

Third, our approach is efficient and cost-effective. We generated an ensemble of 15-day forecasts for a 5^◦ × 5^◦ domain in about 90 seconds using a modest virtual machine (with 4 Intel Xeon Gold 6430 cores and 14 GB RAM under VMware virtualization). By contrast, OceanMAPS4.0i, the operational forecast system used for benchmarking, requires approximately 7,000 CPU hours (or 420,000 CPU minutes) to generate a global 6-day forecast, as reported in the 2023 OceanPredict National Report. The results in this study, for example, used 25 global OceanMAPS forecasts that were produced on a high-performance super-computer and required almost 20 CPU years to complete; and 300 regional analog forecasts (12 domains by 25 forecasts), required about 7.5 CPU hours. The computational cost of our analog approach is therefore about 25,000 times less than a traditional ocean forecasting system for this study.

Finally, the analog method can function as a multi-system tool. The model archive can incorporate fields from different models, configurations, and resolutions (Table 1). Since all models exhibit some systematic biases (e.g., excessive mixing, misplaced boundary currents, deep mixed layers), including multiple models allows the system to preferentially select NNs that minimise systematic errors.

Despite its advantages, analog forecasting has several limitations. Every analog forecast is merely selected from the model archive. The system can only predict events that have analogs in the archive, meaning that extreme events (e.g., marine heatwaves) can only be forecasted if similar cases exist in the archive of model states. This limitation becomes particularly relevant under climate change, where unprecedented conditions may emerge. One potential mitigation is to incorporate model simulations forced by future climate projections (e.g., CM85, C100; Table 1).

Another limitation is the dependence on observation quality and availability. If observations are sparse or lack adequate spatiotemporal coverage, selected NNs may be suboptimal, degrading forecast skill. To mitigate this, we used a 10-day sequence of observations for NN selection. Incorporating additional observation types (e.g., satellite SST or Argo profiles) is another avenue for improvement, as is combining analog forecasting with other model-data synthesis techniques (e.g., Rykova, 2023).

Moreover, analog forecasting is best suited for regional domains. Applying it at global or even basin scale is impractical due to the high dimensionality of ocean states. kNN methods are known to become unstable in high-dimensional spaces (e.g., Pestov, 2013). For this study, we performed experiments to systematically test domain sizes (from 3^◦ × 3^◦ to 10^◦ × 10^◦) and found that smaller domains generally yield more accurate forecasts. A 5^◦ × 5^◦ domain was chosen to balance accuracy and usability (i.e., we expect that a forecast on 500x500 km domain will be useful for many end users).

Global applications of analog forecasting would seem out of reach. Such an approach would require mesoscale features (e.g., eddies) across multiple ocean regimes to align, which is unlikely. Lorenz (1969) faced a similar issue when searching for analogs in atmospheric forecasting that included the entire northern hemisphere, noting that “there are numerous mediocre analogs but no truly good ones”. While global applications seem unrealistic, regional implementations show good potential.

The archive of SLA fields used here (Table 1) are mostly from models, and mostly include two types of model runs - simulations with IAF and climatological, or RYF. Our results (Table 2) suggest that interannual forcing yields better performance for analog forecasting, though fields from climatological runs remain valuable given the chaotic nature of mesoscale circulation.

Currently, all models in our archive are based on the MOM. Expanding the archive to include simulations from other models such as NEMO (Madec, 2016), the MIT GCM (Adcroft et al., 2004), or HYCOM (Chassignet et al., 2007) could further improve forecast skill.

Although most of the results presented in this study use SLA fields from free-running models, with no data assimilation, we have also performed forecasts using fields from ocean reanalyses and observation-based analyses (Tables 1-3). When these additional sources are included, we find that many NNs are identified from these datasets (Table 2) and the forecast skill typically improved (Table 3). Close inspection of the ensuring forecasts from these sources (not shown) expose the problems with these datasets - often with noticeable discontinuities in SLA over the forecast period. The tradeoff between improved forecast skill and dynamical inconsistent forecasts is something that needs to be considered for any application.

In our experiments, we chose to use ACC between along-track altimetry and model SLA fields to select NNs. We also ran a comprehensive set of experiments using MAD, which yielded similar selections. However, forecasts using ACC to select NNs resulted in better agreement with verifying analyses.

For SLA forecasting, we think ACC is the most appropriate metric because the primary goal of the system is to provide qualitative guidance. In this context, accurately predicting the location, size, and shape of mesoscale features is more important than precisely matching their amplitude, or the amplitude of SLA over a broader area.

The most commonly used error metrics are MAD, Mean Squared Difference (MSD), and ACC. These metrics are mathematically related. For a Gaussian sample, MSD ≈ π/2 MAD². Furthermore, MSD can be formally expressed in terms of ACC. The MSD (Equation 1) between two vectors X and Y can be expressed as:

where MSD is the mean squared difference, X−Y¯ is the bias, and σX and σY are standard deviations of X and Y (e.g., Murphy, 1996; Oke et al., 2002a).

These relationships indicate that selecting NNs based on MAD (or MSD) implicitly incorporates comparisons of mean fields (close to zero in this case, since SLA is an anomaly); the amplitude of mesoscale features (quantified by standard deviation); and the location and shape of mesoscale features (quantified by the ACC). Since our objective is to produce forecasts that qualitatively inform decision making, we prioritise ACC for selecting NNs, as it better captures the spatial structure of mesoscale variability.

While kNN and analog forecasting are not widely used in ocean forecasting, the approach shares similarities with ensemble data assimilation methods. It can be shown that a field of increments from EnOI- or EnKF-based systems can be expressed as a linear combination of ensemble anomalies (e.g., Oke et al., 2007, 2021). Moreover, one of the key advantages of EnKF over EnOI is that the EnKF generates ensemble members that are relevant to the current forecast (e.g., Sakov and Oke, 2008), sometimes called the “errors of the day” (e.g., O’Kane et al., 2011). An analog system also produces an ensemble mean of NNs (ensemble members) that are relevant to the current state. Here, we combine these NNs to produce an ensemble mean – a linear combination of members. We can therefore see that in this way, this approach is similar to EnOI and EnKF, but at a fraction of the computational cost. However, unlike EnOI or the EnKF, analog forecasting does not exploit this ensemble to assimilate data. An ensemble of analog forecasts are merely intended to produce multiple estimates of possible future conditions.

When ensembles are used by traditional forecast systems, the ensemble is used both as an ensemble prediction, and to quantify the system’s error covariance, for data assimilation. When used for data assimilation, the system generally performs better when a large ensemble is used (e.g., Mitchell et al., 2002). While this may also be useful for ensemble prediction, the quantity of forecasts may become unmanageable for users to exploit.

In ensemble forecasting, the ensemble mean is often considered the primary output (e.g., Counillon et al., 2009; Sakov and Sandery, 2015; Aijaz et al., 2023). However, as shown in Figures 9-13, individual ensemble members can provide valuable insight, particularly for mesoscale feature prediction.

End users may benefit from considering both the ensemble mean and individual forecasts. A key question is whether the ensemble mean is always the best metric for assessing forecast skill. The goal of an ensemble forecast is to represent a range of plausible future states. Even if the ensemble mean is inaccurate, the system succeeds if at least one ensemble member closely matches reality. This perspective is common in atmospheric forecasting and could be more widely adopted in oceanography.

For many of the examples presented here (Figures 9-13 and in the Supplementary Material), the most skilful forecast is often one of the NNs, not the ensemble mean of NNs. This may be true for traditional ensemble data assimilation systems, and is something that warrants consideration.

In this study, we demonstrate that analog forecasts can outperform traditional forecasts for a subset of variables. Applications to NWP drew similar conclusions (e.g., Langmack et al., 2012; Nagarajan et al., 2015). Here, we consider the potential value of forecasting only a subset of the ocean state, compared to the entire state. Many marine industries rely on ocean forecasts to guide decision-making (e.g., Kourafalou et al., 2015; Schiller et al., 2020; Rautenbach and Blair, 2021; Ciliberti et al., 2023; Spillman et al., 2025). Most forecasts provided by national weather centres offer a comprehensive suite of data and products (e.g., Brassington et al., 2023). However, many end users do not fully exploit these forecasts. Instead, they tend to make decisions based on qualitative assessments – simply analysing images of predicted conditions (e.g., Höllt et al., 2013; Kuonen et al., 2019; Spillman et al., 2025) or considering indices that summarise broad-scale conditions (e.g., Kumar et al., 2014; Siedlecki et al., 2023), and making modest adjustments to their plans. This is especially true for industries sensitive to mesoscale ocean variability, including oil and gas, fisheries, search and rescue, shipping, and Defence.

For example, the operator of a fishing vessel might target a specific ocean feature, such as a cyclonic eddy. Their interest is likely in the location and intensity of a specific eddy (e.g., Xing et al., 2023). A shipping company might seek deviations from a direct port-to-port route, modestly adjusting their path only if favourable conditions are likely (e.g., Dickson et al., 2019). In search and rescue, identifying the current and near-future direction of surface currents is likely most crucial (e.g., Rosebrock et al., 2015). For defence, the presence or absence of strong mesoscale features could present risks or opportunities in undersea warfare (e.g., Bub et al., 2014). Oil and gas operations may need to anticipate strong currents that could jeopardise exploratory activities or to support incidents like oil spill response (e.g., Lubchenco et al., 2012). Most of these considerations may be associated with strong mesoscale eddies in their area of operation. For all these applications, a reliable forecast of SLA, from which estimates of geostrophic currents can readily be calculated, is likely to be valuable. Moreover, for many applications, a detailed forecast of subsurface properties may be unnecessary. Analog forecasting offers an approach that allows a subset of variables to be predicted, and may be a sensible option for many applications.