The change in the ITCZ location caused by CO2, which is the main greenhouse gas produced by human activity and one of the key drivers of past and present climate change, is of particular interest. Previous CO2 doubling or quadrupling experiments using an atmosphere general circulation model in conjunction with an aquaplanet slab-ocean essentially show a northward ITCZ shift because of radiative feedback from clouds and water vapor [12][13][14] . However, when the realistic distributions of continents and sea ice are taken into account, models project that the ITCZ could move either northward or southward in response to an increase in CO2 7 . The large uncertainty of ITCZ location change is mostly related to the uncertainty in the responses of clouds and sea ice. This ITCZ uncertainty persists even after ocean dynamics are included in models, where a dynamic ocean may mediate the extratropical influences on the ITCZ through changes in ocean heat transport 15,16 . For instance, the Coupled Model Intercomparison Project phase 5 (CMIP5) climate models show diverse ITCZ responses to a quadrupling of atmospheric CO2 concentration 17,18 .
Notably, the aforementioned ITCZ location changes in the CMIP5 models are based on a simulation of CO2 quadrupling primarily over about one and a half centuries. Beyond this time frame, there is a gap in our knowledge of the long-term, toward the millennial evolution of the ITCZ, or in line with climate sensitivity, the equilibrium ITCZ response to CO2 radiative forcing. This gap will be made even more clear by the fact that it will take thousands of years for the Earth system to return to equilibrium following a CO2 perturbation 19,20 . On the other hand, elucidating the equilibrium ITCZ response to increasing CO2 will help us better understand hydrological changes in future centuries 21 and past warm climates such as the warm Miocene and Pliocene Epochs 22 and Early Eocene 23 . Therefore, bridging the aforementioned gap using climate model simulations serves as the focus of the current study. Following that, we will show distinct transient and (quasi-)equilibrium responses of the ITCZ to rising CO2 as simulated by a broad range of fully coupled climate models.
We begin by examining the changes in the ITCZ location from the perspective of tropical precipitation centroid in the CMIP5 and CMIP6 CO2 quadrupling experiments, in which the atmospheric CO2 concertation in the model is abruptly increased from the preindustrial level to four times that level (Method). In comparison to preindustrial times, the multi-model mean exhibits a rapid northward shift of the annual and zonal mean ITCZ over the first one to two decades of CO2 increases (Fig. 1a), which is consistent with previous results [12][13][14] . Seen from tropical precipitation centroid (Method), the zonal-mean ITCZ moves northward by 0.18 ± 0.21 degree (multi-model mean ± one standard derivation among models) during the first 20 years.
The displacement of the rain belt to the north is particularly pronounced over the Indian Ocean, with anomalous decreases and increases in the rainfall maximum and to the north, respectively (Fig. 2a). The ITCZ deepens and narrows over the Pacific [24][25][26] , which is likely due to a strengthening of the Hadley circulation manifested as a <deep-tropics squeeze= 27 . This fast precipitation response in the first 20 years reflects some characteristics of a rapid adjustment of the climate system to abrupt CO2 forcing 28,29 as reported by the Precipitation Driver and Response Model Intercomparison Project (PDRMIP) 29 . For example, precipitation decreases over Central America, the eastern North Pacific, the Caribbean Sea, northern South America, the equatorial South Atlantic and Indian Oceans but increases over the tropical region of Africa and northern Australia. However, our first 20-year precipitation response is based on coupled model simulations, which differs from the PDRMIP fixed sea surface temperature (SST) experiment. As a result, the precipitation response also includes slower SST-mediated changes 29 , such as increased precipitation over the equatorial Pacific. After 20 years, the ITCZ starts to move south.
By the end of the 150-year CMIP5 and CMIP6 simulations, it is relatively close to its preindustrial location (Fig. 1a). We also observe a high level of model uncertainty in the CO2induced change in the ITCZ, which is consistent with previous findings 17,18 .
To elucidate the ITCZ evolution over a period longer than one century or two, we investigate CO2 quadrupling simulations by an eight-model ensemble including seven climate models from the original LongRunMIP 30 with simulation lengths of at least 1000 years and one CMIP6 model with simulation length of 1000 years (referred to as LongRunMIP for the convenience of discussion, Method). The ensemble mean of the LongRunMIP shows a strong northward ITCZ shift during the first few decades (Fig. 1c), particularly over the Indian Ocean (Fig. 2b), which is consistent with previous CMIP5 and CMIP6 model results (Fig. 1a). Tropical precipitation centroid suggests a northward ITCZ migration of 0.16 ± 0.12 degrees (multi-model mean ± one standard derivation among models) in the first 20 years. After that, the ITCZ shows a trend of southward migration, especially 100 years after the CO2 increase, at a rate of about 0.02 degrees per century between 100 and 1000 years (Fig. 1c). The southward ITCZ migration is robust over both the Atlantic and Pacific Oceans (Fig. 2c and d).
The atmospheric energy-flux theory 7,[31][32][33] , which connects the zonal-mean ITCZ location to atmospheric cross-equatorial energy transport and the net energy input to the atmosphere at the equator (Method), can help understand the non-monotonic zonal-mean ITCZ movement. We apply the atmospheric energy-flux theory to the CO2 quadrupling experiments, with a particular emphasis on the LongRunMIP and the multi-model mean result. To indicate the location of the zonal-mean ITCZ, we calculate the latitude of the energy flux equator that is determined by atmospheric cross-equatorial energy transport and the net energy input to the atmosphere at the equator (Method). We find that both metrics, the energy flux equator and tropical precipitation centroid, show a generally consistent pattern of ITCZ movement (Fig. 1). For the first two decades, the energy flux equator moves northward by 0.69 degrees (Fig. 1d) primarily related to an anomalous southward atmospheric cross-equatorial energy transport (Fig. 1c) generated by rising CO2 relative to preindustrial times, given that the contribution of the net energy input change to ITCZ movement (a southward shift by 0.07 degrees) is about one order smaller (Fig. 1d). The anomalous southward energy transport is caused by inter-hemispheric asymmetry of top of atmosphere (TOA) radiation and surface energy (Fig. 3a, Fig. 4a). CO2 increases bring about dramatic global changes in the TOA radiation feedback (Method) of water vapor, temperature, albedo, and clouds (Extended Data Fig. 1), with these changes offset between hemispheres and individual feedback (Fig. 3b). Water vapor, cloud, and albedo feedback, in particular, contributes the most to the inter-hemispheric asymmetry and results in a slightly less TOA radiation increase in the Southern than Northern Hemisphere (Fig. 3a, Fig. 4a), which is consistent with the results from previous studies [12][13][14]34 . On the other hand, the rising CO2 causes ocean heat uptake in global oceans, particularly where the ocean mixed layer is deep (Fig. 5b). The net change in inter-hemispheric surface energy asymmetry indicates that the Southern Ocean absorbs more heat than the northern ones. Note that the CO2 quadrupling simulations from CMIP5 and CMIP6 produce a consistent result (Figs. 1b and 3, Extended Data Fig. 2) on changes in cross-equatorial energy transport and inter-hemispheric asymmetry of TOA radiation, with the exception that the majority of these models prefer more heat uptake in the northern oceans (Fig. 3a).
After the first two decades, the CO2-induced southward atmospheric energy transport diminishes and even shifts northward, which is anti-correlated with the southward migration of the ITCZ 35- 37 (Fig. 1c and d). The energy flux equator and tropical precipitation centroid show significant southward migration trends between years 100 and 1000, of 0.13 degrees per century (p<0.01) and 0.02 degrees per century (p<0.01), respectively (Method). Herein we depict how atmospheric cross-equatorial energy transport varies in the CO2 quadrupling experiment during the first two decades and the final 1000 years. In the latter period, we find a northward energy transport anomaly relative to preindustrial times, which is primarily caused by an increased interhemispheric asymmetry of surface energy flux. Compared to the first two decades, the subpolar North Atlantic absorbs more heat from the atmosphere while the Arctic and North Pacific take less heat by the end of 1000 years (Fig. 5b, d and f), resulting in a smaller reduction of ocean heat uptake in the northern oceans (Fig. 4). In comparison to the northern oceans, the Southern Ocean absorbs even less heat from the atmosphere (Fig. 5b, d and f), especially between 40 o S and 60 o S (Fig. 4), which leads to an anomalous inter-hemispheric asymmetry of ocean heat uptake4the atmosphere losing less heat in the Southern than Northern Hemisphere4and thus an anomalous northward atmospheric cross-equatorial energy transport by the end of 1000 years (Fig. 3a). Note that changes in surface turbulent (sensible and latent) heat flux are primarily responsible for changes in ocean heat uptake over the Southern Ocean (Fig. 5f and h, Extended Data Fig. 3). Compared to the first two decades, despite surface warming enhances over global oceans, the delayed surface warming is especially strong over the Southern Ocean (Fig. 6) due to deep vertical mixing of water 38 and wind-driven upwelling of water from depth 39 . The delayed Southern Ocean warming 20 reduces downward turbulent heat flux at the ocean surface because of a negative turbulent heat flux feedback 37,40 , indicating that surface heat flux response acts to dampen SST anomalies.
However, compared to the inter-hemispheric asymmetry of ocean heat uptake, changes in interhemispheric asymmetry of TOA radiation have a much smaller effect on the anomalous transport of atmosphere energy (Fig. 3a). This could be due to fact that atmospheric processes modulate the atmospheric energy budget more quickly than ocean processes. Relative to the first two decades, TOA radiation has decreased globally, with the exception of a few areas such as the central and eastern tropical Pacific by the end of 1000 years (Fig. 5a, c and e). The Southern Hemisphere experiences a similar TOA radiation reduction to the Northern Hemisphere (Fig. 4c).
In contrast to the Northern Hemisphere, the water vapor and cloud feedback processes result in anomalous positive radiation entering the Southern Hemisphere via the TOA, but their effects are mostly counteracted by the albedo and temperature feedback (Fig. 3b, Extended Data Fig. 1). It is worth noting the hemispheric asymmetries of planetary albedo 34 : further declines in Arctic and Antarctic sea ice by the end of 1000 years cause large increases in TOA radiation in both polar regions via the albedo feedback (Extended Data Fig. 1i). However, both large radiation increases cancel out so that the albedo feedback contributes far less to the inter-hemispheric asymmetry of TOA radiation than other feedback.
In addition, we notice that changes in the net energy input to the atmosphere at the equator have an increasing contribution to ITCZ movements after the first two decades (Fig. 1d). To estimate this contribution, we compare the latitudes of the energy flux equator in two cases, one with the net energy input from the CO2 quadrupling simulation and the other with the net energy input fixed at its preindustrial level (Method). We find that changes in the net energy input to the atmosphere at the equator can reduce the southward ITCZ shift by 39.5% over years 100-1000.
We further investigate a subset of LongRunMIP of three models with simulation times of at least 4000 years (referred to as LongRunMIP_sub, Method), which is sufficient for the Earth9s climate system to achieve a new (quasi-)equilibrium state following CO2. perturbation. We find that the global ITCZ for LongRunMIP_sub multi-model mean changes similarly to that of LongRunMIP over the first millennium, followed by a persistent southward shift (Fig. 1e and f). Between years 100 and 4000, the energy flux equator and tropical precipitation centroid exhibit significant southward migration trends of 0.25 degrees per millennium (p<0.01) and 0.03 degrees per millennium (p<0.01), respectively (Method). Both metrics suggest that the (quasi-)equilibrium ITCZ response by 4000 years appears as a southward shift from preindustrial levels, which differs from the transient ITCZ response during the first few decades (Fig. 1e and f, Fig. 2e).
The southward displacement of the (quasi-)equilibrium ITCZ response can also be explained using atmospheric energy-flux theory. By the end of 4000 years, the quadrupled CO2 has caused an anomalous northward atmospheric cross-equatorial energy transport, primarily due to the anomalous inter-hemispheric asymmetry of ocean heat uptake (Extended Data Fig. 4). Compared (quasi-)equilibrium ITCZ response over 4000 years, which shows a southward shift from preindustrial levels, in contrast to the northward shift of the transient ITCZ response during the first two decades. It merits attention that the time-dependent ITCZ response discussed here is different from that to volcanic eruptions 41 from the perspectives of both forcing scenario and time scale.
Our findings shed light on the role of AMOC change in ITCZ shifts as a result of global warming. Previous freshwater hosing experiments 4 show that an AMOC slowdown caused by ice sheet melt into the North Atlantic can give rise to a southward displacement of the ITCZ owing to abated northward oceanic heat transport across the Atlantic. Our CO2 forcing scenario, however, differs from this freshwater forcing scenario. The LongRunMIP ensemble mean simulates a CO2-induced AMOC deceleration in the first century but a subsequent AMOC recovery 20,42 (Extended Data Fig. 7a). This strengthened AMOC over the next 900 years coincides with a southward ITCZ migration, which differs from the results of freshwater hosing experiments. The underlying reason is that, rather than the AMOC, delayed Southern Ocean warming and reduced heat uptake dominate the southward ITCZ shift during this time. An additional analysis reveals that AMOC recovery leads to a trend of ocean heat transport convergence and hence a decline trend of ocean heat uptake 43 in the North Atlantic (Extended Data Fig. 7b). This decline in Atlantic Ocean heat uptake contributes to a decrease in surface energy flux in the Northern Hemisphere within 30 o N-65 o N (Extend Data Fig. 7b). Nonetheless, the reduced Southern Ocean heat uptake (30 o S-65 o S) is even faster and stronger than its counterpart (30 o N-65 o N), and thus essentially controls the change in interhemispheric asymmetry of surface energy flux over years 100-1000 (Extend Data Fig. 7c). Such dominant role of Southern Ocean heat uptake is robust across models, regardless of AMOC recovery speed uncertainty among models 42 . Our findings underline the significance of the Southern Ocean heat uptake 44 in the long-term ITCZ evolution under climate change.
Our study shows a non-monotonic ITCZ migration under the simple atmospheric CO2 forcing.
The ITCZ migration may become more complex in future scenarios of representative concentration pathways and shared socioeconomic pathways that include other forcings such as anthropogenic aerosols and stratospheric ozone 6 , or in scenarios of CO2 ramp-up and rampdown 45 , or with nonlinearities in ocean warming patterns on century to millennium time scales 46 .
For example, the AMOC exhibits a clear hysteresis under CO2 ramp-up and ramp-down forcings, which contributes to an ITCZ hysteresis 45 . This is because, following the CO2 forcing turnabout, the AMOC weakens further and reaches its minimum value, causing a Northern Hemisphere cooling and a negative atmospheric net energy input, which promotes the Northern Hemisphere poleward atmospheric energy transport and amplifies the inter-hemispheric energy transport contrast. These changes in the AMOC and ITCZ systems, however, include either their direct responses to the varying CO2 forcing or the adjustments due to feedback in both systems. The constant CO2 forcing in our study, on the other hand, excludes the influence of changes in CO2 forcing and thus allows for a comprehensive analysis of the adjustments within the climate system on different timescales. times) in the tropical precipitation centroid (& , Method, multi-model mean, green; intermodel spread, one standard derivation among models, light green) and atmospheric crossequatorial energy transport (& , multi-model mean, purple; inter-model spread, light purple; 1 PW = 10 15 Watt) in the CO2 quadrupling simulations by (a) CMIP5/6, (c) LongRunMIP and (e) LongRunMIP_sub models. (b,d,f) Same as (a,c,e) but for changes in the energy flux times) in the zonal mean (weighted) CO2-induced energy flux changes (multi-model mean, line; inter-model spread, one standard derivation among models, shading) at the TOA (red) and surface (blue), and their difference (TOA minus surface, blue) in the CO2 quadrupling simulations by (a) CMIP5/6 and (b) LongRunMIP over years 1-20. (c) Same as (b) but for LongRunMIP models for the difference between years 981-1000 and years 1-20. preindustrial times, in units of W/m 2 ) in the net (a) TOA radiation and (b) surface energy flux in the CO2 quadrupling simulation for the multi-model mean of LongRunMIP models over years 1-20. (c,d) Same as (a,b) but for years 981-1000. (e,f) The differences between the two periods for the net (e) TOA radiation and (f) surface energy flux (years 981-1000 minus years 1-20). (g,h) Same as (e,f) but for surface (shortwave plus longwave) radiation energy flux and surface turbulent (sensible plus latent) heat flux. Stippling refers to the region where at least 33 of 43 CMIP5/6 models, or 6 of 8 LongRunMIP models agree with the sign changes. units of K) in the CO2 quadrupling simulation for the multi-model mean of LongRunMIP models over (a) years 1-20 and (b) years 981-1000, respectively. (c) Same as (a) but for the difference between years 981-1000 and years 1-20. Stippling refers to the region where at least 6 of 7
LongRunMIP models agree with the sign changes (ECHAM5-MPIOM is not included since SST data are not available for its CO2 quadrupling simulation).
We use preindustrial and CO2 quadrupling simulations with 43 CMIP5/6 models (Supplementary Table 1). These models are chosen primarily due to the availability of data for the kernel calculation. The length of the CO2 quadrupling simulation varies between models but is at least 150 years; we use 150-year simulation outputs for all models. To ensure that each model receives an equal amount of weight in the inter-model analysis, only one ensemble member is chosen from each model.
Furthermore, we use preindustrial and CO2 quadrupling simulations from an eight-model ensemble (Extended Data Table 1), which includes seven LongRunMIP climate models and one CMIP6 model (ACCESS-ESM1-5) not included in the aforementioned CMIP5/6 models. The length of the CO2 quadrupling simulation varies between models but is at least 1000 years; we use 1000-year simulation outputs for all models. There is also a LongRunMIP subset of three models (CESM104, GISS-E2-R, and MPI-ESM1-1) with CO2 quadrupling simulations lasting more than 4000 years (referred to as LongRunMIP_sub, Extended Data Table 1). For these three models, we use 4459-year simulation outputs.
We realize that the double ITCZ bias remains an issue in several generations of climate models [47][48][49][50][51] . This ITCZ bias has been suggested to be linked to Southern Ocean cloud bias 48,50 , however, the teleconnection between the Southern Ocean and tropical precipitation biases is muted by adjustments in energy transports in the coupled climate system 15,16 . Furthermore, a direct relationship between the mean-state double ITCZ bias and ITCZ changes is not statistically significant 52 .
The overturning Hadley circulation transports moist static energy in the direction of its upper branches, that is away from the ITCZ. Since the eddy contribution to the tropical atmospheric energy transport is negligible in comparison to the overturning Hadley circulation contribution, the zonal-mean ITCZ should lie near the <energy flux equator= where the atmospheric meridional energy transport alters sign 7,[31][32][33] . According to the atmospheric energy balance, the energy flux equator ( ) can be expressed as
where denotes the radius of Earth. Eq. ( 1) states that, to the first order, the energy flux equator is determined by the atmospheric cross-equatorial energy transport ( ) and the net energy input to the atmosphere at the equator ( ). When the temporal change of atmospheric energy storage is neglectable on decadal or longer timescales, can be calculated as
The change in TOA radiation can then be divided into components caused by temperature, water vapor, albedo, cloud feedback, and a residual term. Planck and lapse rate feedback are included in the temperature feedback, and both shortwave and longwave cloud feedback are included in the cloud feedback. Due to data availability for the kernel calculation, the above decomposition of TOA radiation change is only applied to ACCESS-ESM1-5 for the LongRunMIP.
We adopt two metrics to estimate the location of the zonal-mean ITCZ. The first one is the latitudinal centroid of tropical precipitation:
where = 20 and = 20 are the latitudinal integration bounds, and is zonal mean precipitation 32 . The second one is the energy flux equator ( ) 31,32 . We compute the change in each metric with respect to its preindustrial control in the CO2 quadrupling simulation. For instance, the changes in the energy flux equator and tropical precipitation centroid are represented by & and & , respectively. We further quantify the contributions of and changes to ITCZ shifts by defining
where the net energy input to the atmosphere at the equator is fixed at its preindustrial level ( ). The difference between & and & reverals the effect of changes in the net energy input on ITCZ shifts.