This work belongs to a line of research on quantum algorithms and numerical methods for complex systems, where computational cost is a pressing issue. Neighboring works deal with El Niño prediction, accelerating power systems, and quantum methods, but this paper stands out by transferring the sampling method from theoretical models to a realistic predictive model, the Zebiak–Cane model, with dimension on the order of $10^4$–$10^5$. This is the first test of the approach in climate modeling.
In the Zebiak–Cane (ZC) model, sea surface temperature (SST) and thermocline depth anomalies are considered. The optimal precursor is sought as the solution to the conditional nonlinear optimal perturbation (CNOP) problem: maximize the SST norm after 9 months subject to a constraint on the initial perturbation. Instead of the exact gradient, which requires an adjoint model, the sampling method uses Stokes' formula: the gradient of the objective function is expressed as the expectation of function values over random directions. In practice, the gradient is approximated by the average over n random vectors uniformly distributed on a sphere. Since samples are independent, computations parallelize easily: on a supercomputer with Intel Xeon Gold 6132 processors, parallel model runs were launched.
The spatial patterns of optimal El Niño precursors obtained by the sampling method with 1000 probes practically coincide with those from the adjoint method: a warm patch of about 0.2°C in the eastern tropical Pacific and a cold one in the central Pacific, as well as a thermocline deepening of about 100 m. With 200 probes, the main features are preserved, though a small amount of noise appears. Quantitatively, the objective function (SST norm after 9 months) is 16.8441 for the adjoint method, 16.6307 for sampling with 1000 probes (98.73%), and 15.4193 for 200 probes (91.54%). The dynamics of the Niño 3.4 index over the year are also close: the relative error is on the order of 0.01°C for 1000 probes and 0.1°C for 200 probes, which is small compared to the typical amplitude of the index (about 1°C). Computation time drops from 15 seconds (50 iterations of the adjoint method) to 3 seconds with parallel sampling, that is, a factor of 5.
The sampling method eliminates the need for an adjoint model, which is difficult to develop for coupled ocean–atmosphere models and models with parameterization. This opens the possibility of computing CNOPs for realistic climate models such as WRF or FGOALS-s2, where the adjoint approach is practically infeasible. The work also shows that accuracy is maintained even when the number of probes is significantly reduced, which is important for practical applications.
Further research could apply the sampling method to more complex forecast models, including regional atmospheric models and global climate models. Also promising is the use of the method in four-dimensional variational data assimilation (4D-Var) for coupled systems where an adjoint model is absent. It may be possible to improve the method's stability with a small number of probes.
The method will affect numerical climate and weather forecasting, especially in centers that use models without adjoint versions. It may also be useful in data assimilation and predictability assessment.
Test the method on the WRF model for mesoscale forecasting and on the FGOALS-s2 model for decadal climate prediction. Investigate stability with fewer than 200 probes and optimize the parallel implementation.
In numerical climate modeling, the problem of computing optimal perturbations for high-dimensional models without an adjoint code remains unsolved. This work directly addresses the open question of whether an adjoint-free algorithm exists for models of dimension 10^4–10^5, posed earlier.
The sampling method is based on Stokes' theorem, which relates the integral over a volume to the integral over its boundary. Here, it allows replacing the gradient with an average of function values on a sphere: an elegant mathematical trick.
In Isaac Asimov's story 'The Last Question', a supercomputer searches for a way to reverse entropy by trying possibilities. The sampling method also tries random directions, but instead of an exhaustive search, it uses statistics to quickly find the optimal perturbation.