The rapid decline of Arctic sea ice, observed by space-borne microwave satellites beginning in 1979, is one of the most robust indicators of anthropogenic climate change (Cavalieri & Parkinson, 2012;Serreze & Stroeve, 2015;Stroeve & Notz, 2018). As a result, the radiative effect of Arctic sea ice has been estimated to weaken by 0.04-0.05 W⋅m -2 ⋅ decade -1 (Duspayev et al., 2024), affecting the global energy balance. The Special Report on Ocean and Cryosphere in a Changing Climate (Pörtner et al., 2019) and the Sixth Assessment Report of the Intergovernmental Panel on Climate Change (IPCC, 2021) both highlighted that such dramatic sea-ice changes not only have caused detrimental influences on weather, ecosystems, and socioeconomics within the Arctic, but could also have exerted far-reaching impacts on mid-latitude atmospheric circulation and extreme events. Furthermore, in all future warming scenarios, state-of-the-art global climate models project the Arctic to be likely ice-free in September by 2050 (Jahn et al., 2024), and the first ice-free day could even occur before 2030 (Heuzé & Jahn, 2024).
The loss of Arctic sea ice, via a series of feedback processes, amplifies warming of Arctic near-surface air temperature by a factor of 3-4 times the global or tropical average (Chylek et al., 2022;Rantanen et al., 2022). This phenomenon is called Arctic amplification (AA). In response to sea-ice loss and the accompanied AA, which intrinsically manifests itself as the reduced meridional temperature gradient in the lower troposphere, a weakening of mid-latitude jet stream due to the thermal wind adjustment is expected. Some studies further suggested that the strength and propagation of tropospheric planetary waves can be altered, and subsequently affect the occurrence of extreme events in mid-latitudes (Cohen et al., 2014(Cohen et al., , 2018;;J. A. Francis & Vavrus, 2012, 2015;J. Francis et al., 2018). Another branch of studies highlighted the role of the stratosphere (Cohen et al., 2014;Jaiser et al., 2016;Kim et al., 2014;Liang et al., 2024;Nakamura et al., 2016;Sun et al., 2015;Xu et al., 2021Xu et al., , 2023)): the sea-ice loss and the associated diabatic heating excite vertically-propagating of planetary waves that disrupt the stratospheric polar vortex via eddy-mean flow interactions, and then the stratospheric anomalies gradually move downward to the troposphere, projecting onto a negative Northern Annular Mode (NAM, Thompson & Wallace, 2000) or a negative North Atlantic Oscillation (NAO, Barnston & Livezey, 1987). These proposed dynamical pathways link Arctic sea-ice loss to large-scale atmospheric circulation changes in midlatitudes. Climate models exhibit inconsistent atmospheric responses to Arctic sea-ice loss, hindering the interpretation of the underlying mechanisms (Cohen et al., 2020;Screen et al., 2018). For example, a group of studies using atmospheric general circulation models (AGCMs) forced by sea-ice loss conditions showed that the largescale atmospheric circulation response in the Northern Hemisphere resembles a negative NAM-like or NAO-like pattern (Liang et al., 2021;Peings & Magnusdottir, 2014;Seierstad & Bader, 2009), whereas another group presented an opposite or muted response (Cassano et al., 2014;Screen et al., 2014;Singarayer et al., 2006;Strey et al., 2010). This inconsistency not only reveals that large internal variability in the atmosphere can overwhelm the sea ice-forced circulation responses, but also reflects varied temporal and spatial distributions of sea-ice forcings prescribed to drive models, different experimental configurations, or model biases.
To address these issues, several coordinated multi-model projects have been performed (e.g., Liang et al., 2020;Ogawa et al., 2018;D. M. Smith et al., 2019). Among them, the Polar Amplification Model Intercomparison Project (PAMIP, Danabasoglu, 2019;D. M. Smith et al., 2019) of the Coupled Model Intercomparison Project phase 6 (CMIP6, Eyring et al., 2016) provides simulated results using multiple climate models forced by the same future and pre-industrial sea-ice concentration (SIC) conditions. Following the protocol of PAMIP, modeling centers simulated at least 100 ensemble members for each experiment to address the effects of internal variability. The products of PAMIP serve as a perfect testbed for studying Arctic-mid-latitude linkages and can shed insights into the origin of the inconsistencies. In this study, we use the PAMIP simulation sets from the Whole Atmosphere Community Climate Model version 6 (WACCM6, Gettelman et al., 2019), developed by the National Center for Atmospheric Research, as a reference simulation set of fully physics-based modeling. WACCM6 has a high-top configuration that includes more vertical layers in the stratosphere, presumably better capturing the stratospheric pathway linking Arctic sea-ice loss to mid-latitude circulation changes.
More recently, machine learning (ML) techniques have been applied to forecasting weather and projecting future climates (Barnes et al., 2020;Beucler et al., 2023;Eyring et al., 2024;Lai et al., 2024;Li et al., 2023;Mamalakis et al., 2020;Reichstein et al., 2019). Data-driven approaches have achieved remarkable success in forecasting weather on time-scales of a few days to a week (Bi et al., 2023;Keisler, 2022;Lam et al., 2023;Rasp et al., 2024;Schreck et al., 2024;Watt-Meyer et al., 2024). However, simulations of these ML models on longer timescale beyond the weather prediction remain challenging and have not demonstrated better performance than existing dynamical climate models. This hinders their applicability to study subseasonal-to-seasonal forecasts, interannual and decadal variability, and climate change (Brenowitz & Bretherton, 2019;Watt-Meyer et al., 2023).
Instead of using purely data-driven approaches, recent studies have begun to adopt a "hybrid" modeling approach to construct climate models that can be integrated over an extended period of time (Kochkov et al., 2024;Rasp et al., 2018;Yuval & O'Gorman, 2020). The "hybrid" configuration typically combines two components: (a) a physics-based dynamical core that solves the primitive equations for atmospheric circulation and thermodynamics, and (b) a blending of trained parametrization schemes using neural networks that emulate radiation, convective processes, subgrid-scale dynamics, etc. The coupling of dynamical and ML components with additional training efforts can alleviate the instabilities and model drift that previous data-driven models encountered (Brenowitz & Bretherton, 2019;Kochkov et al., 2024), thus allowing for longer, more stable integrations. Furthermore, hybrid models are more interpretable than pure data-driven models, as the large-scale dynamics and thermodynamics are physically resolved in the dynamical core component.
NeuralGCM, a recently developed hybrid climate model by Google Research (Kochkov et al., 2024), represents a significant advancement in atmospheric modeling. In addition to its enhanced skill in ensemble weather forecasting, NeuralGCM exhibits strong capability in simulating atmospheric variability across a broad range of time-scales-from subseasonal to decadal. Owing to its flexible architecture, which permits modification of lower boundary conditions such as SIC, NeuralGCM could be an atmospheric model used to explore atmospheric responses to Arctic sea-ice loss.
In this study, we use NeuralGCM to perform, for the first time, a set of pre-industrial and future sea-ice experiments and compare the results with simulations conducted using WACCM6. Section 2 provides an overview of the model structure and experimental design. Section 3 presents the large-scale atmospheric responses to Arctic sea-ice loss in both NeuralGCM and WACCM6, including analyses of atmospheric blocking and storm tracks over the North Atlantic and its underlying mechanisms. To test the physical-based climate model dependence, simulations from another GCM-Taiwan Earth System Model Version1 (TaiESM1)-are also analyzed. A broader discussion of the findings, future outlook, and concluding remarks are provided in Section 4. We anticipate that our analysis framework could pave a new avenue for future climate modeling studies using both new generation hybrid or ML models and physics-based models.
NeuralGCM, recently released by Google Research, aims to provide a deep-learning based atmospheric model for both weather forecasting and climate projections. As reviewed in Section 1, hybrid models, such as NeuralGCM, integrate data-driven approaches with physical modeling, offering enhanced numerical stability and reduced model drift, thereby enabling simulations that extend beyond short-term weather forecasts (Brenowitz & Bretherton, 2019;Kochkov et al., 2024). The model architecture consists of a differentiable dynamical core (Gelbrecht et al., 2023), which solves a set of hydrostatic primitive equations with moisture for atmospheric dynamics and thermodynamics (Bourke, 1974). This dynamical core is coupled to a neural network-based physics module that parameterizes unresolved subgrid-scale processes, including radiative transfer, moist convection, precipitation, and turbulent mixing (see Figure 1 in Kochkov et al., 2024).
In this study, we utilize the deterministic configuration of NeuralGCM at horizontal resolutions of 2.8°and 1.4°, with 37 interpolated vertical levels extending from 1,000 to 1 hPa (in the dynamical core, 32 equidistant sigma levels for vertical discretization from the surface to the top of atmosphere are implemented), following the configuration described in Kochkov et al. (2024). Therefore, NeuralGCM can be categorized as the low-top atmospheric model. Importantly, NeuralGCM's encoder-decoder interface permits modification of surface boundary conditions, such as SIC and sea-surface temperature (SST). Leveraging this flexibility, we perform sea ice-forced experiments by prescribing varying SIC conditions, as detailed in Section 2.3. We show the results of NeuralGCM with 2.8°resolution in the main text and presents the results with 1.4°resolution in the Supporting Information.
For weather forecasting with a 10-day lead time, NeuralGCM demonstrates prediction skill comparable to that of purely data-driven models such as Pangu (Bi et al., 2023) and GraphCast (Lam et al., 2023). For medium-range forecasts extending up to 15 days, NeuralGCM exhibits superior predictive performance relative to the ensemble prediction products from the European Center for Medium-Range Weather Forecasts (ECMWF). In terms of seasonal prediction, NeuralGCM realistically reproduces the seasonal cycle of global mean temperature at 850 hPa, consistent with both the ERA5 reanalysis data set (Hersbach et al., 2020) and simulations from the global cloud-resolving climate model X-SHiELD (Cheng et al., 2022). On decadal time-scales, Atmospheric Model Intercomparison Project (AMIP)-style simulations for the period 1981-2014 (Eyring et al., 2016) indicate that NeuralGCM effectively captures the observed global warming trends. Moreover, NeuralGCM exhibits reduced spatial biases compared to AMIP simulations from the CMIP6 archive. These findings, as comprehensively documented in Kochkov et al. (2024), reinforce the potential of NeuralGCM as a reliable tool for a wide range of climate modeling applications. More recent studies further showed that NeuralGCM has the ability to capture the Arctic amplification (Rucker et al., 2025) and extreme heatwaves (Duan et al., 2025), and investigated the atmospheric responses to uniform SST warming (Zhang & Merlis, 2025).
To provide a fully dynamical modeling reference for the evaluation of NeuralGCM's performance, we use WACCM6 configuration of the Community Earth System Model version 2 (CESM2, Danabasoglu et al., 2020) with a horizontal resolution of 0.9°latitude × 1.25°longitude and 70 vertical layers extending from the surface to 6x10 -6 hPa (140 km) (Gettelman et al., 2019). The dynamical core in WACCM6 is based on the Godunov-type finite volume scheme in which noise-free numerical solutions can be obtained (Lin & Rood, 1997). The physical parameterization schemes of WACCM6 are the same as those in the Community Atmosphere Model version 6 (CAM6, Danabasoglu et al., 2020), except for the nonorographic gravity wave drag parameterization scheme (Gettelman et al., 2019). The updated gravity wave drag scheme improves the simulation of frontogenesis, and gives rise to a realistic self-generated quasi-biennial oscillation (QBO). Recent studies reported that WACCM6 simulates more realistic sudden warming events and QBO in the stratosphere, as well as storm tracks, stationary waves, and blockings in mid-latitude troposphere, compared with those simulated in CAM6 (Ayarzagüena et al., 2020;Gettelman et al., 2019;Simpson et al., 2020).
In addition, WACCM6 is developed to resolve a full range of chemical reactions and the interactions between chemical processes and atmospheric dynamics and thermodynamics (Gettelman et al., 2019). However, we do not activate the interactive chemistry module in WACCM6 due to its high computational costs (see Table 2 in Gettelman et al., 2019), instead we use the version with specified chemistry. In this configuration, for example, the ozone profile is prescribed; therefore, interactions between ozone and atmospheric circulation are not permitted. Thorough documentation and discussions on the effect of chemistry-atmosphere interactions can be found in Gettelman et al. (2019).
We also use the atmospheric component of Taiwan Earth System Model Version 1 (TaiESM1), developed by the Research Center for Environmental Change, Academia Sinica, to test the physical model dependence. TaiESM1 shares the same spatial resolution as WACCM6, except its vertical layer only contains 30 layers with a model top at 2 hPa. Thus, TaiESM1 represents a low-top model with respect to the high-top configuration of WACCM6. The model is based on CAM version 5.3 (Neale et al., 2010) with a finite-volume dynamical core (Lin, 2004). A newly-developed triggering functions for the deep convection scheme (Y.-C. Wang et al., 2015), a revised cloud macrophysics scheme (Shiu et al., 2021), and a three-moment aerosol scheme (Chen et al., 2013) were implemented into CAM version 5.3. The model also considers the three-dimensional radiation-topography interactions (Lee et al., 2013). More details and model performance were documented in Lee et al. (2020). The results of TaiESM1 are shown in Supporting Information S1 and discussed throughout the manuscript.
We follow the experimental protocol designed by the PAMIP (D. M. Smith et al., 2019) to perform time-slice simulations using NeuralGCM, WACCM6, and TaiESM1. Each simulation is integrated from April 1st of the initial year through 31 May of the following year, yielding 14-month simulations for analysis. For the WACCM6 and TaiESM1 simulations, the CMIP6 radiative forcing at year 2000 is fixed during the integration period. To be consistent, the NeuralGCM is initialized with conditions from ERA5 for 1 April 2000. The simulations use prescribed SIC fields that represent pre-industrial and future conditions, specified as lower boundary conditions. The pre-industrial SIC fields are derived from the Coupled Model Intercomparison Project Phase 5 (CMIP5) historical simulations, selected based on their consistency with the observed global-mean surface temperature (about 13.67°C) during the pre-industrial era (Haustein et al., 2017). Future SIC fields, in contrast, are obtained from CMIP5 projections under the Representative Concentration Pathway 8.5 (RCP8.5) scenario (Hausfather, 2020), corresponding to a climate state in which global-mean surface temperature exceeds the preindustrial level by approximately 2°C. In the Southern Hemisphere, SIC is held fixed at present-day values, estimated using the 1979-2008 mean SIC from the Hadley Center Sea Ice and Sea Surface Temperature (HadISST) data set (Rayner et al., 2003). The spatial distributions of SIC for March and September, as well as the temporal evolution of Arctic-mean monthly SIC, are illustrated in Figure 1 of Liang et al. (2024). The preindustrial and future experiments are also referred to as experiment 1.5 pdSST-piArcSIC and 1.6 pdSST-futArcSIC listed in Tabel 1 in D. M. Smith et al. (2019).
To account for the influence of internal variability, the PAMIP protocol recommends that participating models generate 100 ensemble members for each simulation set (D. M. Smith et al., 2019). However, recent studies have indicated that even larger ensemble sizes may be necessary to robustly quantify the uncertainties associated with internal variability (Labe et al., 2019;Liang et al., 2020;Peings et al., 2021;Sun et al., 2022). In light of this, we expand the ensemble size to 200 members for the NeuralGCM, WACCM6, TaiESM1 simulations. To generate ensemble members, we add small perturbations to one grid point of the temperature fields in both models when starting integration. The ensemble mean is subsequently computed across all members to effectively suppress the noise arising from internal variability, thereby enhancing the signal-to-noise ratio of the sea-ice forced response (Deser et al., 2020;Maher et al., 2021). The atmospheric response to Arctic sea-ice loss is quantified by taking the difference between the ensemble-mean fields from the future sea-ice simulation and those from the pre-industrial counterpart.
Using GPUs (NVIDIA GeForce RTX 4060 Ti), to complete one 2.8 degree NeuralGCM simulation requires about 8 min, while using CPU 224 cores on supercomputer Forerunner1, where we ran WACCM6 and TaiESM1 simulations, to complete one WACCM6 simulation takes about 8 hr. It is noted that NeuralGCM with 2.8°r esolution successfully generates 200 ensemble members without the stability issue. However, the version with 1.4°suffers from the stability problem: about one-third of the ensemble members crashes during the integration.
We use the ERA5 reanalysis data set (Hersbach et al., 2020) in 1 April of 2000 as the initial conditions to integrate NeuralGCM. We also analyze daily ERA5 500-hPa geopotential height (hereafter Z500) for observational estimates of blocking during the 1979-2023 period. A global mean cold bias in the lower stratosphere during the 2000-2006 period was documented (Simmons et al., 2020), so caution is needed when analyzing stratospheric fields.
Blocking, a key atmospheric phenomenon examined in this study, is characterized by a persistent high-pressure system that disrupts the prevailing jet stream and alters synoptic weather patterns in mid-latitudes (Arnheim et al., 2025;Booth et al., 2017;Rex, 1950;Woollings et al., 2018). Given its strong association with extreme weather events-including cold spells and heatwaves-blocking has wide socioeconomic implications (Kautz et al., 2021;Woollings et al., 2018). However, previous studies have reported inconsistent responses of blocking to Arctic warming. Some studies suggested an increase in blocking frequency (J. A. Francis & Vavrus, 2015;Liu et al., 2012;Mori et al., 2014), while others indicated minimal changes or even a decrease in blocking occurrence (Kennedy et al., 2016).
To investigate the response of blocking in the North Atlantic sector (90°W-60°E and 35°N-75°N) to Arctic seaice loss, we detect a blocking event using Z500. Specifically, we follow the blocking detection algorithm developed by Scherrer et al. (2006) to calculate the Z500 meridional gradients to the north and south of a target longitude and latitude (x 0 , y 0 ):
where y S = y 0 -16°and y N = y 0 + 16°. If the Z500 gradients satisfy two criteria, we identify it as an instantaneous blocking: (a) ∇Z500 S > 0 and (b) ∇Z500 N < -10 m/°. The first criterion characterizes the reversal of the Z500 gradient to the south, which occurs when a block emerges (Lejenäs & Økland, 1983); whereas the second criterion ensures that a westerly wind to the north of the blocking occurs as the surrounding jet stream is split (Barnes et al., 2012;Tibaldi & Molteni, 1990). Finally, if instantaneous blocking lasts for five or more consecutive days, we count those days as the blocking days.
To define the jet latitude associated with blocks, we take the zonal average of the daily 850-hPa zonal wind (hereafter U850) across different latitudes over the longitudinal domain of 75°W to 15°E, following previous
10.1029/2025MS005264
studies (Davini & Cagnazzo, 2014;Kwon et al., 2020;Woollings et al., 2010). A 5-day running average is then performed over the time dimension, and the latitude with highest wind speed is determined as the jet latitude.
We also investigate the storm track response to Arctic sea-ice loss. The storm track is defined as the transient eddy heat flux (<v ʹTʹ> ) at 850 hPa following Kwon et al. (2020). We first apply the Lanczos high-pass filter at 8 days with 15 weights to the daily meridional wind (v) and the air temperature t. The strength of the storm track is calculated as the covariance of filtered v and T during extended winter from December to March.
To enhance our understanding of the underlying causes of storm track changes, we calculate the Eady growth rate (Eady, 1949;Lindzen & Farrell, 1980;Simmonds & Lim, 2009), which characterizes the development of baroclinic instability (Kwon et al., 2020). The Eady growth rate can be formulated as:
where f is the Coriolis parameter, z the vertical coordinate, U is the wind component, N the Brunt-Väısälä frequency (N 2 = g θ ∂θ ∂z , g is the acceleration due to gravity and θ is the potential temperature). We multiply the rate by 86,400 (second per day) to convert to units of 1/day.
We briefly investigate the climate mean states and extremes in the three models, as the mean state differences could influence the responses to Arctic sea-ice loss. We first look into the 60°N zonal-mean zonal wind at 10 hPa (U10) and the number of sudden stratospheric warmings (SSWs), which are identified following the algorithm from Charlton and Polvani (2007). NeuralGCM has higher U10 during October-January and similar U10 compared with WACCM6, while TaiESM1 simulates stronger U10 throughout the analysis period (Figure S1a in Supporting Information S1). Correspondingly, NeuralGCM produces similar SSWs as WACCM6 in most months, except November, while TaiESM1 underestimates the occurrence of SSWs (Figure S1b in Supporting Information S1). These results suggest NeuralGCM is as good at capturing the strength of mean stratospheric circulation and extreme events, as WACCM6.
We next analyze the zonal-mean zonal wind and temperature in height-latitude coordinates. In the troposphere, the three models generally agree with each other (Figure S2 in Supporting Information S1). Evident difference occurs in the stratosphere where the temperatures between equator and 30°N, between 10 and 1 hPa in Neu-ralGCM and WACCM6 are stronger and extend further to the north than TaiESM1 (Figures S3a-S3c in Supporting Information S1). The stratospheric zonal winds in mid-latitudes are strongest in TaiESM1, followed by WACCM6 and NeuralGCM (Figures S3e and S3f in Supporting Information S1).
To further examine the longitude-latitude spatial structure of tropospheric and stratospheric dynamical and thermodynamical variables, we select near-surface air temperature at 1,000, 500, 10 hPa (hereafter referred to as T1000, T500, T10), and geopotential height at 1,000, 500, 10 hPa (hereafter referred to as Z1000, Z500, Z10). The spatial patterns of temperature are similar in the three models, but the magnitude of temperature in NeuralGCM is overall weaker than those in the other two models, particularly at 10 and 1,000 hPa (Figure S4 in Supporting Information S1). The spatial structures of geopotential height are similar in the three models, for example, the low height values above Barents-Kara Seas at 10 hPa, but NeuralGCM has overall lower values at 10 hPa and higher values at 500 and 1,000 hPa (Figure S5 in Supporting Information S1).
To investigate extreme events simulated in the three models, We follow Johnson et al. (2018) to identify cold and warm extremes. Without the hourly data output, we instead use daily data and consider the temperature anomalies falling and above below the tenth percentile during consecutive 5 days. We show in Figure S6 in Supporting Information S1 the spatial distributions of cold and warm extreme frequency during the December-January-February (DJF) period. For the warm extremes, the three models present similar large-scale patterns. Neu-ralGCM somehow overestimates the low frequency regions, including the western coasts of the North America, North Europe, and north of Iceland. For the cold extremes, NeuralGCM generally underestimates the extreme frequency. Therefore, NeuralGCM may find it challenging to capture cold extremes.
10.1029/2025MS005264
We begin by examining the prescribed Arctic sea-ice area (SIA) in both the pre-industrial and future simulations. As shown in Figure 1a, the seasonal cycle of SIA in both scenarios displays a minimum in September and a maximum in March. As designed, the SIA in the future scenario is lower than that in the pre-industrial scenario throughout the year, with the largest difference (approximately 9.3 × 10 6 km 2 ) in October (gray bars in Figure 1a).
In NeuralGCM, WACCM6, and TaiESM1, the Arctic T1000 exhibits a delayed seasonal cycle, with minimum values occurring in January or February and maximum values in July (Figure 1b). This phase lag between SIA and T1000 reflects the influence of ocean-atmosphere heat exchange and associated local feedback processes, as highlighted in previous studies (Bintanja & Van der Linden, 2013;Boeke & Taylor, 2018;Chung et al., 2021;Hahn et al., 2022;Liang et al., 2022;D. M. Smith et al., 2017;Wu et al., 2023). The T1000 values simulated by NeuralGCM are higher during the cold season compared to those simulated by WACCM6 and TaiESM1, and a reduction in the amplitude of the near-surface air temperature seasonal cycle in NeuralGCM. We also notice that the ensemble spread, estimated as one standard deviation across 200 members (color shadings in Figure 1b), in the three models are similar, suggesting that NeuralGCM produces internal variability comparable to WACCM6 and TaiESM1.
Focusing on the atmospheric response to Arctic sea-ice loss, the T1000 response simulated by NeuralGCM begins to increase during the warm season, peaks in November-approximately one month after the maximum sea-ice forcing shown in October-and subsequently declines into the cold season (Figure 1c). This lagged response indicates that NeuralGCM effectively captures the delayed atmospheric adjustment to Arctic sea-ice loss, consistent with findings from previous studies (Liang et al., 2021;Peings & Magnusdottir, 2014). A similar seasonal evolution in T1000 responses are also simulated by WACCM6 and TaiESM1. The T1000 responses in WACCM6 and TaiESM1 are slightly stronger than that of NeuralGCM after November. In terms of internal variability of the temperature response, the ensemble spreads tend to increase during middle to late winter (post-January), indicating enhanced atmospheric variability in the deep cold season in both models. The ensemble spread range of NeuralGCM is mostly within the range of WACCM6 and TaiESM1, with the former smaller than the latter.
We next examine the relationship between Arctic and global T1000 responses to assess the extent of AA in NeuralGCM, WACCM6, and TaiESM1 (Figure 1d). In all models, the Arctic T1000 response is consistently stronger than the global-mean response, presenting AA. While NeuralGCM simulates a slightly smaller Arctic T1000 response and a similar global response compared to WACCM6 and TaiESM1, the magnitude of ensemble spread is overall similar. This result confirms again that the effect of internal variability is similar in the three models.
The resolution of the model has been discussed to affect the simulated atmosphere and ocean (Chang et al., 2020;Scaife et al., 2019). To test the sensitivity of the response to model spatial resolution, we conduct the same NeuralGCM experiments with a higher horizontal resolution (1.4°). The seasonal evolution remains largely consistent, with the response magnitude in NeuralGCM weakened and the variability more comparable to WACCM6's (Figure S7c in Supporting Information S1). Notably, the higher-resolution version of NeuralGCM enlarges the internal variability of the global T1000 response (red dots in Figure 1d). However, the Arctic and global responses become somewhat weaker. We further discuss the impact of resolution in the Discussion section.
We next look into the zonal-mean temperature and zonal wind responses simulated by NeuralGCM. Two prominent warming regions emerge in the mid-to high-latitude troposphere and stratosphere (Figure 2a). The first likely corresponds to the direct thermal response to Arctic sea-ice loss, characterized by a strong near-surface warming whose amplitude decreases both southward and upward. A similar bottom-heavy warming structure is also evident in the WACCM6 simulation (Figure 2c). The associated zonal-mean zonal wind responses in both models exhibit a weakening centered around 60°N (Figures 3a and 3c), broadly consistent with the thermal-wind balance adjustment arising from hydrostatic and geostrophic balance of forces. The second warming region in NeuralGCM appears in the high-latitude stratosphere, centered near 50 hPa (Figure 2b). In contrast, WACCM6 shows a different response, with stratospheric cooling north of 60°N and warming to the south, resulting in opposing stratospheric wind responses between the two models (Figures 3b and 3d). We note that the zonal-mean temperature and wind responses in the higher-resolution version of NeuralGCM are generally stronger than those in the lower-resolution counterpart (Figures S8 and S10 in Supporting Information S1). We also find that TaiESM1 simulates a warming and weakening of zonal wind in the stratosphere (Figures S9 and S11 in Supporting Information S1), but the weakening center is different from NeuralGCM's. The tropospheric responses in TaiESM1 are similar to those in NeuralGCM. Taken together, these results suggest that NeuralGCM produces a reasonable simulation of the large-scale atmospheric temperature and circulation responses to Arctic sea-ice loss in the troposphere when compared with WACCM6 and TaiESM1. However, its representation of stratospheric responses appears to deviate from that of WACCM6 and TaiESM1, highlighting potential limitations in simulating stratospheric processes.
In addition to the thermal wind adjustment, the propagation of Rossby waves and the associated eddy-mean flow interactions also modulate the tropospheric and stratospheric jets (Andrews et al., 1987). To assess this mechanism, we calculate the Eliassen-Palm (EP) flux and its divergence (proportional to positive zonal-wind acceleration) (Edmon et al., 1980). In both models, the EP flux responses predominantly propagate upward from the surface between 40°N and 50°N, turn poleward near 300 hPa, and descend around 60°N (Figures 4a and 4c). In NeuralGCM, this downward propagation is particularly pronounced, resulting in a positive zonal wind tendency due to EP flux divergence between 800 and 500 hPa over the middle-to-high latitudes. In contrast, WACCM6 exhibits an upward EP flux response in the lower atmosphere that counteracts the descending flux, leading to EP flux convergence and a negative zonal wind tendency. This likely contributes to the stronger negative zonal wind response in WACCM6 compared to NeuralGCM (Figures 3c and 3a).
In the upper troposphere and stratosphere, EP flux responses display upward behaviors in the two models, but the opposite wind tendencies (Figures 4b and 4d): in NeuralGCM a negative wind tendency occurs above 50 hPa to the north of 60°but there is a positive tendency in WACCM6. To investigate how the opposite wind tendency occurs, we calculate the total EP fluxes entering the high-latitude stratosphere region (the blue box in Figures 4b and 4d). We find that large EP fluxes leave the southern boundary of the selected region in NeuralGCM (the blue number, which is calculated following Sigmond and Scinocca (2010)), favoring an enhanced zonal wind tendency. But this is not seen in the wind tendency. In WACCM6, the net loss of EP fluxes, although smaller than that in NeuralGCM, corresponds to enhanced zonal-wind tendency. These results indicate the EP fluxes may not be the main process in NeuralGCM to modulate the zonal wind, or the eddy-mean flow interaction is not well represented. In addition, the high-resolution simulation of NeuralGCM shows opposite EP fluxes and zonal wind tendency responses with low-resolution version in the stratosphere (Figures S12a and S12b in Supporting Information S1). TaiESM1 also shows different EP flux responses from NeuralGCM and WACCM6 in the Journal of Advances in Modeling Earth Systems 10.1029/2025MS005264
stratosphere, but the enhanced wind tendency locates at levels higher than 5 hPa (Figures S13a and S13b in Supporting Information S1). These discrepancies, again, highlight the limitations of stratospheric representation in NeuralGCM.
To investigate stratosphere-troposphere interactions in response to Arctic sea-ice forcing during the cold season in the Northern Hemisphere, we analyze the polar-cap (65°N-90°N) averaged geopotential height (hereafter Z) responses from 1 September to 31 March of the following year. NeuralGCM exhibits strong vertical coupling between the troposphere and the entire stratosphere, characterized by persistent positive height anomalies extending from November through March, except for a period of weak anomalies from the middle of December to the middle of January, (Figure 5a). In contrast, WACCM6 displays weaker stratosphere-troposphere coupling, with positive height responses largely confined to the lower stratosphere during this period (Figure 5b). As a result, notable discrepancies between the two models emerge above 50 hPa, with the largest differences occurring in October, November, February, and March (Figure 5c), a pattern also shown in the high-resolution version of NeuralGCM (Figure S14c in Supporting Information S1). Similarly, TaiESM1 presents a weaker stratospheretroposphere coupling (Figure S15 in Supporting Information S1), compared with NeuralGCM. These Journal of Advances in Modeling Earth Systems The zonal-mean and polar-cap averages offer an informative dynamical perspective on the large-scale atmospheric circulation responses to Arctic sea-ice loss. However, zonally asymmetric features also deserve examination, as they reveal horizontal patterns that may relate to specific weather regimes and regional weather phenomena. We look at the geographic distribution of T1000 responses in the two models. The warming amplitude of NeuralGCM is overall weaker than that of WACCM6. NeuralGCM captures the warming hotspots associated with large sea-ice loss in the Barents-Kara Seas, the Bering Strait, the Sea of Okhotsk, and the Hudson Bay (Figure 6c), which are also shown in WACCM6 (Figure 6f), although the temperature response is generally weaker in NeuralGCM than in WACCM6. At 500 hPa, a dipolar temperature response with warming above Greenland and North America, and cooling over eastern Asia and North Pacific appears in both NeuralGCM and WACCM6, although the former has some sporadic cooling to the north of Canadian archipelago and less cooling in the North Pacific (Figures 6d and 6e). Moving upward to stratosphere, the discrepancy between the two models becomes evident. NeuralGCM exhibits a warming response over Siberia, while WACCM6 shows dominant Journal of Advances in Modeling Earth Systems 6a and 6d). In TaiESM1, tropospheric temperature responses are similar, except the warming over land is weaker (Figure S17f in Supporting Information S1). The stratospheric responses, on the other hand, present a dipole pattern with positive center in Eurasia and negative one in North Alantic, distinct from NeuralGCM and WACCM6 (Figure S17d in Supporting Information S1).
Similar to the temperature responses, the geopotential height responses in NeuralGCM and WACCM6 exhibit similar spatial distributions, except in the stratosphere. At 1,000 hPa, NeuralGCM effectively captures the characteristic heat-low circulation response to local warming induced by sea-ice loss (Deser et al., 2010;Peings & Magnusdottir, 2014), along with the elevated height extending from Greenland and the North Atlantic to Siberia (Figure 7c), consistent with the WACCM6 simulations (Figure 7f). Notably, in the North Atlantic sector, the height response resembles a negative NAO pattern. At 500 hPa, both models show an overall increase in geopotential height, with the strongest anomalies centered over the Barents-Kara Seas, where sea-ice loss is most pronounced. NeuralGCM simulates a slightly weaker Z500 response over the Barents-Kara Seas but exhibits a pronounced height increase over eastern Siberia (Figure 7b), a feature absent in WACCM6 (Figure 7e). Corresponding to these geopotential height responses, a general weakening of zonal wind is presented over midlatitudes in the lower and middle troposphere (Figure 8). Compared to WACCM6, NeuralGCM produces a stronger zonal wind weakening over the North Atlantic but a weaker response over the North Pacific (Figures 8b, In panels (a) and (b), the contour line indicates that the response is statistically significant with a twotailed Student's t test at 95% confidence level, while that in (c) presents the difference of the responses between the two models is statistically significant.
10.1029/2025MS005264 8c, 8e, and 8f). As expected, the geopotential height and zonal wind responses in the stratosphere differ substantially between NeuralGCM, WACCM6, and TaiESM1 (Figures 7a, 7d, 8a, and 8d and Figures S19a, S19d, S21a, and S21d in Supporting Information S1).
In summary, NeuralGCM effectively simulates key atmospheric responses to Arctic sea-ice loss in the troposphere, broadly consistent with WACCM6 and TaiESM1 simulations. These include localized warming to large
10.1029/2025MS005264 sea-ice decline, the heat-low circulation patterns, and a hemispheric negative NAM response or a negative NAO response in the North Atlantic, accompanied by a general weakening of mid-latitude jets. However, substantial discrepancies emerge in the stratospheric responses between the (the three models). The NeuralGCM may be improved in representing stratospheric dynamical and thermodynamical processes, although a careful assessment of its response within the inter-model range or not is needed. Overall, NeuralGCM serves as a valuable tool for investigating Arctic-midlatitude linkages through tropospheric pathways.
We now examine the response of atmospheric blocking in the North Atlantic sector to Arctic sea-ice loss. We first analyze the DJF ensemble-mean number of blocking days and the zonal wind at 200 hPa simulated by Neu-ralGCM under future and pre-industrial sea-ice conditions (Figures 9a and 9b). Blocking events predominantly occur along the periphery of the prevailing 200-hPa jet, forming an arc-shaped pattern with centers over Greenland, the British Isles, southern Norway, and the eastern Atlantic. Under future sea-ice conditions, blocking frequency increases north of this arc while decreasing to the south, relative to the pre-industrial simulations (Figures 9a and 9b). The blocking changes reflect the deceleration to the north and acceleration to the south of the North Atlantic jet (dashed contour lines in Figure 9c). In WACCM6 and TaiESM1, the spatial distribution of blocking events exhibits a similar arc pattern, but with the high-frequency region (Figures 9d and 9e and Figures S23d and S23e in Supporting Information S1) is different from NeuralGCM's, a well-known bias in AGCMs (Kwon et al., 2018). WACCM6 and TaiESM1 show comparable blocking and zonal wind responses, but their magnitudes are notably weaker (Figure 9f and Figure S23f in Supporting Information S1). In terms of the spatial pattern and magnitude of blocking days in ERA5, it is clear that NeuralGCM simulates more comparable features with ERA5 than WACCM6 and TaiESM1 (cf., Figures 9a and 9d, Figures S23 and S28 in Supporting Information S1). In particular, NeuralGCM simulates more frequent blockings above the British Isles and above southern Greenland, which are closer to those in ERA5. This result suggests that NeuralGCM, even using coarser spatial resolution, may be more reasonably capturing the spatiotemporal structure of extremes than traditional physics-based climate models. It is encouraging to see that the NeuralGCM exhibits much more realistic blocking pattern than typical AGCMs.
To investigate the mechanisms driving the blocking response to sea-ice loss, we examine the meridional temperature gradient (-dT /dy), the Eady growth rate, and the transient eddy heat flux (<v ʹTʹ> ) at 850 hPa simulated by NeuralGCM (Figure 10). A weakened meridional temperature gradient extends from northeastern Canada to Iceland and the Barents-Kara Seas (Figure 10a), indicating a reduction in baroclinic instability under sea-ice forcing. These regions largely coincide with or are adjacent to areas of lower tropospheric warming induced by regional Arctic sea-ice loss. Consistently, the Eady growth rate (Figure 10b) is also reduced in these regions, leading to a weakening of transient eddy heat fluxes (Figure 10c) and suggesting suppression of storm track activity. WACCM6 and TaiESM1 exhibit similar responses (Figures 10d,10e,and 10f,and Figures S25d,S25e,and S25f in Supporting Information S1), indicating robust dynamical consistency across these models. Journal of Advances in Modeling Earth Systems 10.1029/2025MS005264
To further assess the response of the eddy-driven jet, we analyze the 850-hPa jet latitude in NeuralGCM simulations. The jet latitude of the ensemble-mean wind fields during DJF exhibits a southward shift in the future seaice simulations compared to the pre-industrial simulations, as indicated by the shifted histogram (Figure 11a). The DJF-averaged ensemble-mean jet latitude in the future simulations is located at 47.64°N (red vertical line in Figure 11a), whereas in the pre-industrial simulations, it is positioned at 48.68°N (blue vertical line in Figure 11a). In contrast, WACCM6 simulates jet latitudes that are generally positioned farther north, with less latitudinal variability and a smaller difference between future and pre-industrial simulations (Figure 11c). These findings suggest that NeuralGCM produces a stronger jet latitude response than WACCM6, consistent with its more pronounced blocking and stronger T responses (Figures 9c and 9f). We further examine the frequency distribution of jet latitude occurrence across DJF days and 200 ensemble members in both models (Figures 11b and 11d). Both future and pre-industrial simulations exhibit two distinct peaks in jet latitude, consistent with reanalysis data sets (Kwon et al., 2018(Kwon et al., , 2020)). Under future sea-ice conditions, the southern peak becomes more pronounced, reflecting the southward shift of the ensemble-mean jet latitude. Additionally, the separation between the two peaks is more distinct in NeuralGCM than in WACCM6, highlighting differences in the ability of these models to simulate the eddy-driven jet in the North Atlantic sector. Also, the more realistic southerly jet position in Neu-ralGCM is consistent with more realistic number of blocking days over Greenland, British Isle and the North Sea. In TaiESM1, the jet position is biased northwards, but the difference between future and pre-industrial simulations is closer to NeuralGCM (Figure S27c in Supporting Information S1).
In summary, the above results reveal a mechanistic link between Arctic sea-ice loss and blocking responses in the North Atlantic sector can be simulated by NeuralGCM. The sea-ice loss, particularly in Hudson Bay and the periphery of Greenland, induces localized warming, which weakens the meridional temperature gradient and reduces baroclinicity from northeastern North America to the Barents-Kara Seas. This weakening of baroclinicity suppresses storm track activity, as indicated by a reduction in synoptic transient eddy heat flux, leading to a southward shift of the eddy-driven jet through eddy-mean flow interactions (Novak et al., 2015;Shaw, 2019). The resulting atmospheric circulation regime favors Rossby wave breaking on the poleward flank of the jet (Barnes & Hartmann, 2010;Davini & Cagnazzo, 2014;Kwon et al., 2018;Rivière, 2009;Woollings et al., 2008), increasing the frequency of blocking events over Greenland and Iceland while decreasing blocking occurrences along the western European coast. All these processes are well-represented in NeuralGCM that give rise to more realistic blockings.
In this study, we utilize NeuralGCM, a hybrid ML atmospheric model, to investigate atmospheric responses to Arctic sea-ice loss for the first time. Following the PAMIP protocol, we conduct simulations forced by preindustrial and future sea-ice conditions, generating 200 ensemble members. Although NeuralGCM simulates an overall warmer near-surface temperature than WACCM6, a fully dynamical climate model, the internal variability (estimated by the spread across 200 ensemble members) in both models are comparable. NeuralGCM simulates Arctic near-surface warming by 2.64 K during boreal winter in response to sea-ice loss, aligning with the warming by 2.90 K produced by WACCM6. In line with the reduced meridional temperature gradient, the tropospheric jet stream weakens in both models. The spatial distribution of the geopotential height response presents a localized heat-low circulation in the lower troposphere, corresponding to regional Arctic warming hotspots, and a negative NAO-like pattern in the North Atlantic sector in the lower and middle tropospheres, and a weakened midlatitude jet stream. These findings demonstrate the capability of NeuralGCM as a viable hybrid atmospheric model for studying the tropospheric impacts of Arctic climate change. Journal of Advances in Modeling Earth Systems 10.1029/2025MS005264
Substantial discrepancies emerge between NeuralGCM and WACCM6 in the stratosphere, however. In high latitudes above 50 hPa, the temperature and the associated circulation responses exhibit opposite patterns in the two models. Analysis of EP flux divergence further reveals that it also differs in sign, due to more meridional fluxes leaving the high-latitude stratosphere in NeuralGCM, contributing to the differences in zonal wind response. Additionally, these inconsistencies in stratospheric responses lead to distinct stratosphere-troposphere coupling characteristics, with NeuralGCM simulating an overly strong interaction between the two atmospheric layers. These pronounced differences underscore the need for caution when applying NeuralGCM to study atmospheric responses to Arctic sea-ice loss via stratospheric pathways and highlight the necessity of improving its stratospheric representation.
A key question immediately arises: what factors contribute to the discrepancies between NeuralGCM and WACCM6? To address this question, we first examine the influence of spatial and vertical resolution, as Neu-ralGCM operates at a horizontal resolution of 2.8°in both latitude and longitude with a low-top vertical configuration, whereas WACCM6 employs a higher resolution of 0.9°latitude × 1.25°longitude and a high-top configuration. Previous studies have suggested that model resolution significantly affects stratospheric representation and extratropical atmospheric variability (e.g., Serva et al., 2024). To investigate this effect, we conduct parallel NeuralGCM simulations with an increased horizontal resolution of 1.4°in both latitude and longitude. This refinement leads to a decrease in the Arctic and global ensemble-mean T1000 responses and an increase in the internal variability of the global T1000 response (red dots in Figure 1d), making them more comparable to the WACCM6 results. The effect of internal variability on other dynamical variables require further examinations. In the troposphere, the overall response weakens slightly (Figures S16, S18, and S20 in Supporting Information S1). Conversely, in the stratosphere, the zonal-mean atmospheric responses become more pronounced (Figures S8 and S10 in Supporting Information S1), accompanied by enhanced stratosphere-troposphere coupling (Figure S14 in Supporting Information S1). We hypothesize that the differences could stem from (a) the low-top configuration of NeuralGCM dynamical core, (b) the absence of a gravity-wave parameterization scheme, (c) differences in basic mean state, and (d) the regionality of sea-ice forcing.
To assess the first factor, we conduct similar experiments using TaiESM1 with a low-top configuration. Similar and consistent tropospheric responses are simulated by TaiESM1 (Figures S17f, S19f, and S21f in Supporting Information S1), but the stratospheric responses are still inconsistent with NeuralGCM (Figures S9, S11, and S15 in Supporting Information S1). Regarding the second factor, incorporating a neural network or ML-based gravitywave parameterization scheme (e.g., Espinosa et al., 2022;Hardiman et al., 2023) or a well-trained stratosphere in future generations of NeuralGCM could help address this limitation. Indeed, the second factor can also be applied to any other parameterization scheme absent in NeuralGCM. For the third factor, comparing results from Neu-ralGCM to those from WACCM6 should be considered with caution, as WACCM6 could be biased in simulating the mean state and atmospheric responses to sea-ice loss. For example, the zonal-mean structure of temperature persistence in WACCM6 is smaller than the other models (Lewis, Seviour, et al., 2024). To address this issue, we also look into the results from TaiESM1. Indeed, the difference between the mean states and responses in the stratosphere in WACCM6 and TaiESM1 can be similar to or larger than the difference between these variables in WACCM6 and NeuralGCM. The tropospheric responses, on the other hand, are mostly consistent across these models. This again illustrates that the stratospheric responses to Arctic sea-ice loss are distinct in the three models. In future study, it would be informative and useful to use PAMIP simulations from different GCMs to provide a range of model structural uncertainty. For the last factor, we notice that the surface warming in the Pacific sector is stronger in WACCM6 than in NeuralGCM (Figures 6c and 6f), which could potentially lead to a stronger stratospheric polar vortex. Sun et al. (2015) showed this connection between warming in the Pacific side and the stratospheric polar vortex using previous version of CESM. An extended study could use the regional sea-ice experiments designed in PAMIP to address this factor more carefully.
In NeuralGCM simulations, we find that wintertime blocking frequency increases in the northern branch of North Atlantic while decreasing in the southern part of the basin. This pattern is driven by Arctic sea-ice loss, which reduces baroclinic instability and weakens the eddy-driven jet. Additionally, a positive feedback mechanism between blocking and enhanced Rossby wave breaking may further reinforce these responses (Michel & Rivière, 2011). Notably, both the blocking response and the eddy-driven jet shift are stronger in NeuralGCM than in WACCM6 and in TaiESM1, and increasing the horizontal resolution in NeuralGCM weakens both the blocking response and the latitudinal shift of the jet (Figures S22-S25 in Supporting Information S1). One of the most promising findings from NeuralGCM is that its simulated blocking frequency more closely resembles that of Journal of Advances in Modeling Earth Systems 10.1029/2025MS005264 ERA5 data, compared to WACCM6 andTaiESM1 (cf. Figures 9a, 9d, Figures S22, S23, andS28 in Supporting Information S1). This result suggests that NeuralGCM may be more capable of capturing highly impactful weather regimes, whereas WACCM6 and TaiESM1 might underestimate the these weather patterns and their responses to sea-ice loss. However, further rigorous analysis of other AI-based and traditional climate models is required to substantiate this conclusion.
Our results demonstrate that NeuralGCM is largely capable of simulating tropospheric responses to Arctic sea-ice loss. In the current NeuralGCM configuration, SST and SIC are prescribed as lower boundary conditions, thermodynamically forcing the atmosphere. However, previous studies have suggested that the absence of threeway coupling between sea ice, the ocean, and the atmosphere in AMIP-style simulations could influence atmospheric responses to Arctic sea-ice loss (see Figure 1 in Screen & Blackport, 2019). This highlights the importance of investigating the role of atmosphere-ocean-sea ice coupling in Arctic-midlatitude linkages. The PAMIP framework includes coupled atmosphere-ocean simulations with sea ice constrained to pre-industrial and future states, similar to those used in this study. However, NeuralGCM does not yet support a fully coupled configuration. We are also aware that prescribed SIC could cause the spurious warming documented in recent literature (Lewis, England, et al., 2024). Caution is needed for this spurious warming in the sea-ice experiments using NeuralGCM in the future. Recent advances in ML Earth system models have enabled the development of data-driven coupled atmosphere-ocean configurations that can simulate climate dynamics and variability on seasonal to centennial timescales (Cresswell-Clay et al., 2025;C. Wang et al., 2024). Future sea-ice experiments using these coupled AI-based models hold promise for advancing our understanding of atmosphere-ocean-sea ice interactions and their role in Arctic-midlatitude linkages. Such efforts could provide valuable contributions to the PAMIP research community and stimulate novel experiment designs.
19422466, 2025, 11, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2025MS005264 by Yu-Chiao Liang -National Taiwan University , Wiley Online Library on [11/11/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License