Paper deep dive
Optimizing Geoengineering Interventions Using Differentiable Climate Models
Pulkit Dubey, Dorian S. Abbot, Ashesh Chattopadhyay
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 92%
Last extracted: 9/16/2026, 8:56:57 AM
Summary
This paper introduces a novel geoengineering control strategy using differentiable climate models, specifically JAX-GCM, to optimize sea-surface temperature (SST) cooling patterns. By employing gradient-based optimization with a receding-horizon approach, the method designs spatiotemporal forcing patterns (via marine cloud brightening or stratospheric aerosol injection analogs) to counteract uniform ocean warming. The learned strategy successfully removes over 92% of land warming and restores hydrological variables like precipitation and evaporation when using a spatial pattern objective, demonstrating robustness when transferred to AI emulators like LUCIE and NeuralGCM.
Entities (8)
Relation Signals (6)
JAX-GCM → uses → Differentiable Climate Models
confidence 95% · In the differentiable primitive-equation atmospheric model JAX-GCM...
Receding-Horizon Control → isusedin → JAX-GCM
confidence 92% · ...we optimize greedily over segments of 8 to 14 days, following receding-horizon control.
JAX-GCM → iscomparedwith → LUCIE
confidence 90% · The learned strategy from JAX-GCM replayed in the AI emulators LUCIE and NeuralGCM without re-optimization is successful...
Marine Cloud Brightening → isanalogousto → SST Cooling Pattern
confidence 88% · This idealized set-up represents a cooling pattern that could be delivered physically either by marine cloud brightening or stratospheric aerosol injection.
Stratospheric Aerosol Injection → isanalogousto → SST Cooling Pattern
confidence 88% · This idealized set-up represents a cooling pattern that could be delivered physically either by marine cloud brightening or stratospheric aerosol injection.
JAX-GCM → issimilarto → NeuralGCM
confidence 85% · JAX-GCM [27], a differentiable physics-based atmospheric model which couples the spectral primitive-equation core of NeuralGCM [19]...
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:The deployment of a geoengineering program to cool Earth's climate may be imminent. It is crucial that tools be developed to ensure that such a program would achieve its objectives while minimizing disruption. Here we exploit recently developed differentiable atmospheric models to demonstrate a novel geoengineering control strategy. In the differentiable primitive-equation atmospheric model JAX-GCM we impose a uniform $+4$\,K ocean warming and ask what pattern of sea-surface temperature cooling -- in five ocean-masked zonal bands of prescribed SST forcings whose amplitudes are free -- returns land near-surface air temperature closest to the model's own unwarmed climatology. This idealized set-up represents a cooling pattern that could be delivered physically either by marine cloud brightening or stratospheric aerosol injection. Gradients through chaotic dynamics decorrelate from the true sensitivity beyond the Lyapunov horizon, so we optimize greedily over segments of 8 to 14 days, following receding-horizon control. The learned strategy removes $92.3 \pm 0.4\%$ of the realized land warming across a ten-member ensemble of two-year rollouts, and a three-year run sustains it. If we use the spatial pattern of land temperature as the optimization objective, the distributions of precipitation, evaporation, and specific humidity over land are restored as well, even though they are not included in the objective function. The learned strategy from JAX-GCM replayed in the AI emulators LUCIE and NeuralGCM without re-optimization is successful, suggesting robustness. These promising results demonstrate a strategy for designing optimal climate interventions that can be applied broadly for geoengineering scenarios under consideration.
Tags
Links
- Source: https://arxiv.org/abs/2609.12528v1
- Canonical: https://arxiv.org/abs/2609.12528v1
Trouble viewing inline? Open PDF directly →
Full Text
87,269 characters extracted from source content.
Expand or collapse full text
Optimizing Geoengineering Interventions Using Differentiable Climate Models Pulkit K. Dubey Affiliation: Department of Applied Mathematics, University of California, Santa Cruz, CA 95064 Dorian S. Abbot Affiliation: Department of the Geophysical Sciences, University of Chicago, Chicago, IL 60637 Ashesh Chattopadhyay Affiliation: Department of Applied Mathematics, University of California, Santa Cruz, CA 95064 Affiliation: Corresponding author. E-mail: aschatto@ucsc.edu Abstract The deployment of a geoengineering program to cool Earth’s climate may be imminent. It is crucial that tools be developed to ensure that such a program would achieve its objectives while minimizing disruption. Here we exploit recently developed differentiable atmospheric models to demonstrate a novel geoengineering control strategy. In the differentiable primitive-equation atmospheric model JAX-GCM we impose a uniform +4+4 K ocean warming and ask what pattern of sea-surface temperature cooling – in five ocean-masked zonal bands of prescribed SST forcings whose amplitudes are free – returns land near-surface air temperature closest to the model’s own unwarmed climatology. This idealized set-up represents a cooling pattern that could be delivered physically either by marine cloud brightening or stratospheric aerosol injection. Gradients through chaotic dynamics decorrelate from the true sensitivity beyond the Lyapunov horizon, so we optimize greedily over segments of 8 to 14 days, following receding-horizon control. The learned strategy removes 92.3±0.4%92.3± 0.4\% of the realized land warming across a ten-member ensemble of two-year rollouts, and a three-year run sustains it. If we use the spatial pattern of land temperature as the optimization objective, the distributions of precipitation, evaporation, and specific humidity over land are restored as well, even though they are not included in the objective function. The learned strategy from JAX-GCM replayed in the AI emulators LUCIE and NeuralGCM without re-optimization is successful, suggesting robustness. These promising results demonstrate a strategy for designing optimal climate interventions that can be applied broadly for geoengineering scenarios under consideration. Keywords: geoengineering || differentiable climate models || gradient-based optimization || climate emulators 1 Introduction Proposals to offset greenhouse warming by engineering the climate system, e.g., marine cloud brightening [1] and stratospheric aerosol injection [2], pose an inverse design problem: which spatiotemporal pattern of forcing best restores a target climate? The inverse problem has been treated as an optimization problem since the work of Ban-Weiss and Caldeira [3]. The existing solution families superpose precomputed linear response patterns [3, 4, 5], tune feedback gains on a few zonal moments [6, 7, 8], invert Green’s-function response operators [9, 10], or search by sampling with reinforcement learning [11, 12] (Section S6). Each assumes a linear response or requires an additional simulation for every added forcing parameter. Gradient-based optimization yields sensitivities to a large number of forcing parameters in a single backpropagation pass, circumventing the simulation expense. Such methods are established in geophysical modeling. Variational data assimilation minimizes forecast error by adjoint descent over a window of hours to days [13, 14]. The ECCO state estimate extends the window to decades and admits the surface forcing fields as control variables, in a coarse, non-eddying ocean where the adjoint remains stable [15]. In chaotic regimes the adjoint grows exponentially with time, and the usable optimization window shrinks to days or weeks. Conditional nonlinear optimal perturbation operates within that horizon, optimizing an initial or boundary perturbation for maximum error growth [16, 17]. Differentiable climate models supply these gradients by construction [18, 19] and have been used to produce worst-case hurricanes [20] and heatwaves [21, 22] and global emissions trajectories [23]. Each of these optimizes toward observed history or maximum error growth, and the ECCO forcing fields are constrained to stay near the observed reanalysis forcing rather than designed (Section S6). To our knowledge, no published work uses these gradients to design a spatiotemporal forcing that steers a climate simulation to a prescribed counterfactual target and holds it there over multiple years. That is the problem addressed here. Figure 1: Problem setup, optimization scheme, and outcomes. Gradient-based optimization through a differentiable climate model learns a time-varying pattern of ocean cooling that returns a warmed climate to its unwarmed target. (A) A uniform +4+4 K SST offset warms the climate; five ocean-masked zonal bands with fixed shapes and free amplitudes a1,…,a5a_1,…,a_5 supply the corrective SST forcing; the target is the model’s own unwarmed climatology. (B) Greedy segmented optimization. Amplitudes for each 8-day segment are found by backpropagation through time, applied, and the state advances to the next segment. No gradient crosses a segment boundary, yet land-mean T2mT_2m converges to the target in ten weeks and holds it. (C) Time-mean forcing learned under the pattern and land-mean objectives. (D) Zonal-mean precipitation residual for the two objectives; the pattern objective restores precipitation, which never enters the loss. (E) The learned strategy, replayed without re-optimization in LUCIE and NeuralGCM, tracks each model’s unwarmed target. Differentiability alone, however, does not make the problem tractable. A gradient backpropagated through chaotic dynamics grows as eΛTe T, where Λ is the leading Lyapunov exponent, even as the true sensitivity of a time-averaged objective remains bounded [24, 25]. Beyond a few multiples of the Lyapunov horizon 1/Λ1/ (about one week for the mid-latitude atmosphere) a single-trajectory gradient carries no usable signal. We therefore design a greedy optimizer (Fig. 1). The forcing amplitudes are optimized over one short segment by exact backpropagation through time, the converged forcing is applied, the state advances, and the loop repeats. No gradient crosses a segment boundary or extends beyond two weeks, yet the concatenated strategy steers a perturbed climate to its target and sustains it over the following years (Sec. 2.1). This scheme is called the receding-horizon approach of adjoint-based turbulence control [26]. The gradient provides the sensitivity to every forcing parameter in one backward pass, at a cost independent of the parameter count, without a precomputed response library and without assuming that the responses to the individual bands superpose linearly. To demonstrate our optimizer’s success on designing strategies for climate intervention, we use JAX-GCM [27], a differentiable physics-based atmospheric model which couples the spectral primitive-equation core of NeuralGCM [19] to SPEEDY subgrid-scale physics [28]. We impose a uniform +4+4 K sea surface temperature (SST) warming in JAX-GCM and ask the optimizer to design the time-varying amplitudes of five forcings shaped as Gaussian zonal bands of perturbation to SST, with the objective of returning near-surface air temperature over land to the model’s own unwarmed climatology. We run the same optimization in two AI emulators, LUCIE [29, 30], built on a spherical Fourier neural operator, and NeuralGCM, a hybrid model with physics-based dynamics and learned subgrid-scale physics [19]. This setup isolates four questions that an intervention-design pipeline must answer. (a) Do backpropagated gradients reflect the true sensitivity of the objective? (b) Does stitching short-horizon solutions produce long-horizon control? (c) What does the optimizer learn beyond a uniform offset, and at what cost to the rest of the climate? (d) Does a strategy (the spatio-temporal variations of the Gaussian-shaped forcings) learned in one model apply to another? To address (a), we compare the backpropagated gradient against central finite differences across segment lengths, seasons, and operating points, and we record an ensemble noise-to-signal statistic at every segment of every run (Section S1). To address (b) we run the greedy trajectory over two and three years (Sec. 2.1). For (c) we optimize two objectives, one on the area-weighted mean of the temperature over land (land-mean objective) and one on its spatial pattern (pattern objective), and measure the response of precipitation, evaporation, and specific humidity, none of which are used in either objective (Sec. 2.2). An objective that uses a global-mean temperature target is known to distort the hydrological cycle [31, 32, 33]. We also compare the learned strategy against flat and linear-response baselines at equal forcing (Sec. 2.3). For (d) we optimize LUCIE and NeuralGCM independently and also replay the JAX-GCM strategy in each without re-optimization (Sec. 2.4, Section S2). 2 Results 2.1 Gradient-based intervention design sustains multi-year control We impose a uniform +4+4 K sea-surface temperature (SST) warming in JAX-GCM and seek to design an intervention, an additional SST perturbation added on top of the warming, that returns the near-surface air temperature over land to the model’s own unwarmed climatology (target). The intervention consists of five zonal bands of SST perturbation centered at 45∘N45 N, 20∘N20 N, 0∘0 , 20∘S20 S, and 45∘S45 S. Within each band the perturbation follows a Gaussian profile in latitude of width σ=10∘σ=10 , and is applied only over the ocean. The value of the SST perturbation at the center of each Gaussian band j=1,…,5j=1,…,5 is referred to as the amplitude aj(t)a_j(t) of the respective zonal band (Fig. 1A). Designing the intervention amounts to determining the individual band amplitudes as a function of time in order to achieve the target. However, the amplitudes are not adjusted continuously. Instead, the run is divided into consecutive 8-day blocks, which we call segments, and the amplitudes are optimized and held constant within each segment (Fig. 1B). Each aj(t)a_j(t) is therefore a piecewise-constant function of time, taking one value per segment. The optimizer selects, for each segment, the five amplitudes that minimize L=α(⟨d⟩w)2+β⟨(d−⟨d⟩w)2⟩w+μaj2¯+λ(aj−aj−)2¯,L\;=\;α\, ( d _w )^2\;+\;β\, (d- d _w )^2 _w\;+\;μ\, a_j^2\;+\;λ\, (a_j-a_j^- )^2, (1) where ⟨⋅⟩w · _w is the area-weighted mean over land and d()=Topt()−Ttgt()d(x)=T_opt(x)-T_tgt(x) is the difference between the near-surface air temperature of the optimized run and that of the unwarmed target, each time-averaged over the segment. The overbar is the mean over the five bands, aja_j are the current segment’s constant amplitude values, and aj−a_j^- are the converged amplitudes from the previous segment. The optimizer determines the amplitudes, evolves the model over a segment’s duration with the selected amplitudes held constant, and updates the amplitudes at the beginning of the following segment. The amplitudes (aj(t)a_j(t)) prescribed by the optimizer differ from the forcing applied on the climate. A band perturbs only the region of the ocean it occupies, and each latitude contributes in proportion to its surface area, which decreases toward the poles as the cosine of latitude (cosφ ). The applied forcing profile is therefore the summed prescribed profile weighted by the zonal-mean ocean fraction and by cosφ (Eq. (S3)). Band j’s contribution to the global-mean applied SST perturbation is thus given by Fj(t)=aj(t)AjF_j(t)=a_j(t)A_j, where AjA_j is the band’s effective ocean area as a fraction of global area (Eq. (S4), Section S7). The five amplitude histories together with the fixed spatial pattern of the bands yield the spatio-temporal variation of the applied forcing, which we call the strategy. The squared land-mean bias and the spatial variance of the residual, d, measure the distance of the optimized state from the target. For α=β=1α=β=1 they sum to the land-mean squared error ⟨d2⟩w d^2 _w (Section S10). The land-mean objective (β=0β=0) retains only the bias term, while the pattern objective (α=1α=1, β=0.5β=0.5) penalizes the spatial structure of the residual as well. The last two terms regularize the intervention, penalizing large amplitudes (μ) and their change between segments (λ). The segments are optimized in succession over multi-year runs (Materials and Methods). Fig. 1 summarizes the method. The following conventions are followed unless stated otherwise. The unwarmed control and the uncontrolled +4+4 K reference are five-member ensemble means. Pattern-objective results are averages over ten ensemble members that differ in their optimizer initialization, land-mean-objective over five, and shading or error bars span ±1± 1 standard deviation across members. All statistics are time means over the final 364 days, one seasonal cycle, of runs two years or longer. Figure 2 shows the temperature reduction achieved with 8-day optimization windows over a two-year run (91 segments). The optimized run converges toward the target climatology within about 10 segments, then tracks the unwarmed seasonal cycle for two years, removing over 90% of the land-mean warming realized in the uncontrolled +4+4 K run. The land-mean trajectories of the ensemble members differ only in the initial 10 weeks and converge beyond it. The ensemble members are generated by varying the initial band amplitudes. Ensembles generated by varying the initial condition instead show a statistically indistinguishable spread (Section S8). The optimized zonal-mean near-surface temperature profile follows the target at all latitudes (Fig. 2(b)). The outcome does not depend on the optimization horizon up to a limit (Section S1). A 14-day-segment ensemble over three years (Fig. S9) converges to a similar spatiotemporal variation of band amplitudes and removes a comparable fraction of the warming (Section S8). The accuracy and precision of the backpropagated gradients, which underlie the optimization, were measured using finite-difference tests and noise-to-signal ratios, respectively. Both degrade beyond a 14-day horizon and fail by about 28 to 30 days (Section S1). The time-mean strategy the optimization produces is shown in Fig. 2(d). Each filled Gaussian profile denotes a zonal band’s latitudinal extent, with its height indicating the band’s time-averaged prescribed amplitude. The thick black line is the corresponding applied forcing profile (Eq. (S3)). The prescribed amplitude and the applied forcing differ most for the 45∘N45 N band, which has the highest prescribed amplitude but the smallest applied forcing, as it acts on the smallest ocean area. Figure 2: Segment-wise control with the temperature pattern based objective via sea-surface-temperature (SST) forcing in JAX-GCM. Each segment (optimization window) is 8 days, but the optimizer reduces the mean temperature to the target over 10 weeks and sustains it over two years. All time series plots are shown as a moving average over 5 segments (40 days). (a) 2m temperature averaged over all land areas, T2mT_2m (K), vs time (months) for the warmed (+4+4 K), target, and optimized trajectories. The optimized curve shows every fifth segment with a dot. (b) Zonal-mean T2mT_2m (K) against latitude. The optimized zonal mean is close to that of the target (unwarmed) state. (c) Temporal variations in the learned strategy. Each band’s cooling amplitude is shown as the deviation from its time-averaged value in (d), drawn about the band’s own latitude; the bar at the upper right gives the scale, ±1± 1 K. (d) Time-averaged values of the applied cooling against latitude, with the average taken over the full two-year run. The filled curves denote the five Gaussian bands with the amplitudes indicated by their heights. Thick black line shows the ocean-fraction and area-weighted forcing profile of Eq. (S3). 2.2 A temperature pattern objective restores precipitation and evaporation as well Both objectives of Eq. 1 act only on near-surface temperature over land. The land-mean objective (α=1α=1, β=0β=0) penalizes only the difference between the segment-averaged land-mean temperature of the optimized run and that of the target, thus targeting only for the spatiotemporal average of near-surface temperature over land return to the target value. The pattern objective (α=1α=1, β=0.5β=0.5) is a mean-squared error modified to weigh the land-mean bias and the spatial structure of the land temperature residual separately, the bias at full weight and the spatial variance of the residual at half weight (Section S10). Results are robust to minor changes in their relative weights. While both strategies remove a comparable fraction of the realized land-mean warming (Section S8), they achieve it via different strategies. The land-mean optimizer concentrates its cooling on the equator (Fig. 3(a)). Its 0∘0 band takes the deepest time-mean amplitude and carries 34% of the total applied forcing (Fig. 5(b), Eq. (S4)). The pattern optimizer applies shallower cooling at the equator and spreads the rest across the bands (Fig. 3(a)), with an effort about 10% larger than the land-mean strategy’s (Eq. 4). The two strategies thus reach the same land-mean temperature while placing the forcing at different latitudes. The consequences of the differences in forcing structure between the two objectives extend to variables unconstrained by the objective (loss) function. Out-of-loss variables test the extent to which a strategy restores the full climate as opposed to merely transferring residuals to unconstrained fields. The in-loss variable in both objectives is only the near-surface temperature over land, averaged in different ways. Since neither objective constrains the hydrological cycle, we evaluate precipitation, evaporation, and specific humidity as out-of-loss variables (Fig. 3(b), absolute profiles in the top row, departures from the unwarmed target in the bottom row). Land-mean based control leaves its largest residuals on the equator. It replaces the warmed run’s excess tropical precipitation with an equatorial deficit of comparable magnitude and opposite sign, with similar effects on specific humidity and evaporation (Fig. 3(b), bottom row). This deficit may be attributed to the heavy equatorial forcing under the land-mean objective (Fig. 3(a)). The 0∘0 band lies under the Intertropical Convergence Zone (ITCZ), and the strong time-mean cooling applied there suppresses tropical convection and dries the ITCZ. Because the objective sees only the land-mean temperature, the intervention removes the warming while disrupting the hydrological cycle. Pattern temperature control, by contrast, drives the zonal-mean residuals of all three hydrological variables being near zero, with an equatorial precipitation residual almost an order of magnitude smaller than under the land-mean temperature objective. (a) Applied forcing (b) Hydrological response Figure 3: Applied forcing and zonal-mean hydrological response under the two land-temperature objectives. The hydrological variables are out-of-loss (not explicitly restored by the optimization scheme). The pattern objective applies less cooling at the equator and restores the hydrological variables better than the land-mean objective. (a) Ocean-fraction and area-weighted forcing profile of Eq. (S3) (K) against latitude, time-averaged over the final 364 days, for the runs optimized against the temperature pattern (Pattern, blue) and against the land-mean temperature (Land-mean, red). (b) Top row: zonal-mean profiles of precipitation (m day-1), evaporation (m day-1), and specific humidity (g kg-1) against latitude for the unwarmed target (dashed), the uncontrolled +4+4 K warmed system (orange), and the same two optimized runs. Bottom row: departures of the warmed and optimized systems from the unwarmed target, with the +4+4 K departure drawn as the shaded envelope. The land-mean objective leads to a drying of the equatorial region, whereas the pattern objective keeps residuals near zero. 2.3 Optimization performance compared with baselines We compare the performance of the learned strategies using two metrics – the restoration gain G, which is the factor by which an intervention reduces the mean-square error of the time-averaged fields relative to the uncontrolled +4+4 K run, and the effort E (in K) which is the overall forcing prescribed by the intervention, measured as the spatiotemporal average of the magnitude of the forcing summed over all bands. A good strategy should deliver high gain at low effort. Fig. 4 compares the two learned strategies against four baselines on a gain (G)(G) vs effort (E)(E) ((4) and (3), Materials and Methods) plane. The four baselines fall into two groups, ablations of the learned strategy (with pattern objective) and independent references. The first reference is the uniform baseline which applies −4-4 K on all five bands, as a naive solution to the 44 K ocean warming. The second reference is the linear baseline that solves for constant amplitudes under the linear-response assumption of the optimal-pattern methods [3, 9] and runs the full nonlinear model at those amplitudes (Section S4). The two ablations of the learned pattern strategy are created to isolate the effects of spatial and temporal variation of the applied forcing. The matched-effort baseline applies one constant amplitude to all five bands (−3.72-3.72 K), chosen such that its effort equals the learned strategy’s, E=3.44E=3.44 K. Its gain measures what the learned effort achieves without either the spatial distribution or the time variation. The matched-pattern baseline constructs the amplitudes as time averages of individual amplitudes from the learned strategy with pattern objective, so it matches the learned spatial distribution of the forcing pattern and its effort but not its temporal variation. The gain difference between the learned strategy and matched-effort is the contribution of the spatial distribution. The difference between the learned strategy and matched-pattern is the contribution of the time variation. All baselines are evaluated as five-member ensembles and hold their band amplitudes constant in time. The learned strategy with pattern objective and its constant-in-time matched-pattern baseline have the highest gain for each variable. The uniform baseline spends more effort than the pattern strategy (E=3.70E=3.70 K against 3.443.44 K) and gains less on every variable, so a larger forcing does not imply better restoration. At equal effort the pattern strategy beats the matched-effort baseline on 2 m temperature (G=27.4G=27.4 against 20.820.8), evaporation (49.249.2 against 23.623.6), specific humidity (7474 against 3434), and precipitation (9.8±1.79.8± 1.7 against 6.9±0.76.9± 0.7). The matched-pattern baseline reproduces the learned gains to within the ensemble scatter, indicating that the contribution of time variation of the amplitudes (Fig. S10) is not resolved by these time-mean metrics. The land-mean strategy spends the least effort and gains the least. The linear baseline spends less effort than others (except the land-mean strategy) and is not dominated for any variable, but gains less than the learned (pattern) strategy for all variables. Figure 4: Skill and effort of the learned strategy against baselines. The learned strategy and its constant-in-time matched-pattern baseline reach the highest gain on every variable and beat the flat forcing matched to the same prescribed effort, so the spatial structure of the forcing contributes beyond the amount of cooling. Each panel shows one variable, with the restoration gain G of (3) against the prescribed effort E of (4). G is the ratio of the area-averaged mean-square error of the +4+4 K warmed run to that of the optimized run, both measured on time-mean fields over land against the unwarmed target. G=10G=10, for instance, indicates that the error between the optimized and target runs is one tenth of that between the +4+4 K warmed and the same target runs. G=1G=1 indicates no improvement. The effort E is the global-area and time mean of the magnitude of the summed prescribed cooling profile, in K. Learned (pattern) and Learned (land-mean) are the optimized strategies under the two objectives. Matched effort applies one constant amplitude on all bands with the effort of Learned (pattern). Matched pattern applies the time-mean amplitudes of Learned (pattern). Flat −4-4 K applies −4-4 K on every band. Linear solves for constant amplitudes under a linear-response assumption and runs the full model with solved amplitudes. The dashed line joins the strategies with the largest gain for a given effort. 2.4 Cross-model agreement and transfer The strategy learned in JAX-GCM is not specific to this model. We replay the 8-day band strategy of Fig. 2 through LUCIE and NeuralGCM, each running forward only with no re-optimization, alongside the same replay in JAX-GCM and compare the gain (G)(G) to those via the uniform −4-4 K baseline (Fig. 6). Since NeuralGCM becomes unstable beyond roughly one year, we restrict the runs to one-year in every model, and treat the first 60 days as transients. Fig. 6 shows the land-mean variation from the replay in different models and gain comparison between the replay and the −4-4 K baseline. Although the flat baseline is successful in restoring the land-mean temperature, the learned strategy outperforms the flat −4-4 K baseline on both temperature and precipitation pattern restoration in every model (Fig. 6). The zonal-mean profiles behind these scores, and a zonal-profile restoration skill computed from them, are shown in Section S3 (Fig. S7). Note that near-surface temperature is a different quantity on a different grid in each model, so the bar heights of Fig. 6 are computed against the model’s own unwarmed control. The transfer is consistent with part of what LUCIE and NeuralGCM learn on their own. Figure 5 compares the time-mean forcing each model learns under the two objectives. Under the land-mean objective all three concentrate their forcing on the equatorial band (Fig. 5(b)), the placement that dries the ITCZ in JAX-GCM (Fig. 3(b)). Under the pattern objective all three spread the forcing across the bands (Fig. 5(a)), although the spatial structure of forcing in LUCIE differs from that in the other two models. NeuralGCM converges on the structure found in JAX-GCM, applying the largest time-mean cooling on the 45∘N45 N band with a secondary maximum at 45∘S45 S. LUCIE places its largest cooling on the subtropical bands and applies its smallest forcing at 45∘N45 N. (a) Pattern objective (b) Land-mean objective Figure 5: Time-mean forcing in the three models under (a) the pattern objective and (b) the land-mean objective. Each model’s subplot shows the five Gaussian bands filled at that run’s time-mean amplitudes (K), the raw summed amplitude profile (dotted), and the ocean-fraction- and area-weighted forcing profile of Eq. (S3) (thick black), against latitude on a forcing scale shared across each subfigure. Models left to right in each group: JAX-GCM, NeuralGCM, and LUCIE. The pattern-objective runs are the JAX-GCM run of Fig. 2 (two years), the run of Fig. S5 (NeuralGCM, one year), and that of Fig. S6 (LUCIE, one year); the land-mean runs are a two-year JAX-GCM run identical except for its objective, a 42-segment NeuralGCM run, and a five-year LUCIE run. In (b) all three models concentrate the forcing on the equatorial band; in (a) the forcing is spread, with deepest amplitude at 45∘N45 N in JAX-GCM and NeuralGCM. Figure 6: The JAX-GCM 8-day band strategy of Fig. 2, replayed calendar-aligned through all three models (JAX-GCM, LUCIE, NeuralGCM) running forward only, alongside a flat −4-4 K baseline. Left: land-mean near-surface temperature (K) against time (months) over the replay year, with the JAX-GCM unwarmed target (grey dashed), the +4+4 K warmed reference (orange), the flat baseline in JAX-GCM (black), and the controlled trajectory in each model. Each model’s trajectory is drawn as its own anomaly against its own unwarmed control, added back onto the reference model’s curve. All time series are shown as a moving average over 15 days. Right: restoration skill G over land, the mean-square error of the time-mean field against the model’s own unwarmed control for the uncontrolled +4+4 K run, divided by the same quantity for the run shown. Dark bars are the T2mT_2m field, light bars the precipitation field. Each model carries one pair of bars for the replayed strategy and one for its own flat −4-4 K baseline; the dotted line at G=1G=1 is no improvement over the uncorrected warming. In every model and both variables the strategy scores above the baseline. 3 Discussion Our results demonstrate differentiable climate control as an inverse-design method and emphasize the importance of the choice of objective. Interestingly, matching the spatial pattern of land temperature is enough to achieve good matching of hydrological variables, even though they do not appear in the loss function. In feedback-regulated multi-latitude stratospheric aerosol injection, temperature objectives built from the global mean and two zonal moments did not by themselves restore precipitation, which was restored only when included explicitly as a control objective [33]. The temperature target here is the full land T2mT_2m pattern instead of moments, and the forcing is applied on the surface directly rather than through a constrained aerosol loading. The improvement of precipitation under this temperature-only pattern objective suggests that resolving the temperature structure associated with hydrology may provide a parsimonious objective. A lesson from this work is that objectives for intervention design should be evaluated on out-of-loss variables as well as on the quantities they optimize. Backpropagation based gradient computation scales well to higher- dimensional problems. Optimal-pattern methods estimate the response to each basis function separately and assume that the responses combine linearly [3, 4]. Feedback approaches regulate a small number of climate moments with a linearized design model [6, 8]. Sampling-based approaches incur search costs that grow with the action space [11, 12] dimension. Reverse-mode differentiation, in contrast, returns the sensitivity to every forcing parameter in one backward pass, without a parameter-wise response ensemble or a linearized model. As an example of a higher dimensional forcing experiment, we consider a forcing space of 30 two-dimensional Gaussian SST patches (instead of 5 zonal bands) such that their superposition yields the zonal basis (Section S4). The computational cost of the 30-dimensional gradient computation remains nearly equal to that for the 5-dimensional problem. However, the patch basis shows less restoration skill than the zonal basis, potentially due to a more complex optimization landscape in the higher-dimensional forcing space. The learned strategies may depend on the model that provides the gradients. Under the pattern objective LUCIE places its deepest cooling on the subtropical bands, whereas JAX-GCM and NeuralGCM place theirs at 45∘N45 N (Fig. 5(a)). A coordinated intercomparison of AI models under the same prescribed +4+4 K SST warming finds that purely data-driven emulators depart from physics-based solvers under out-of-training-distribution forcing and attributes the departure to their treatment of land cells [34]. Our objectives act on land temperature under out-of-training-distribution forcing, so LUCIE’s departure is consistent with this failure mode. The matched-pattern strategy scores nearly the same as the learned one (Sec. 2.3), even though a clear seasonal cycle exists for some band amplitudes (Section S5). In the final year the learned amplitude vector departs from the matched constant by a seasonal cycle with a peak-to-peak swing of 2.72.7 K in the 45∘N45 N band. While the cycle is clear, the gain in skill over the matched-pattern strategy is minimal because the loss is insensitive to a displacement of this size about its minimum. Broadly, the loss, as a function of the amplitude vector, is a shallow parabola in amplitude space, and the optimizer brings the amplitudes close to its bottom. The learned (seasonal) amplitude vector differs from that of the matched-pattern by a seasonal displacement of about 11 K, roughly 10%10\% of the vector’s norm. Because the amplitudes sit near a quadratic minimum, that change in the amplitude vector leads to a much smaller change in the loss terms, about 0.0070.007 K2 per segment against a weather noise of 0.50.5 K2 per segment and 0.040.04 K2 after ten members and a year. As a result, the seasonal cycle of the optimized amplitudes provides a very small additional gain relative to applying their time-mean values. The moving-horizon procedure yields a feasible rather than a provably optimal strategy. Prescribed SST cooling omits ocean adjustment, radiative forcing, aerosol microphysics, energy-budget closure, and deployment constraints. Hence, these results represent an algorithmic demonstration in an idealized inverse-design problem rather than a deployable intervention. These limitations suggest future work. The objectives can be restricted to specific areas of the land or specific time of a calendar year. Seasonally resolved objectives would show if the temporal variation of the forcing contributes to the optimization. The greedy segment-by-segment update can be replaced by model predictive control, which optimizes several segments ahead and commits only the first. The forcing dimensionality can be increased from five fixed bands to finely resolved two-dimensional patches (Section S4) or bands with widths and centers as free parameters. For computational efficiency, an L-BFGS optimizer and early truncation of the per-segment iterations could reduce the per-segment optimization cost. These extensions would move differentiable climate control from a low-dimensional demonstration toward efficient exploration of higher-dimensional intervention design spaces for which existing parameter-wise response calculations become prohibitively expensive. 4 Materials and Methods Models JAX-GCM [27] is the spectral primitive-equation core of NeuralGCM [19] coupled to SPEEDY physics [28], run at T31 resolution (96×4896× 48) with eight vertical levels, a 30-minute step, and prescribed SST. LUCIE [29, 30] is a spherical Fourier neural operator [35] trained on ERA5 [36, 29], on a 96×4896× 48 grid with a 6-hour step. NeuralGCM couples the same class of spectral core to learned physics [19] and runs at 128×64128× 64. All three are differentiable end to end in their native frameworks, and the same optimization loop runs in each (Section S2). JAX-GCM has no 2-m diagnostics. Its T2mT_2m is a fixed blend of the skin temperature and SPEEDY’s surface air temperature t0t_0, the lowest-level air temperature extrapolated to the surface along the local lapse rate, T2m=Tskin+0.585(t0−Tskin)T_2m=T_skin+0.585\,(t_0-T_skin). LUCIE outputs 2 m temperature directly, as a variable trained on the ERA5 2 m temperature. NeuralGCM has no 2-m diagnostics, and T2mT_2m is its air temperature at the lowest pressure level, 1000 hPa. Each model is optimized and evaluated against its own unwarmed control of the same variable, so absolute temperatures are not compared across models. Warming, forcing, and target The perturbed climate is a uniform +4+4 K increase of the prescribed SST. The forcing is a set of five zonal SST bands (Sec. 2.1), Gaussian in latitude with centers 45∘N45 N, 20∘N20 N, 0∘0 , 20∘S20 S, 45∘S45 S and width σj=10∘ _j=10 , applied over ocean only, each carrying an amplitude aj(t)a_j(t) held constant within an optimization segment. The target is the model’s own unwarmed control climatology, so a perfect controller returns the warmed model to its own unforced climate. Greedy optimization Each segment minimizes the objective of (1), with the land weights w normalized to sum to one. The amplitude regularizer weight is μ=0.01μ=0.01 and the movement penalty is λ=0.1λ=0.1. Each segment is optimized with Adam at learning rate 0.10.1 for 15 iterations, with the gradient averaged over a three-member ensemble of perturbed initial states, by exact backpropagation-through-time through the model. The converged amplitudes are applied, the model integrates through the segment, and the final state initializes the next segment. No gradient crosses a segment boundary. The production runs use 91 segments of 8 days (728 days); the seasonality runs (Section S5) use 78 segments of 14 days. Skill and effort Runs are scored on four variables – near-surface temperature, precipitation, evaporation, and specific humidity, based on their time-mean fields of the final 364 days of a run. For variable v the restoration skill is Rv= 1−∥x¯opt−x¯tgt∥x¯warm−x¯tgt∥,R_v\;=\;1\;-\; x_opt- x_tgt x_warm- x_tgt , (2) where the overbar is the time mean over the final 364 days, ∥⋅∥ · is the cosine-latitude-weighted spatial RMS over land (grid cells with land fraction above one half), and the three trajectories are the optimized run, the target, and the uncontrolled +4+4 K run. Rv=1R_v=1 is exact restoration of the time-mean field and Rv=0R_v=0 is no improvement. However, restoration is not linear in R. The residual error ratio is 1−Rv1-R_v, so raising RvR_v from 0.900.90 to 0.950.95 halves the remaining error, as large a multiplicative step as raising it from 0.800.80 to 0.900.90. The primary skill measure is therefore the gain Gv=1(1−Rv)2=∥x¯warm−x¯tgt∥2∥x¯opt−x¯tgt∥2,G_v\;=\; 1(1-R_v)^2\;=\; x_warm- x_tgt ^2 x_opt- x_tgt ^2, (3) the factor by which the strategy reduces the mean-square error of the time-mean field relative to the uncontrolled warmed run. Gv=10G_v=10 means the error variance is ten times smaller than under no control, and Gv=1G_v=1 is no improvement. Equal ratios of G are equal multiplicative reductions of the residual error, which is the reward structure of the restoration problem. The effort is the cooling a strategy prescribes, E=⟨|∑jgj(φ)aj(t)|⟩A,t,E\;=\; \, |\, _jg_j( )\,a_j(t)\, | _A,\,t, (4) the area-weighted global mean and time mean, over the same final 364 days, of the magnitude of the summed prescribed profile, in K, with gj(φ)g_j( ) being the unit-amplitude Gaussian profiles of Eq. (S3) before the ocean mask is applied. E measures what the controller asks for, not what the ocean receives or what a deployment would cost. The five bands do not tile the sphere, and the area mean of their summed unit profile is 0.9250.925, so prescribing a flat −4-4 K on all five bands spends E=3.70E=3.70 K rather than 44 K; only a forcing of −4-4 K at every point of the globe would yield E=4E=4. Cross-model replay The strategy learned in JAX-GCM is replayed into each model, which runs forward without re-optimization and evaluates the gain ((3)) over the final 300 days of a one-year replay. The amplitudes drive each model’s own implementation of the same SST bands of cooling. In LUCIE the SST perturbation is smoothed before entry (Gaussian smudging, σ=1.7σ=1.7 grid cells). Note that the near-surface temperature variable differs for different models (Section S3). Data, Materials, and Software Availability Optimization code, learned forcing strategies, and analysis scripts will be deposited in a public repository upon publication. Acknowledgments PD and AC acknowledge the support from the National Science Foundation (grant no. 2425667), the Sloan Foundation, and Schmidt Sciences, LLC. DA and AC acknowledges funding from the University of Chicago Climate Systems Engineering Initiative (CSEi). Computational resources were provided by NSF ACCESS MTH240019, MTH250006 and NCAR CISL UCSC0008, and UCSC0009. Author Contributions P.K.D., D.S.A., and A.C. designed research; P.K.D. performed research and analyzed data; and P.K.D., D.S.A., and A.C. wrote the paper. Competing Interests The authors declare no competing interest. Supplementary Information S1 Gradient validation A gradient backpropagated through chaotic dynamics grows as eΛTe T with the horizon T, while the true sensitivity of a time-averaged objective stays bounded, so beyond horizons of order 1/Λ1/ the gradient direction decouples from the sensitivity [24, 25]. In this section we show that the backpropagated gradients match those from finite differences (Fig. S1) for a suitable step ϵε (Fig. S2). There is a weather-dependent uncertainty in the backpropagated gradients which makes the gradients of longer segments noisy (Fig. S3). As a result the noise-to-signal ratio becomes unfavorable for runs with long segments (Fig. S4). Accuracy against finite differences The segment loss J()J(a) is the objective L of Eq. (1) evaluated for one segment run from a fixed initial condition with the five band amplitudes =(a1,…,a5)a=(a_1,…,a_5) held constant. It is the squared area-weighted land mean of the segment-mean near-surface temperature difference from the target, plus β times the spatial variance of that difference over land, plus the amplitude penalty μaj2¯μ\, a_j^2, in K2. Throughout this section α=1α=1 and μ=0.01μ=0.01, the pattern objective has β=0.5β=0.5 and the land-mean objective β=0β=0, and the movement penalty is absent (λ=0λ=0) because each case is a single segment with no predecessor. We compared the backpropagated gradient of J with respect to the five band amplitudes against central finite differences on single deterministic trajectories, at segment lengths of 7, 14, and 28 days, from four initial conditions (January 1, April 1, July 1, and October 1), at two operating points (band amplitudes all zero and all −2-2 K), and under both objectives, for 16 cases per segment length. For each case we perturbed one band amplitude at a time by ±1± 1 K about the operating point, held the other four fixed, re-ran the segment from the same initial condition, and formed the five-component central-difference gradient (gFD)j=J(+ϵj)−J(−ϵj)2ϵ,ϵ=1K, (g_FD )_j= J(a+ _j)-J(a- _j)2ε, ε=1\ K, where je_j is the unit vector of band j. Fig. S1 plots all five components of all 16 cases against the backpropagated gradient, gBPTTg_BPTT, at each segment length. At 7 and 14 days the components lie along the 11:11 line with correlations of 0.910.91 and 0.870.87, and a least-squares fit through the origin has slopes of 0.920.92 and 0.830.83, so the backpropagated gradient points the right way with a magnitude 88 to 17%17\% smaller than finite difference. The median angle between the two gradient vectors is 5.2∘5.2 and 7.6∘7.6 under the pattern objective and 5.7∘5.7 and 9.0∘9.0 under the land mean objective. At 28 days the points spread out and correlation falls to 0.490.49. The median angle between FD and BPTT gradient vectors rises to 45∘45 under the pattern objective and 41∘41 under the land mean, the cases range from 9∘9 to 86∘86 , and sign disagrees at least sometmes in three components out of five. Figure S1: Backpropagated gradient (gBPTTg_BPTT) vs. finite-difference gradient (gFDg_FD) of the JAX-GCM segment loss. The two agree in direction and to within 17%17\% in magnitude at 7 and 14 days (correlation R=0.91R=0.91 and 0.870.87) and decorrelate at 28 days (R=0.49R=0.49). Each panel plots the components of the gradient of the loss, J, with respect to the five band amplitudes, one point per band per case and 80 per panel. The cases are four initial conditions at two operating points, band amplitudes all zero and all −2-2 K, under both objectives, with ϵ=1ε=1 K on single deterministic trajectories. Colors denote bands; circles, pattern objective; triangles, land-mean objective. Grey line, 11:11; dashed line, least-squares fit through the origin; R, correlation of the 80 points. Units, K. Choice of the finite-difference step A central difference of a smooth function converges as ϵ→0ε→ 0, but here the step is bounded from below by weather noise. Fig. S2 (left) shows the segment loss, J, along the equatorial band amplitude, a3a_3, at 0.050.05 K resolution for an 8-day segment under the land-masked pattern objective, with the other four amplitudes set to zero. The loss follows a nearly linear trend that the backpropagated tangent tracks, with a scatter of σJ=0.014 _J=0.014 K2 about it (1.2%1.2\% of the range of J) because each evaluation runs a different trajectory. A finite difference over a step Δa a inherits that scatter as 2σJ/Δa 2\, _J/ a. At Δa=0.1 a=0.1 K the finite-difference slopes scatter by ±0.18± 0.18 K about a mean of 0.280.28 K such 5 of 79 have the opposite sign (negative) (Fig. S2, middle); at Δa=1 a=1 K they scatter by ±0.02± 0.02 K and converge onto a much smoother curve 12%12\% above the backpropagated gradient (Fig. S2, right), the same shortfall as the fit slopes of Fig. S1 at 7 and 14 days. The 11 K step of Fig. S1 is set by this noise floor and not by the behavior of the gradient. Figure S2: Finite-difference slope of the segment loss against the differencing step in JAX-GCM. Left: segment loss, J, of the land-masked pattern objective as a function of the equator-band amplitude a3a_3. Here we plot 81 evaluations between 00 to −4-4 K in 0.050.05 K steps for an 8-day segment with the other four amplitudes set to zero. We also plot the backpropagated tangent for at every tenth point as a line. Middle and right: points denote finite-difference slopes gFD=[J(m+Δa/2)−J(m−Δa/2)]/Δag_FD=[J(m+ a/2)-J(m- a/2)]/ a over pairs of sampled amplitudes Δa=0.1 a=0.1 K and 11 K apart, plotted at their midpoint m; line denotes the backpropagated dJ/da3dJ/da_3. Dependence on the weather realization Fig. S1 compares backpropagated and finite-difference gradients of the same trajectory, so the weather dependence of the gradient is not relevant. To measure the effect of weather we compare the backpropagated gradient of the pattern objective, with all band amplitudes at zero, across twelve initial conditions drawn from the +4+4 K warmed band-free run (Fig. S3). The first is the spun-up January 1 state that also starts the January cases of Fig. S1. Each subsequent one is the previous state advanced 7 days with the band amplitudes held at zero, so the twelve are consecutive states of one trajectory, 7 days apart and 77 days end to end. Each faint line is one backpropagated gradient vector gBPTT=(∂J/∂a1,…,∂J/∂a5)g_BPTT=(∂ J/∂ a_1,…,∂ J/∂ a_5) from one weather realization, drawn against the band centre, and the thick black line is the mean of the twelve. At 7 and 14 days every realization returns a positive gradient on all five bands, and the standard deviation across realizations is 99 to 23%23\% of the mean component at 7 days and 1818 to 40%40\% at 14 days. As an angle, a single realization’s gradient sits 5.5∘5.5 from the twelve-member mean at 7 days and 9.0∘9.0 at 14 days, after removal of a linear drift over the 77-day span, and 41∘41 at 28 days, where 8 of the 12 realizations return a negative component on at least one band. A descent step along a unit vector at angle θ to the true gradient reduces the loss by cosθ θ of the achievable amount, so 9.0∘9.0 at 14 days costs 1.2%1.2\% of the available reduction per segment and 41∘41 at 28 days costs 24%24\%. Figure S3: Weather dependence of the backpropagated segment gradient in JAX-GCM. The gradient of a single weather realization stays within 9∘9 of the twelve-realization mean up to 14 days and departs by 41∘41 at 28 days. Backpropagated gradient of the pattern objective with respect to the five band amplitudes, at all amplitudes set to zero, from twelve weather realizations (faint lines) and their mean (black), at segment lengths of 7, 14, and 28 days from left to right on a common vertical axis. Initial conditions are sampled every 7 days over a 77-day span of the +4+4 K warmed band-free run. Units, K. Noise-to-signal ratio in the production runs For segment lengths of 7, 14, and 30 days over two years (104, 52, and 24 segments) the optimizer evaluates the gradient of each segment from three initial conditions separated by a small perturbation. Each member returns its own gradient of the segment loss, and the run logs the norm of the ensemble-mean gradient (the signal) and the norm of the component-wise standard deviation across members (the noise). Fig. S4 plots the ratio of the two averaged over the 15 iterations of each segment. At 7 and 14 days every segment stays below one (medians 0.030.03 and 0.070.07, maxima 0.190.19 and 0.700.70). At 30 days the median is 1.21.2 and noise exceeds signal in 20 of 24 segments. The crossing between 14 and 30 days matches the breakdown between 14 and 28 days in Figs. S1 and S3. Figure S4: Noise-to-signal ratio of the production gradient against segment length in JAX-GCM. Every segment stays below one at 7 and 14 days and most exceed one at 30 days. Per-segment ratio of the norm of the cross-member standard deviation of the gradient to the norm of the ensemble-mean gradient, averaged over the segment’s iterations, for three runs of the five-band zonal control problem over two years at segment lengths of 7, 14, and 30 days (104, 52, and 24 segments; land-masked pattern objective with α=1α=1 and β=0.5β=0.5, 0.050.05, and 0.10.1; amplitude regularizer weight 0.01; three-member ensembles; 15 Adam iterations per segment). Blue, below one; orange, above; dashed line, noise equal to signal. The cooling degrades later than the gradient. Over the final 364 days of each run the optimizer removes 89%89\% of the land-mean warming at 7 days, 95%95\% at 14 days, and 95%95\% at 30 days (the 7-day run has β=0.5β=0.5 against 0.050.05 and 0.10.1 for the other two, so its fraction is not directly comparable), and a 50-day run not shown in Fig. S4, with every segment above one and a median ratio of 1.81.8, removes 46%46\%. The 30-day run therefore reaches its target while the noise exceeds the signal in 20 of its 24 segments. S2 Validation in LUCIE and NeuralGCM The optimization of Sec. 2.1 runs unchanged in LUCIE and NeuralGCM. Each model optimizes its own five band amplitudes through its own dynamics, toward its own unwarmed climate as target, with the land-masked objective of JAX-GCM (the pattern loss with β=0.5β=0.5, or the land-mean loss), 7- or 8-day segments, and the same amplitude regularizer. The runs shown are single realizations with zero movement penalty. Our method controls both models. Under the pattern objective each model’s optimized run converges within a few segments and then tracks its own target seasonal cycle (Figs. S5 and S6), removing 92% (NeuralGCM) and 82% (LUCIE) of the land-mean warming realized in the corresponding +4+4 K run over one year. The land-mean runs remove 90% and 91% of the land-mean warming (not shown). Under the pattern objective NeuralGCM reproduces the JAX-GCM strategy, placing its largest time-mean cooling on the 45∘N45 N band with a secondary maximum at 45∘S45 S and weaker forcing at lower latitudes (Fig. 5(a)). In contrast, LUCIE puts its largest time-mean amplitudes on the subtropical bands, −4.4-4.4 K at 20∘N20 N and −4.2-4.2 K at 20∘S20 S, and its 45∘N45 N amplitude is the shallowest of its five, −2.7-2.7 K (Fig. 5(a)). LUCIE’s disagreement may be the result of poor behavior for out-of-distribution SST forcing [34]. Figure S5: NeuralGCM optimizing its own band amplitudes under the pattern loss over one year (52 segments of 7 days, three-member gradient ensemble, no movement penalty). Left: land-mean T2mT_2m (K) against time (months) for the warmed (+4+4 K), target, and optimized runs; dots mark every fifth segment. Middle: zonal-mean temperature at 1000 hPa (K) against latitude for the same three runs, time means over the final 364 days. Right: time-mean applied cooling (K) against latitude for the ocean-fraction- and area-weighted forcing profile of (S3) (thick black), the raw amplitude profile (thin dotted), and the five Gaussian bands (filled curves). Figure S6: LUCIE optimizing its own band amplitudes under the pattern loss over one year (45 segments of 8 days, three-member gradient ensemble, no movement penalty). Panels as in Fig. S5, with the middle panel showing zonal-mean surface temperature (K). S3 Cross-model transfer of the learned strategy Figure S7 shows the zonal-mean structure behind the bars of Fig. 6 computed from the same replay. The 8-day JAX-GCM strategy of Fig. 2 is run forward in JAX-GCM, LUCIE, and NeuralGCM with no re-optimization, and each column compares that model’s controlled run against its own unwarmed control and its own uncontrolled +4+4 K run. The profiles are time-averaged over 300 days, with the first 60 days discarded as transient. Near-surface temperature is a different quantity in each model. JAX-GCM has no 2-m diagnostic and interpolates it from the skin temperature and SPEEDY’s surface air temperature t0t_0 as T2m=Tskin+0.585(t0−Tskin)T_2m=T_skin+0.585\,(t_0-T_skin), LUCIE carries the ERA5 2-m temperature as a model variable, NeuralGCM reports its lowest-level air temperature at 1000 hPa. Absolute temperatures should therefore not be compared across columns, and when NeuralGCM optimizes for itself (Section S2) its objective is its 1000-hPa temperature rather than the JAX-GCM T2mT_2m. Precipitation and evaporation are in m day-1 in all models. The annotated R is the restoration skill of the time-mean zonal profile, R=1−RMS(controlled−control)/RMS(warmed−control)R=1-RMS(controlled-control)/RMS(warmed-control), with a cosine-latitude-weighted RMS over latitude, the zonal-profile counterpart of the field skill RvR_v of the Materials and Methods. R=1R=1 is full restoration and R=0R=0 is no restoration with respect to the warmed initial state. Only T2mT_2m over land entered the loss that produced the strategy, and only in JAX-GCM. In all three models the replayed strategy restores the zonal temperature profile as well as the the precipitation and evaporation profiles, which were not included in the loss function. Restoration across models including variables not included in the loss function indicates that the strategy corrects the forced response and not the temperature field alone. Figure S7: Zonal-mean, time-mean profiles over days 60–360 of the replay year, one column per host, for near-surface temperature (top, K), precipitation (middle, m day-1), and evaporation (bottom, m day-1). Lines represent the host model’s own unwarmed control (grey dashed), the uncontrolled +4+4 K run (orange) and the frozen JAX-GCM strategy replayed in the host (blue). The annotated R is the zonal-profile restoration skill defined in the text. S4 Optimizer variants and baseline construction This section constructs the baselines of Fig. 4 (Sec. 2.3) and compares variants of the optimizer on the same effort–gain plane (Fig. S8). The effort E is the prescribed cooling ∑jgj(φ)aj(t) _jg_j( )\,a_j(t) over the final 364 days with gjg_j being the Gaussian band profiles and aj(t)a_j(t) the amplitudes (Eq. (4)). The gain G is the measure of Fig. 4 ((3)) evaluated per variable over the final 364 days. RvR_v is the fractional reduction of the cosine-latitude-weighted spatial RMS error over land of the time-mean field relative to the uncontrolled +4+4 K run, Gv=1/(1−Rv)2G_v=1/(1-R_v)^2 is the factor by which an intervention reduces the mean-square error of the time-mean field (Eqs. (2) and (3)). The 8-Day pattern strategy (ten-member ensemble, pattern-loss configuration, Sec. 2.1) and the land-mean strategy (five-member ensemble, land-mean-objective Sec. 2.2) are learnt. The −4-4 K baseline applies −4-4 K on all five bands. The matched-effort baseline has an amplitude of −3.72-3.72 K for all bands, chosen so that its effort equals the learned ensemble’s (E=3.44E=3.44 K in both). Its amplitude is set by the learned result, so it is not an independent design; at fixed effort it isolates what the spatial structure of the forcing contributes. The matched-pattern baseline holds the learned pattern strategy’s time-mean amplitudes constant and isolates what the time variation contributes. Each flat baseline is a five-member ensemble generated by varying the spin-up length from one to five years under identical forcing. The linear baseline replaces the search over amplitudes with a single inversion of a response operator, and stands in for the optimal-pattern and response-function methods of Section S6 [3, 9]. Collect the residual field d of (1) at the n=4608n=4608 grid points in d, and the land weights, normalized to sum to one, in w. The two temperature terms of (1) are then a single quadratic form ⊤Md^\! Md, with M=βdiag()+(α−β)⊤,M\;=\;β\,diag(w)\;+\;(α-β)\,ww^\! , (S1) the rank-one term carrying the squared land-mean bias and the diagonal term the spatial variance of the residual. Modeling the residual as linear in the amplitudes about the uncontrolled +4+4 K state, ()≈0+d(a) _0+G\,a, where 0d_0 is the residual of that uncontrolled run and column j of the n×kn× k operator G is the time-mean T2mT_2m response per unit amplitude on band j, the minimizer of (1) follows in closed form, ⋆=−(⊤M+μkI)−1⊤M0,a \;=\;- (G^\! M\,G+ μk\,I )^-1G^\! M\,d_0, (S2) where G is the response operator and not the gain G of (3). The matrix inverted in (S2) is k×k× k with k=5k=5 bands. The amplitude regularizer of (1) is a mean over the five bands, so it enters (S2) as μ/kμ/k, and the movement penalty λ vanishes for amplitudes held constant. Each column of G was estimated from forward runs by central differences, j=[(+εj)−(−εj)]/2εG_j=[d(+ _j)-d(- _j)]/2 with ε=4 =4 K, following the Green’s-function construction [37]. Each entry is a difference of two 364-day time means, and each perturbation run is a 24-member ensemble over perturbed initial states in which member i of the +ε+ run shares its initial state with member i of the −ε- run. The operator took 288 forward integrations of 728 days and no gradient. The full nonlinear model was then run at the fixed ⋆a , as a five-member spin-up ensemble like the flat baselines. The response is not linear across [−ε,+ε][- ,+ ]. Warming the tropical ocean produces about twice the land-weighted mean response per K that cooling it produces, so ⋆a depends on which secant estimates G. This baseline is a linear-response reference rather than a converged optimum. A comparison of the gain obtained via the optimizer against the baselines discussed above is shown in Fig. 4. Fig. S8 varies the optimizer rather than the baseline, with the segment length at 4, 8, 14, and 28 days, the land-mean objective, a three-band actuator (centers 45∘N45 N, 0∘0 , 45∘S45 S), a seven-band actuator with centers every 20∘20 from 60∘S60 S to 60∘N60 N, and a patch actuator that opens the five bands to longitudinal variation, six patches Among the segment lengths, overall skill (the four-variable mean of RvR_v) is highest for the 8-day variant. On 2 m temperature the 8-, 14-, and 28-day segments do nearly equally well, but the longer segments fall off on the other variables and the 8-day segment leads consistently across all variables. Raising the amplitude regularizer (μ) far above its adopted weight shrinks the prescribed amplitudes and the associated restoration skill toward zero. At μ=0.5μ=0.5, fifty times the adopted weight (with 14-day segments), the prescribed effort falls to E=1.1E=1.1 K and the gain to G=2.2G=2.2 on 2 m temperature and G=2.2G=2.2 on precipitation, averaged over three ensemble members (not shown in Fig. S8). Since the gradients of the loss function with respect to all control parameters are obtainable in a single backward pass, the optimization algorithm scales well to high-dimensional control and the computational cost of the 30-patch problem remains similar to the 5-band problem. However, the added degrees of freedom does not aid the optimization as the patch basis shows less skill than the 5-band basis across all variables. This may be due to the design of the objective function or convergence toward a local solution. While the increased degrees of freedom in the patch basis do not improve the optimizer skill, increasing the number of zonal bands has an outsized effect. For the 7-band basis, although the precipitation and evaporation gains saturate, the 2 m temperature gain increases by a factor of 2.7 and the specific humidity gain by 1.3. Increasing the effort in the 5-band basis (via decreasing μ in the loss function) to match that of the 7-band solution yields only marginal improvements in the restoration skill (not shown). Thus the increased band count and greater coverage over the ocean area are the primary drivers for high gain in the 7-band basis. Conversely, replacing the five bands with three has a large adverse effect as the precipitation gain falls to 1 (no improvement) and the evaporation and humidity gains to 3.5 and 5, while still removing about three quarters of the land-mean warming. This suggests that latitudinal structure as detailed as the physical cooling mechanism (marine cloud brightening or stratospheric aerosols) allows should be exploited. Sweeps not shown in Fig. S8 leave the result unchanged. Across the optimizer (Adam against L-BFGS-B, and best-iterate against last-iterate selection), the movement penalty in weight and norm, and the amplitude regularizer within its adopted range, the five-band amplitude structure keeps a similar shape, and the ensemble-mean restoration skill is unchanged to within the member-to-member scatter. The robustness statement applies to the neighborhood of the adopted operating point, not to arbitrary parameter values. Figure S8: Restoration gain G against prescribed effort E (K) for variants of the optimizer in JAX-GCM, one panel per variable, computed over the final 364 days of each 728-day run. The variants are the segment lengths (4, 8, 14, and 28 days) under the pattern objective, the land-mean objective, three- and seven-band actuators, and an actuator of 30 two-dimensional ocean Gaussian patches, the five bands split into six longitudes each. Markers are ensemble means where more than one member exists (8- and 14-day ten members, 28-day four, Land-mean and 3 Band five, 7 Band and Patch three) with bars one standard deviation; the 4 Day variant is a single run. The dashed line joins the strategies not dominated in that panel. Figure S9: Three-year control in JAX-GCM at the 14-day optimization horizon, a ten-member ensemble with the loss and penalty of the Fig. 2 ensemble, 78 segments of 14 days. All time series plots are shown as a moving average over 5 segments (70 days). (a) Land-mean T2mT_2m (K) against time (months) for the warmed (+4+4 K), target, and optimized runs. The optimized curve shows every fifth segment with a dot. (b) Zonal-mean T2mT_2m (K) against latitude for the same three runs. (c) Temporal variations in the learned strategy. Each band’s cooling amplitude is shown as the deviation from its time-averaged value in (d), drawn about the band’s own latitude; the bar at the upper right gives the scale, ±1± 1 K. (d) Time-averaged values of the applied cooling against latitude, with the average taken over the full three-year run. The filled curves denote the five Gaussian bands with the amplitudes indicated by their heights. Thick black line shows the ocean-fraction and area-weighted forcing profile of Eq. S3. S5 Learned amplitudes and seasonality The 14-day ensemble of Fig. S9 has ten members initialized with different band amplitudes, with 78 segments of 14 days, the pattern loss, and the movement penalty λ=0.1λ=0.1. Fig. S10 shows the learned band amplitudes and forcings of this ensemble. The optimizer does not prescribe equal amplitudes across bands (Fig. 10(a)) as the 45∘N45 N band maintains an amplitude about 1-2 K larger than every other band. The realized applied forcing Fj(t)=aj(t)AjF_j(t)=a_j(t)A_j of (S4) (Fig. 10(b)) has a lower spread between time-averaged amplitude means. The 45∘N45 N band with the largest amplitude has the lowest applied forcing and the 45∘S45 S band with greatest ocean fraction has the highest applied forcing. The 45∘N45 N and the 45∘S45 S bands carry a clear seasonal cycle and are anti-phase with each other, but for the current objective its effect is not resolvable beyond the noise. The learned strategy’s gains over uniform cooling therefore come from its spatial structure, not from its time variation. A constant-in-time strategy built from the time-mean amplitudes restores temperature and the auxiliary variables as well as the full time-varying strategy does. The absence of a time-variation gain does not show that the seasonality is dynamically irrelevant. The cycle is reproducible across members, but the loss is a shallow quadratic in the amplitudes near its minimum, so the seasonal displacement of about 11 K changes the segment loss by about 0.0070.007 K2, below the weather noise of 0.50.5 K2 per segment and 0.040.04 K2 after ten members and a year (Sec. 3). The tie therefore measures the flatness of the loss, not the effect of the seasonal forcing on the climate. (a) (b) Figure S10: The learned strategy over three years in JAX-GCM for a ten-member ensemble run with identical settings. The ten thin lines are the individual members, the thick line is the ensemble mean, and both are drawn under a 90-day centered moving average. (a) Prescribed amplitude (K) of each of the five Gaussian cooling bands against time (months). (b) The ocean-area-weighted forcing Fj(t)=aj(t)AjF_j(t)=a_j(t)\,A_j each band delivers to the global-area mean (K, (S4)), where AjA_j is band j’s effective ocean area as a fraction of global area. Note the change of vertical scale between the panels. S6 Prior approaches to the intervention inverse problem The forward problem, the climate response to a prescribed intervention, is well-studied through coordinated model experiments [38]. The solution families for the inverse problem make the search tractable in different ways. Optimal-pattern studies superpose precomputed single-pattern response runs and assume the responses add linearly [3, 4, 5]. Feedback controllers adjust a few zonal moments of the forcing with gains tuned on a linearized design model [6, 7, 8]. Inverse methods estimate a linear response operator from Green’s-function ensembles and invert it [9, 10]. Reinforcement learning treats the climate model as a black box and searches by sampling [11, 12]. In feedback-controlled multi-latitude aerosol injection, temperature objectives built from a few zonal moments do not restore precipitation, and precipitation must be targeted explicitly [33]. Re-optimizing the same aerosol-injection objectives in two versions of one climate model gives substantially different strategies [39], so a strategy learned in one model cannot be assumed to transfer to another. Among gradient methods, variational data assimilation minimizes forecast error by adjoint descent over windows of hours to days. The ECCO ocean state estimate extends the window to decades and includes the surface forcing fields among its adjoint control variables [15], but treats the forcing as a nuisance parameter regularized toward a prior. Conditional nonlinear optimal perturbation finds the initial or boundary perturbation, sea-surface temperature (SST) patterns included, that maximizes short-range error growth [16, 17]. Model gradients have recently been used to optimize initial conditions for worst-case heatwave storylines [21, 22], and those optimizations end within the deterministic forecast range. In each case the optimized quantity is a fit to observed history, a short-horizon amplification, or an initial state, not a multi-year intervention forcing. S7 Area weighting of the band forcing Figs. 2, 5 and S9 show the applied cooling weighted by the ocean area for each band with a black curve. To compute this net applied forcing, let aja_j denote the time-mean amplitude of band j in K, gj(φ)=exp[−(φ−φj)2/2σj2]g_j( )= \! [-( - _j)^2/2 _j^2 ] its Gaussian profile about the center latitude φj _j (σj=10∘ _j=10 ), and o(φ)o( ) the zonal-mean ocean fraction. With cosφ weighting each latitude by its area, the thick black line is F(φ)=∑jajgj(φ)o(φ)cosφ,F( )= _ja_j\,g_j( )\,o( ) , (S3) and the thin dotted line in the same panel is the unweighted profile ∑jajgj(φ) _ja_jg_j( ). Band j’s contribution to the global-area-mean applied SST perturbation is Fj=ajAj,Aj=∑φgj(φ)o(φ)cosφ∑φcosφ,F_j=a_jA_j, A_j= _ g_j( )\,o( ) _ , (S4) where the sum is over the grid latitudes and AjA_j is band j’s effective ocean area as a fraction of the global area, so that ∑jFj _jF_j is the global-area-mean applied SST perturbation in K. Since AjA_j is fixed in time, the relation also holds at each instant, Fj(t)=aj(t)AjF_j(t)=a_j(t)\,A_j. Table S1 lists the time-mean amplitudes, band area fractions, and per-band forcings for the pattern-loss run of Fig. 2 and for the land-mean-loss run (the JAX-GCM panel of Fig. 5(b)). The pattern strategy is close to uniform in delivered forcing; its 45∘N45 N band pairs the deepest prescribed amplitude with the smallest effective ocean area. The land-mean strategy varies meridionally, with its equatorial band carrying 34% of the total applied forcing. The 14-day-horizon pattern run of Fig. S9 has a similar spread as the 8-day run. Thus strategy is thus robust to variation in horizon length. Table S1: Time-mean band amplitudes and delivered forcings for the two loss designs. aja_j is the time-mean prescribed amplitude of band j, AjA_j the band’s effective ocean area as a fraction of global area, and Fj=ajAjF_j=a_jA_j the band’s contribution to the global-area-mean applied SST perturbation ((S4)). Pattern loss, the 8-day run of Fig. 2; land-mean loss, the run of the JAX-GCM panel of Fig. 5(b). Amplitudes are time means over the final 364 days of the five-segment moving-average series. Pattern loss Land-mean loss Band AjA_j (%) aja_j (K) FjF_j (K) aja_j (K) FjF_j (K) 45∘N45 N 7.7 −5.42-5.42 −0.42-0.42 −3.69-3.69 −0.29-0.29 20∘N20 N 13.9 −3.55-3.55 −0.49-0.49 −3.57-3.57 −0.50-0.50 0∘0 16.8 −3.19-3.19 −0.54-0.54 −4.38-4.38 −0.74-0.74 20∘S20 S 16.1 −3.30-3.30 −0.53-0.53 −2.60-2.60 −0.42-0.42 45∘S45 S 14.3 −3.75-3.75 −0.54-0.54 −1.63-1.63 −0.23-0.23 Sum −2.52-2.52 −2.17-2.17 S8 Control-run statistics The ten-member 8-day ensemble of Fig. 2 (per member, 91 segments of 8 days, three-member gradient ensemble, movement penalty λ=0.1λ=0.1) has removal fractions of 92.3±0.4%92.3± 0.4\% and a 2 m temperature gain ((3)) of G=27.4±3.8G=27.4± 3.8. The members of this ensemble differ in the optimizer’s initial amplitudes and share one base weather realization; a six-run test that instead varied the initial condition, with the optimizer seed fixed, gave a statistically indistinguishable spread and ensemble-mean restoration skills agreeing to within 2% for the four scored variables. The 14-day-segment ensemble of Fig. S9, with the same loss and penalty, removes 94.6±1.1%94.6± 1.1\% over three years, with G=22.2±2.6G=22.2± 2.6. The median in-run noise-to-signal ratio over the last ten iterations of each segment, for a single member, is 0.075 in the 8-day run and 0.24 in the 14-day run. All 91 segments of the 8-day run and 77 of 78 segments of the 14-day run have noise-to-signal ratio less than one. Between the 8- and 14-day runs, four of the five time-mean band amplitudes agree to within about 0.2 K, the 45∘S45 S amplitude differs by 0.5 K, and both runs place their deepest cooling on the 45∘N45 N band, where their amplitudes agree to within 0.050.05 K. The five-member land-mean-objective ensemble of Sec. 2.2, otherwise identical, removes 93.2±0.3%93.2± 0.3\% of the land-mean warming against the pattern ensemble’s 92.3±0.4%92.3± 0.4\%, with G=11.9±0.5G=11.9± 0.5 on 2 m temperature against 27.4±3.827.4± 3.8. S9 Replay scores and out-of-distribution warming On the T2mT_2m field the replayed strategy of Fig. 6 scores G=25.3G=25.3, 2.32.3, and 21.821.8 in JAX-GCM, LUCIE, and NeuralGCM, against 16.916.9, 0.80.8, and 10.710.7 for the models’ own uniform −4-4 K baselines, advantages of 1.5×1.5× to 2.7×2.7×. On the precipitation field the strategy scores 7.27.2, 1.81.8, and 2.02.0 against 5.25.2, 1.01.0, and 1.11.1, advantages of 1.4×1.4× to 1.8×1.8×. A 1.5-year repetition, whose 364-day averaging window is a full seasonal cycle containing no ramp-up, reproduces the ordering in the two hosts that completed it, with strategy-over-baseline advantages of 1.29×1.29× on temperature and 1.61×1.61× on precipitation in JAX-GCM and 2.00×2.00× and 1.43×1.43× in LUCIE; the NeuralGCM run of the repetition becomes unstable after about 14 months. Measured against each host’s own unwarmed control and averaged after the 60-day spin-up of the replay year, the flat −4-4 K run in JAX-GCM sits within 0.050.05 K of the target land mean, while the learned strategy retains a warm bias of +0.25+0.25 K in JAX-GCM, +0.04+0.04 K in LUCIE, and +0.23+0.23 K in NeuralGCM. Thus the uniform baseline is more effective in restoring the land-mean bias in JAX-GCM and NeuralGCM, but is much less effective for pattern and out-of-loss variables. Out-of-distribution behavior differs by model class. Under the uniform +4+4 K ocean warming, the land response is +3.8+3.8 K in JAX-GCM and +3.9+3.9 K in NeuralGCM. LUCIE damps the land response to +1.8+1.8 K, half the warming of the other two. This is consistent with a coordinated intercomparison of AI models under the same prescribed +4+4 K SST forcing, in which only a hybrid physics-AI model reorganized regional climates in agreement with physics-based models, while data-driven models are shown to struggle with land response to ocean warming [34]. S10 Decomposition of the temperature loss The two temperature terms of Eq. (1) are the bias and variance parts of the land-weighted mean squared error of the segment-mean temperature. Let Topt()T_opt(x) and Ttgt()T_tgt(x) be the near-surface air temperature of the controlled run and of the target, each averaged in time over one segment, and let d=Topt−Ttgtd=T_opt-T_tgt. The spatial mean ⟨f⟩w=∑w()f() f _w= _xw(x)\,f(x) uses land area weights normalized to ∑w=1 _xw=1, so ⟨⋅⟩w · _w is linear and maps constants to themselves. Splitting d into its land mean and the variance about it, d=⟨d⟩w+(d−⟨d⟩w)d= d _w+(d- d _w), squaring, and averaging gives ⟨d2⟩w=(⟨d⟩w)2+⟨(d−⟨d⟩w)2⟩w, d^2 _w\;=\; ( d _w )^2\;+\; (d- d _w )^2 _w, (S5) since the cross term vanishes, 2⟨d⟩w⟨d−⟨d⟩w⟩w=2⟨d⟩w(⟨d⟩w−⟨d⟩w)=02\, d _w d- d _w _w=2\, d _w ( d _w- d _w )=0. The first term is the squared land-mean bias and, by linearity, ⟨d⟩w=⟨Topt⟩w−⟨Ttgt⟩w d _w= T_opt _w- T_tgt _w, the squared difference of the two land means. The second term is the spatial variance of the residual. The land-mean objective (β=0β=0) therefore compares one number per segment, the pattern objective (α=1α=1, β=0.5β=0.5) keeps the variance part at half weight, and α=β=1α=β=1 recovers the mean squared error exactly. References [1] J. Latham et al. (2012) Marine cloud brightening. Philosophical Transactions of the Royal Society A 370, p. 4217–4262. Cited by: §1. [2] P. J. Crutzen (2006) Albedo enhancement by stratospheric sulfur injections: a contribution to resolve a policy dilemma?. Climatic Change 77, p. 211–220. Cited by: §1. [3] G. A. Ban-Weiss and K. Caldeira (2010) Geoengineering as an optimization problem. Environmental Research Letters 5 (3), p. 034009. Cited by: §1, §2.3, §3, §S4, §S6. [4] D. G. MacMartin, D. W. Keith, B. Kravitz, and K. Caldeira (2013) Management of trade-offs in geoengineering through optimal choice of non-uniform radiative forcing. Nature Climate Change 3, p. 365–368. Cited by: §1, §3, §S6. [5] E. Brody, Y. Zhang, D. G. MacMartin, D. Visioni, B. Kravitz, and E. M. Bednarz (2025) Using optimization tools to explore stratospheric aerosol injection strategies. Earth System Dynamics 16, p. 1325–1341. Cited by: §1, §S6. [6] B. Kravitz, D. G. MacMartin, D. T. Leedal, P. J. Rasch, and A. J. Jarvis (2014) Explicit feedback and the management of uncertainty in meeting climate objectives with solar geoengineering. Environmental Research Letters 9 (4), p. 044006. Cited by: §1, §3, §S6. [7] D. G. MacMartin, B. Kravitz, D. W. Keith, and A. Jarvis (2014) Dynamics of the coupled human–climate system resulting from closed-loop control of solar geoengineering. Climate Dynamics 43, p. 243–258. Cited by: §1, §S6. [8] B. Kravitz, D. G. MacMartin, M. J. Mills, J. H. Richter, S. Tilmes, J. Lamarque, J. J. Tribbia, and F. Vitt (2017) First simulations of designing stratospheric sulfate aerosol geoengineering to meet multiple simultaneous climate objectives. Journal of Geophysical Research: Atmospheres 122, p. 12616–12634. Cited by: §1, §3, §S6. [9] J. Lu, F. Liu, L. R. Leung, et al. (2020) Neutral modes of surface temperature and the optimal ocean thermal forcing for global cooling. npj Climate and Atmospheric Science 3, p. 9. Cited by: §1, §2.3, §S4, §S6. [10] H. Ren, J. Lu, Z. J. Hou, T.-C. Chen, L. R. Leung, and F. Liu (2024) Neural networks to find the optimal forcing for offsetting the anthropogenic climate change effects. Artificial Intelligence for the Earth Systems 3 (3). Cited by: §1, §S6. [11] C. Schroeder de Witt and T. Hornigold (2019) Stratospheric aerosol injection as a deep reinforcement learning problem. Note: arXiv:1905.07366 (ICML 2019 Climate Change AI workshop) Cited by: §1, §3, §S6. [12] H. Quan, D. D. B. Koll, N. Lutsko, and J. Yuval (2025) Solar geoengineering strategies based on reinforcement learning. Journal of Geophysical Research: Atmospheres 130, p. e2025JD044319. Cited by: §1, §3, §S6. [13] O. Talagrand and P. Courtier (1987) Variational assimilation of meteorological observations with the adjoint vorticity equation. I: theory. Quarterly Journal of the Royal Meteorological Society 113, p. 1311–1328. Cited by: §1. [14] F. Rabier, H. Järvinen, E. Klinker, J. Mahfouf, and A. Simmons (2000) The ECMWF operational implementation of four-dimensional variational assimilation. I: experimental results with simplified physics. Quarterly Journal of the Royal Meteorological Society 126, p. 1143–1170. Cited by: §1. [15] G. Forget, J. Campin, P. Heimbach, C. N. Hill, R. M. Ponte, and C. Wunsch (2015) ECCO version 4: an integrated framework for non-linear inverse modeling and global ocean state estimation. Geoscientific Model Development 8, p. 3071–3104. Cited by: §1, §S6. [16] M. Mu, W. Duan, and B. Wang (2003) Conditional nonlinear optimal perturbation and its applications. Nonlinear Processes in Geophysics 10, p. 493–501. Cited by: §1, §S6. [17] X. Song, M. Mu, R. Zhan, and Y. Gao (2026) SST sensitivity of rapid intensification in Typhoon Nanmadol (2022) revealed by CNOP. Journal of Geophysical Research: Atmospheres. Cited by: §1, §S6. [18] M. Gelbrecht, A. White, S. Bathiany, and N. Boers (2023) Differentiable programming for Earth system modeling. Geoscientific Model Development 16 (11), p. 3123–3135. Cited by: §1. [19] D. Kochkov, J. Yuval, I. Langmore, P. Norgaard, J. Smith, G. Mooers, et al. (2024) Neural general circulation models for weather and climate. Nature 632, p. 1060–1066. Cited by: §1, §1, §4. [20] D. A. Plotkin, R. J. Webber, M. E. O’Neill, J. Weare, and D. S. Abbot (2019) Maximizing simulated tropical cyclone intensity with action minimization. Journal of Advances in Modeling Earth Systems 11 (4), p. 863–891. Cited by: §1. [21] P. T. Vonich and G. J. Hakim (2024) Predictability limit of the 2021 Pacific Northwest heatwave from deep-learning sensitivity analysis. Geophysical Research Letters 51, p. e2024GL110651. Cited by: §1, §S6. [22] T. Whittaker and A. Di Luca (2026) Constructing extreme heatwave storylines with differentiable climate models. Weather and Climate Dynamics 7, p. 393–410. Cited by: §1, §S6. [23] C. B. Womack, S. Bouabid, A. Sokolov, P. Salunke, G. Flierl, S. D. Eastham, and N. E. Selin (2026) Optimal scenario design for climate emulation. Note: arXiv:2606.19302 Cited by: §1. [24] D. J. Lea, M. R. Allen, and T. W. N. Haine (2000) Sensitivity analysis of the climate of a chaotic system. Tellus A 52 (5), p. 523–532. Cited by: §1, §S1. [25] L. Metz, C. D. Freeman, S. S. Schoenholz, and T. Kachman (2021) Gradients are not all you need. Note: arXiv:2111.05803 Cited by: §1, §S1. [26] T. R. Bewley, P. Moin, and R. Temam (2001) DNS-based predictive control of turbulence: an optimal benchmark for feedback algorithms. Journal of Fluid Mechanics 447, p. 179–225. Cited by: §1. [27] E. H. Davenport J. V. Madan et al. (2026) JCM v1.1: a differentiable, intermediate-complexity atmospheric model. Geoscientific Model Development 19, p. 6451–6466. Cited by: §1, §4. [28] F. Molteni (2003) Atmospheric simulations using a GCM with simplified physical parametrizations. I: model climatology and variability in multi-decadal experiments. Climate Dynamics 20, p. 175–191. Cited by: §1, §4. [29] H. Guan et al. (2025) LUCIE: a lightweight uncoupled climate emulator with long-term stability and physical consistency. Journal of Advances in Modeling Earth Systems 17, p. e2024MS005152. Cited by: §1, §4. [30] H. Guan, T. Arcomano, A. Chattopadhyay, and R. Maulik (2025) LUCIE-3d: a three-dimensional climate emulator for forced responses. arXiv preprint arXiv:2509.02061. Cited by: §1, §4. [31] G. Bala, P. B. Duffy, and K. E. Taylor (2008) Impact of geoengineering schemes on the global hydrological cycle. Proceedings of the National Academy of Sciences 105, p. 7664–7669. Cited by: §1. [32] K. Ricke, J. S. Wan, M. Saenger, and N. J. Lutsko (2023) Hydrological consequences of solar geoengineering. Annual Review of Earth and Planetary Sciences 51, p. 447–470. Cited by: §1. [33] W. R. Lee, D. G. MacMartin, D. Visioni, and B. Kravitz (2020) Expanding the design space of stratospheric aerosol geoengineering to include precipitation-based objectives and explore trade-offs. Earth System Dynamics 11, p. 1051–1072. Cited by: §1, §3, §S6. [34] C. C. Merchant, M. Klöwer, B. Stanley-Clamp, M. Höver, S. L. L. Michel, E. Groot, and H. M. Christensen (2026) How do AI climate models respond to warming across climate zones?. Note: arXiv:2608.17986 Cited by: §S2, §3, §S9. [35] B. Bonev, T. Kurth, C. Hundt, J. Pathak, M. Baust, K. Kashinath, and A. Anandkumar (2023) Spherical Fourier neural operators: learning stable dynamics on the sphere. In Proceedings of the 40th International Conference on Machine Learning, Cited by: §4. [36] H. Hersbach et al. (2020) The ERA5 global reanalysis. Quarterly Journal of the Royal Meteorological Society 146 (730), p. 1999–2049. Cited by: §4. [37] P. Hassanzadeh and Z. Kuang (2016) The linear response function of an idealized atmosphere. part I: construction using Green’s functions and applications. Journal of the Atmospheric Sciences 73 (9), p. 3423–3439. External Links: Document Cited by: §S4. [38] B. Kravitz, A. Robock, O. Boucher, H. Schmidt, K. E. Taylor, G. Stenchikov, and M. Schulz (2011) The Geoengineering Model Intercomparison Project (GeoMIP). Atmospheric Science Letters 12 (2), p. 162–167. Cited by: §S6. [39] J. T. Fasullo and J. H. Richter (2023) Dependence of strategic solar climate intervention on background scenario and model physics. Atmospheric Chemistry and Physics 23, p. 163–182. Cited by: §S6.