Evaluating the large-scale hydrological cycle response within the Pliocene Model Intercomparison Project Phase 2 (PlioMIP2) ensemble

The mid-Pliocene (∼ 3 Ma) is one of the most recent warm periods with high CO2 concentrations in the atmosphere and resulting high temperatures, and it is often cited as an analog for near-term future climate change. Here, we apply a moisture budget analysis to investigate the response of the large-scale hydrological cycle at low latitudes within a 13-model ensemble from the Pliocene Model Intercomparison Project Phase 2 (PlioMIP2). The results show that increased atmospheric moisture content within the midPliocene ensemble (due to the thermodynamic effect) results in wetter conditions over the deep tropics, i.e., the Pacific intertropical convergence zone (ITCZ) and the Maritime Continent, and drier conditions over the subtropics. Note that the dynamic effect plays a more important role than the thermodynamic effect in regional precipitation minus evaporation (PmE) changes (i.e., northward ITCZ shift and wetter northern Indian Ocean). The thermodynamic effect is offset to some extent by a dynamic effect involving a northward shift of the Hadley circulation that dries the deep tropics and moistens the subtropics in the Northern Hemisphere (i.e., the subtropical Pacific). From the perspective of Earth’s energy budget, the enhanced southward cross-equatorial atmospheric transport (0.22 PW), induced by the hemispheric asymmetries of the atmospheric energy, favors an approximately 1 northward shift of the ITCZ. The shift of the ITCZ reorganizes atmospheric circulation, favoring a northward shift of the Hadley circulation. In addition, the Walker circulation consistently shifts westward within PlioMIP2 models, leading to wetter conditions over the northern Indian Ocean. The PlioMIP2 ensemble highlights that an imbalance of interhemispheric atmospheric energy during the mid-Pliocene could have led to changes in the dynamic effect, offsetting the thermodynamic effect and, hence, altering mid-Pliocene hydroclimate.


Introduction
Global warming can induce regional and global anomalies in the Earth's hydrological cycle, thereby regulating the balance of global water resources (Eltahir and Bras, 1996). Many studies have indicated that pronounced climate change can occur as anthropogenic CO 2 rises, including an increase in surface temperature Long et al., 2014), Arctic amplification (Stuecker et al., 2018;Smith et al., 2019), and impacts on animal and plant populations (Root et al., 2003). Under current global warming, both observations and model simulations suggest a tendency for the "wet regions getting wetter and dry regions getting drier" phenomenon (Held and Soden, 2006;Wentz et al., 2007;Chou et al., 2009;Wang et al., 2012;Li et al., 2013). That is, precipitation minus evaporation (PmE) increases (decreases) in regions of climatological convergence (divergence). Note that this phenomenon is primarily focused on the ocean. A study by Greve et al. (2014) reported that only 10.8 % of the global land area shows the dry gets drier and wet gets wetter pattern. These changes in the large-scale hydrological cycle could induce severe climatic disasters worldwide, leading to considerable impacts on economies, ecosystems, and agriculture (Asokan and Destouni, 2014;Bengtsson, 2014). Therefore, understanding the potential processes responsible for large-scale hydrological cycle changes in a warmer climate is of great importance.
Previous studies have suggested that the thermodynamic effect caused by increased atmospheric moisture content in a warmer climate is one of the primary contributors to a tendency toward wet gets wetter and dry gets drier conditions (Chou et al., 2009;Seager et al., 2010). This mechanism directly follows the nonlinearity of the Clausius-Clapeyron relationship, which acts to increase atmospheric moisture content over regions with the warmest surface temperatures (Allen and Ingram, 2002;Stephens and Ellis, 2008). On the other hand, large-scale atmospheric circulation can change substantially due to nonuniform temperature changes under global warming and, hence, induce changes in the hydrological cycle via the so-called dynamic effect (Han et al., 2019a). The dynamic effect is relatively more complicated than the thermodynamic effect among climate models. Seager et al. (2010) demonstrated that the dynamic component is modulated by the weakening of the Hadley circulation and Walker circulation. An increased CO 2 concentration could directly increase atmospheric static stability over tropical oceans, favoring a slowdown of these atmospheric overturning circulations (Vallis et al., 2015). Other studies have indicated that the local Hadley circulation shifts poleward due to the decreased meridional temperature gradient in response to increased CO 2 concentrations (Sharmila and Walsh, 2018;Y. Hu et al., 2018). These circulation anomalies widen the subtropical dry zones (Previdi and Liepert, 2007;Sun et al., 2013a). In addition, Long et al. (2016) highlighted that model uncertainty in tropical rainfall comes from the discrepancies in the atmospheric circulation anomalies among models. Thus, the spread of circulation changes in response to global warming across climate models leads to a diversity in responses in the hydrological cycle.
Proxy data indicate that the mid-Pliocene (∼ 3 Ma) was one of the most recent warm periods with CO 2 levels similar to the current anthropogenically elevated value of 400 ppm and can be considered an analog for future climate change (Dowsett et al., 2012;Burke et al., 2018;Tierney et al., 2019). Pliocene Model Intercomparison Project Phase 1 (PlioMIP1) simulations have been used to investigate how the climate system responded to mid-Pliocene boundary conditions, including elevated atmospheric CO 2 concentrations. These past warm climate simulations exhibit many similarities with future climate projections. For example, one robust characteristic is increased temperature from 1.8 to 3.6 • C during the Pliocene compared with the preindustrial period (PI) , with Arctic amplification in response to a significant decline in sea-ice extent (Howell et al., 2016;Zheng et al., 2019). These features could have reduced the meridional surface temperature gradient, inducing weaker tropical circulation (i.e., local Hadley circulation) during the Pliocene (Sun et al., 2013b;Li et al., 2015;Corvec and Fletcher, 2017). Additionally, some studies have suggested a weakened zonal sea surface temperature (SST) gradient in the Pacific during the Pliocene (Wara et al., 2005;Scroxton et al., 2011), which would have favored a weaker Walker circulation. These features could have induced large-scale changes in Pliocene hydroclimate. Using a climate simulation that captures the warming patterns seen in early-Pliocene sea surface temperature proxies, Burls and Fedorov (2017) suggested that the dynamic process might play a key role in driving wetter subtropics due to this weaker tropical circulation during the early-Pliocene warm climate compared with the future climate.
Although PlioMIP1 can reproduce similar patterns of the change in surface temperature to the reconstructed SST, models cannot capture the magnitude of warming at higher latitudes. For example, Dowsett et al. (2013) indicated that the ensemble of PlioMIP1 models underestimates the warming in the North Atlantic compared with the reconstructed SST. This might be induced by the uncertainties in PlioMIP1, including the uncertainty in atmospheric CO 2 concentrations Howell et al., 2016) and paleogeography and bathymetry Feng et al., 2017). PlioMIP2 models show the closed Canadian Archipelago and Bering Strait and a reduced Greenland ice sheet relative to PlioMIP1. For one of the PlioMIP2 models, it has been shown that updating the paleogeography to PRISM4 is the major contributor to climate differences from PlioMIP1 to PlioMIP2 . In addition, PlioMIP2 focuses on a specific time slice during the mid-Pliocene at approximately 3.025 Ma, which could reduce the uncertainties in reconstructions (McClymont et al., 2020). Researchers have been investigating the mid-Pliocene climate by using PlioMIP2, including Arctic warming (De Nooijer et al., 2020), Atlantic meridional overturning circulation (Z. , climate sensitivity , global monsoons (Q. , and subtropical rainfall changes (Pontes et al., 2020). However, it is difficult to distinguish the relative impact of the Hadley circulation and Walker circulation on Pliocene hydrological cycling at low latitudes. Fortunately, the three-pattern decomposition of global atmospheric circulation (3P-DGAC; Hu et al., 2017;S. Hu et al., 2018a, b) method can help us to decompose atmospheric circulation into zonal (i.e., local Walker circulation) and meridional (i.e., local Hadley circulation) circulation at low latitudes. We apply this method to develop moisture budget analyses, which might provide some insight into the mechanisms of hydrological cycling during the mid-Pliocene. This paper set is in the framework of updated PlioMIP2 models to quantitatively distinguish the relative contribution from zonal and meridional circulation anomalies to hydro-logical cycle changes. In the following section, we first introduce the PlioMIP2 models and moisture budget decomposition. We then evaluate the simulated large-scale hydroclimate cycle response within the PlioMIP2 ensemble in Sect. 3. Section 4 provides each moisture budget component's relative contribution to investigate the potential mechanisms driving the simulated changes in the mid-Pliocene hydrological cycle. The corresponding mechanisms are discussed in Sect. 5. The last section contains the conclusion and discussion.

Climate model simulations
In this study, we use the simulations from 13 models participating in PlioMIP2 (Table 1). All models include a preindustrial (PI) simulation and a Pliocene climate simulation. In PlioMIP2 models, the boundary conditions have been updated using the new version of the U.S. Geological Survey PRISM4 dataset , including soils, lakes, land-ice cover, vegetation, topography, and bathymetry. The CO 2 levels for the mid-Pliocene and PI simulations are set at 400 and 280 ppmv, respectively. To calculate the ensemble mean, we interpolate all data onto a common grid with a 1 • × 1 • resolution using bilinear interpolation.

Development of moisture budget decomposition
To examine the changes in precipitation (P ) minus evaporation (E) in the PlioMIP2 mid-Pliocene experiments relative to their respective PI simulation, we decompose the moisture budget equation based on Seager et al. (2010), i.e., Here, g is gravity, ρ w is the density of water, V is the horizontal wind, q is the specific humidity, δ(·) is the annual mean difference in variables between the warmer climate state (mid-Pliocene) and the PI simulation, and subscript 0 represents the variables in the PI simulation. In the warmer climate, the change in P minus E (PmE, the left-hand side of Eq. 1) is balanced by the thermodynamic (δTH, induced by increased specific humidity) and dynamic (δMCD, induced by circulation anomalies) contributions and residual term (R, which is mainly involved in the contributions from high-frequency variability of transient eddies, nonlinear effects, and surface boundary terms). As we are interested in understanding the relative contribution from zonal circulation (i.e., local Walker circulation) changes and meridional circulation (i.e., local Hadley circulation) anomalies to the changes in PmE in a warmer climate, we further apply the three-pattern decomposition of global atmospheric circulation (3P-DGAC; Hu et al., 2017;  S. Hu et al., 2018a, b) method in this study. The horizontal, meridional, and zonal circulations that can be viewed as the global generalization of the Rossby wave in the middle-high latitudes and the Hadley and Walker circulations in the low latitudes are defined to decompose the global atmospheric circulation into a superposition of the horizontal, meridional, and zonal circulations by using the 3P-DGAC method. Based on the essential features of the Rossby, Hadley and Walker circulations, Hu et al. (2017) defined the 3D horizontal circulation V R , meridional circulation V M , and zonal circulation V Z in the spherical σ -coordinate system as follows: Here, the following continuity equations are satisfied: Equation (3) is the sufficient condition that the components of V R , V M , and V Z can be represented by the stream functions R(λθ σ ), H (λθ σ ), and W (λθ σ ), respectively, as fol-lows: Because three-pattern circulations (horizontal, meridional, and zonal circulations) exist in both the low latitudes and middle-high latitudes, the global atmospheric circulation can be expressed as the superposition of the horizontal, meridional, and zonal circulations -that is, with the following components: Equation (5) or (6) is called the three-pattern decomposition model. In contrast to the traditional 2D decomposition of the atmospheric motion into vortex and divergent parts, the continuity Eq. (5) cannot guarantee the uniqueness of the stream functions R(λθ σ ), H (λθ σ ), and W (λθ σ ) because the threepattern circulations V R , V M , and V Z have three spatial dimensions, respectively (Hu et al., 2017;S. Hu et al., 2018a, b). The following restriction condition is needed to pick up the correct decomposition (Theorems 1 and 2 in Y. Hu et al., 2018): Equation (7) guarantees both the uniqueness of the stream functions R, H , and W and the physical rationality of the 3P-DGAC method.
Using the 3P-DGAC method, we can rephrase the moisture budget in Eq. (1) to involve the contributions from zonal and meridional circulation. Here, we neglect the relatively smaller terms at low latitudes, including transient eddies, nonlinear effects, and surface boundary terms. Thus, we mainly explore the contributions from δTH and δMCD to changes in PmE in this study. Then, the δTH and δMCD can be rewritten as follows: where subscripts D and A represent the terms that are related to divergence and moisture advection, respectively. In addition, the subscripts R, Z, and M indicate the terms that are related to the horizontal, zonal, and meridional circulations, respectively. Note that V R represents the horizontal vortex winds, which are not divergent, which indicates that the terms that are related to the divergence/convergence of V R (i.e., δTH D_R and δMCD D_R ) are zero. These terms can be clearly seen in Figs. 3h and 4h. In addition, we ignore these two terms in this study.

Changes in precipitation minus evaporation (PmE) in the PlioMIP2 models
The last 100 years of individual PlioMIP2 simulations are used to calculate the multi-model mean (MMM) PmE in Fig. 1 and individual PlioMIP2 models in Fig. 2. Figure 1a shows that most subtropical regions experience reduced PmE in the mid-Pliocene simulations with respect to the PI simulations, including the subtropical Pacific and subtropical Atlantic in both hemispheres and the subtropical Indian Ocean in the Southern Hemisphere (SH). There is also drying over the South Pacific convergence zone (SPCZ), except in the GISS-E2-1-G, COSMOS, and HadGEM3 models ( Fig. 2), consistent with other studies evaluating the hydrological cycle response within the PlioMIP2 simulations (Pontes et al., 2020). Note that there is a moistening signal in the southern part of the SPCZ in the tropical southern Pacific. In contrast, the increased MMM PmE is located in the deep tropics (i.e., Pacific intertropical convergence zone, ITCZ, and northern Indian Ocean) as well as at middle-high latitudes (Fig. 1a). However, some models (i.e., the CESM2, GISS-E2-1-G, COSMOS, and HadGEM3 models) show a drier Maritime Continent (Fig. 2), which might be related to the changes in Walker circulation (we will discuss this latter in Sect. 5.3). In addition, the North African and Southeast Asian monsoon regions also show significant moistening signals, which are consistent with faunal remains and palynological transfer functions (Sanyal et al., 2004;Trauth et al., 2007;Xie et al., 2012) as well as with other modeling studies Li et al., 2020;. Zhang et al. (2016) indicate that the combined influence of SST and CO 2 level, as well as the vegetation changes, play a very important role in changing the atmospheric circulation over North Africa during the mid-Pliocene, owing to the increased net atmospheric energy there. Additionally, the expansion of vegetation into the Sahara region tends to decrease the surface albedo, which can enhance the Saharan heat low and, hence, impact rainfall over West Africa, reflecting the vegetation-albedo feedback (Charney, 1975). Recent studies indicate that the enhanced vegetation in the PlioMIP2 ensemble is likely to have contributed to increased mid-Pliocene West African summer rainfall Berntell et al., 2021). This change over Southeast Asia is robust among PlioMIP2 models, and only the COSMOS model shows a drier change over East Asia (Fig. 2e). Furthermore, the MMM PmE changes over Southeast Asia are mainly focused on the summertime (not shown), suggesting a consequence of strengthened East Asian summer monsoon circulation (Salzmann et al., 2008;Wan et al., 2010;Yan et al., 2012;Zhang et al., 2013;Li et al., 2018;Lu et al., 2021). Note that the mid-to high-latitude North Atlantic becomes drier (Fig. 1a). The response of the hydrological cycle during the mid-Pliocene generally shows a wet regions getting wetter and dry regions getting drier pattern, especially over the ocean. These features are apparent in the zonal average of the PmE change ( Fig. 1b), except in the GISS-E2-1-G model (Fig. 2f). The tropical regions become wetter, and subtropical regions become drier, which are similar to the results from future high-CO 2 scenario experiments (Chou et al., 2009). Earlier studies have indicated that these features of changes in PmE at low latitudes are linked to the increased specific humidity (i.e., changes in the thermodynamic effect). However, there are some opposite phenomenon as well, when looking at the regional changes in PmE (i.e., northern Indian Ocean, North Africa, and the SPCZ). These may suggest that another factor, such as atmospheric circulation anomalies (i.e., changes in dynamic effect), may play an important role in changing the regional PmE pattern at low latitudes.

Previous model-data comparisons of hydrological changes in the PlioMIP2 ensemble
Multi-proxy studies are qualitatively consistent with the results of the PlioMIP2 ensemble . For several studies, proxy reconstructions suggest an expansion of woodland and a higher density of land cover over northern Africa, indicating moistening signals there (Salzmann et al., 2008;Bonnefille, 2010). The sedimentological indicators and pollen data also suggest a more humid climate over the Levant and Arabian peninsulas during the mid-Pliocene (Munoz et al., 2002;Heermance et al., 2013). In addition, the faunal remains or palynological transfer functions show wetter conditions in East and South Asia during the mid-Pliocene (Sanyal et al., 2004;Igarashi and Yoshida, 1988;Kou et al., 2006). However, uncertainties related to the hydroclimate of the mid-Pliocene still remain. For instance, pollen evidence suggests little hydroclimate change during the Pliocene in Qaidam Basin and southwest China's Yuanmou region (Wang et al., 1999;Chang et al., 2010;Heermance et al., 2013). Some proxies even show a drier climate over the Loess Plateau region (Ji et al., 2017;Sun et al., 2010). Note that the relatively low availability of Pliocene hydroclimate proxies makes it difficult to perform a model-proxy comparison. Furthermore, PlioMIP2 modeling experiments are designed to simulate the Marine Isotope Stage KM5c (3.205 Ma) during the mid-Pliocene, and this particular orbital interval likely does not represent the full Pliocene hydroclimate variability, adding uncertainties to model-proxy comparison .

Thermodynamic and dynamic contributions to changes in PmE
Moisture budget analyses are conducted to shed light on the mechanisms driving the changes in PmE during the mid-Pliocene. Based on this decomposition, the changes in PmE are mainly influenced by the changes in humidity with unaltered atmospheric circulation (called the thermodynamic term, δTH) and changes in atmospheric circulation with no change in humidity (called the dynamic term, δMCD) at low latitudes. The thermodynamic term (δTH) and its decomposition are plotted in Fig. 3. It is clear that δTH captures the main features of hydrological cycle change (Figs. 3a vs. 1a) -that is, the positive and negative contributions over the already convergent (i.e., the ITCZ and SPCZ) and divergent (subsidence of local Hadley circulation) regions, respectively. In general, the thermodynamic term increases PmE by ∼ 58.6 % over the tropics, and decreases PmE by ∼ 84.6 % over subtropics (not shown), respectively. This term does not alter the spatial distribution of climatological PmE (contours in Fig. 1a) but amplifies the intensity of the existing pattern of PmE, reflecting the wet getting wetter and dry getting drier mechanisms (Held and Soden, 2006). These results are consistent with future global warming scenarios (Chou et al., 2009;Wang et al., 2012;Li et al., 2013). From the perspective of global atmospheric circulation, previous studies have indicated that global atmospheric circulation can be decomposed into a superposition of horizontal, meridional, and zonal circulations (Hu et al., 2017;S. Hu et al., 2018a, b). δTH is further decomposed using the 3P-DGAC method (Fig. 3c-k). The estimated δTH in Fig. 3b, calculated as the sum of the right-hand side in Eq. (8) of the 3P-DGAC decomposition method, shows a similar distribution to the δTH field shown in Fig. 3a with a pattern correlation coefficient (PCC) of 0.80. This result indicates that the decomposition is representative. At low latitudes, the δTH mainly comes from terms that are related to climate mean meridional and zonal circulation (Fig. 3c, d), whereas at middle-high latitudes, the δTH mainly comes from horizontal circulation (Fig. 3e). It is clear that the thermodynamic changes associated with meridional circulation can explain the large portion of δTH (PCC of 0.9) at low latitudes, which is caused by increased specific humidity within the di- vergence of climate mean meridional circulation (δTH D_M ; Fig. 3f). The zonal circulation can also explain δTH to some extent, with a positive contribution mainly over the Maritime Continent extending eastward to the equatorial central Pacific and eastern coast of North/South America, and a negative contribution over the eastern Pacific extending from the western Indian Ocean to the Greater Horn of Africa and the eastern tropical Atlantic (Fig. 3d). These changes associated with zonal circulation are linked to the increased specific humidity with divergence of the mean zonal circulation (δTH D_Z ; Fig. 3g). At middle-high latitudes, the δTH induced by climate mean horizontal circulation is caused by changes in moisture advection (Fig. 3k), e.g., the western coast of North America, a region extending from the southern tip of South America to the central tropical Pacific Ocean, and southern tip of South Africa.
It is evident that the δTH component does not describe the full contribution to the changes in PmE, especially over the North African and Southeast Asian monsoon regions, the SPCZ, and the northern Indian Ocean, where we must con-sider the dynamic effect. The dynamic effect (δMCD), reflecting the impact of circulation changes, partially offsets the δTH at low latitudes (Fig. 4a). In particular, δMCD reduces PmE in the deep tropics, i.e., the ITCZ, the SPCZ, and the Maritime Continent. In contrast, δMCD can moisten subtropical regions, especially over the subtropical eastern Pacific, southern Indian Ocean, and Atlantic Ocean of both hemispheres. Compared with δTH, the dominating contribution from δMCD to changes in PmE lies adjacent to the northern Indian Ocean, the SPCZ, and the North African and Southeast Asian monsoon regions (Fig. 4a).
The estimated δMCD in Fig. 4b, calculated as the sum of the right-hand side terms in Eq. (9) of the 3P-DGAC decomposition method, is consistent with the δMCD in Fig. 4a with a PCC of 0.93. This result indicates that the decomposition is representative. The anomalous divergence of the meridional circulation component (δMCD D_M ) appears to dry the deep tropics but moisten the Northern hemispheric part of the deep tropics, which is associated with the northward shift of the ITCZ (Fig. 7c). In particular, the northward shift of the ITCZ is clear from 150 • E to the east in the Pacific (Fig. 4f). In addition, the component contributes a large portion to enhance PmE over the North African and Southeast Asian monsoon regions. However, the tropical southern Pacific is even more complicated. The δMCD D_M term contributes to reduced PmE over the SPCZ region but to increased PmE over the southern part of the SPCZ in the tropical southern Pacific. A previous study has suggested that these changes in PmE followed the southward shift of the SPCZ, which was mainly modulated by the intensified and westward shift of the South Pacific subtropical high for the mid-Pliocene compared with the PI simulation (Pontes et al., 2020). For the adjacent northern Indian Ocean, the convergence of zonal circulation anomalies (δMCD D_Z ) is the first-order contribution to strengthen the dynamic effect (by ∼ 45 %) and, hence, enhances the PmE (Fig. 4g).
In summary, the dynamic and thermodynamic terms can explain the largest changes in PmE at low latitudes (Figs. 5 vs. 1). The thermodynamic term induced by the divergence of the mean meridional circulation is the dominant process driving changes in PmE at low latitudes (Fig. 3f). However, the dynamic term partially offsets δTH, especially over the ITCZ, the SPCZ, and the Maritime Continent, via changes in the divergence of meridional circulation. Even the dynamic term overwhelmingly contributes to the increased PmE over the North African and Southeast Asian monsoon regions and the northern Indian Ocean (Figs. 4 vs. 5). Note that the former two are mainly caused by meridional circulation anomalies, but the latter is dominated by zonal circulation anomalies.
We further decompose the meridional moisture transport into terms that reflect the changes in specific humidity (meridional moisture transport induced by the thermodynamic effect; MMTT) and circulation (meridional moisture transport induced by the dynamic effect; MMTD) in Fig. 6a. As expected, all models show that the MMTT is responsible for the wetter tropics and drier subtropics in the mid-Pliocene simulation, indicating a dry gets drier and wet gets wetter mechanism. These features are robust among models and are associated with the increased specific humid- ity combined with the mean meridional circulation from the PI control (Fig. 6b), as mentioned above. This is because the zonal-mean wind depicts southerly (northerly) wind between the Equator and subtropical SH (Northern Hemisphere; NH) for the climate mean meridional circulation in the PI simulations. When the climatological wind is combined with increased specific humidity in the low-level troposphere (Fig. 7b), more moisture is transported from the subtropics to the tropics, resulting in a drier subtropics and wetter tropics. In contrast, the MMTD shows a large spread across PlioMIP2 models. On average, the anomalous MMTD appears to weaken thermodynamic contributions in the subtropical NH but strengthen it in the subtropical SH via meridional circulation anomalies (Fig. 6c). This indicates that the changes in MMTD favor the transport of more (less) moisture from the tropics to the NH (SH) subtropics, which is caused by the northward shift of the meridional circulation (as detailed further in Sect. 5.2). The equatorward moisture transport anomalies of the dynamic component in the SH are due to anomalous southerly winds in the subtropical southern Pacific (Fig. 9). This feature acts to dry the SPCZ and moisten the southern SPCZ and the equatorial central-eastern Pacific (Fig. 4f).

Mechanisms for the changes in moisture budget components
Thus far, we have shown that the anomalous hydroclimate within the mid-Pliocene simulations involves anomalies of both thermodynamic and dynamic effects at low latitudes. In this section, we further examine the corresponding mechanisms in turn. In the MMM, SSTs are between 1 and 6 • C warmer in the mid-Pliocene simulations than in the PI simulations. Note that the SST warming is amplified in the northwest tropical Indian Ocean, whereas it is reduced off the Indonesian coast, showing a pattern similar to the tropical Indian Ocean dipole (IOD). The sharp SST gradients drive strong southeasterly wind anomalies on the Equator (not shown). Xie et al. (2010) suggest that this easterly wind anomaly may shoal the thermocline in the east, helping lower the SST there via upwelling and indicating this SST anomaly over tropical Indian Ocean may be related to the Bjerknes feedback. The simulated North Atlantic warming might be related to an intensified mid-Pliocene Atlantic meridional overturning circulation . However, Z.  suggest that the increased background ocean vertical mixing parameters could also contribute to the warm SSTs there. In addition, the relative smaller SST warming in Southeast Pacific and Atlantic, which is co-located with the intensified southeast trade winds, suggests the role of windevaporation-SST feedback . These SST warming patterns are consistent with current studies Williams et al., 2021). As expected, the specific humidity is increased in the low-level troposphere in the mid-Pliocene warm period (Fig. 7b) (Murray, 1966;Held and Soden, 2006). On the other hand, the sinking branch of meridional circulation in the control climate is located in subtropical regions, showing divergent circulation ∇ · V > 0 in the low-level troposphere; the contrary applies for the regions of deep tropics, i.e., the ITCZ and SPCZ. These two factors contribute the δTH D_M term (i.e., − 1 ρg p s 0 δq∇ · V M0 dp) to the thermodynamic effect (Fig. 3f) and, hence, changes in PmE.

Changes in specific humidity
Although the δTH D_M term is the first-order control on the thermodynamic effect in most regions, the δTH D_Z term contributes to the thermodynamic effect to some extent, especially over the adjacent northern Indian Ocean. The climate mean zonal circulation characterizes ascending motion in the tropical western Pacific, tropical African, and tropical southern American regions, favoring convergent circulation (i.e., −∇ · V Z 0 > 0) there (Fig. 6d;Hastenrath, 1991;Peixoto and Oort, 1992). With increased specific humidity (δq > 0) under a warmer climate, the δTH D_Z term (i.e., − 1 ρg p s 0 δq∇ · V Z0 dp) shows a positive contribution and, hence, increases PmE in these regions (Fig. 3g). On the contrary, the δTH D_Z favors a decrease in PmE over the western Indian Ocean, eastern Pacific, and tropical Atlantic (Fig. 3g), where the climate mean zonal circulation is divergent (Fig. 7d).

Response in meridional circulation
In Sect. 4, we have demonstrated that the primary dynamic contribution to changes in PmE is a consequence of anomalous meridional circulation (the δMCD D_M term). Figure 8a shows the annual mean mass stream function (MSF) of meridional circulation for the PI simulation (contours), which is similar to present-day meridional circulation (Cheng et al., 2020). During the mid-Pliocene, the meridional  (Prahl and Wakeham, 1987) and foraminifera calcite Mg/Ca (Delaney et al., 1985). Stippling in panels (b-d) indicates regions where at least 10 of 13 simulations in the model group agree on the sign of the ensemble mean. circulation changes are characterized by enhanced meridional circulation in the SH tropics and weakened meridional circulation in the NH tropics (shading in Fig. 8), which is caused by the northward shift of meridional circulation in the SH, as indicated in our later discussion. To quantify meridional circulation changes, we further calculated the intensity in Fig. 8b. The intensity is defined as the maximum of the absolute average MSF between 200 and 925 hPa in the range from 30 • S to 30 • N (Oort and Yienger, 1996) in Fig. 8a.
Models simulate a consistently weakened meridional circulation intensity in the NH and a slightly strengthened intensity in the SH (Fig. 8b), which is related to the hemispheric asymmetry of the atmospheric energy budget .
As a result, meridional circulation anomalies could induce divergent/convergent circulation anomalies in the low-level troposphere (Fig. 9). The weakened local meridional circulation leads to anomalous southerly winds spanning northeastern South America eastward to the northwestern Pacific region. These meridional circulation anomalies induce the anomalous divergence (convergence) of circulation over the Indo-Pacific warm pool (adjacent to subtropical regions), resulting in a negative (positive) contribution from δMCD D_M (Figs. 9 vs. 4f). In fact, this anomalous meridional circulation is closely related to the strengthened Asian summer monsoon (not shown), consistent with previous studies Prescott et al., 2019). In addition, anomalous northerly winds exist in the western tropical Pacific, but southerly winds are located in the central Pacific (Fig. 9). These circulation anomalies could induce δMCD D_M , which favors moistening of the equatorial central Pacific and southern part of the SPCZ region but dries the SPCZ (Fig. 4f). Previous studies have indicated that these circulation anomalies are caused by the southward shift of the SPCZ, which is mainly modulated by the intensified and westward shift of the South Pacific subtropical high for the mid-Pliocene compared with the PI simulation (Pontes et al., 2020).
One question arises regarding what causes the meridional circulation changes under mid-Pliocene conditions. At low latitudes, it is worth noting that the ITCZ lies at the foot of the ascending branch of the meridional circulation, which is highly linked to the hemispheric asymmetry of the atmospheric energy budget (Frierson et al., 2013). We further quantify the shift of the ITCZ in Fig. 8c and Earth's energy budget in Fig. 8d and e. The definition of the ITCZ location is the latitude of the maximal annual mean precipitation between 20 and 20 • S (Frierson and Hwang, 2012;Donohoe et al., 2013). On average, ensemble models show that the NH atmosphere receives 1.5 W m −2 more net radiation than the SH (Fig. 8d), which could induce an increased crossequatorial southward energy flux of 0.22 PW (Fig. 8e). Thus, this imbalance in the atmospheric energy budget causes a 1.1 • northward shift in the zonal-mean ITCZ latitude. Consequently, this shift of the ITCZ reorganizes atmospheric circulation (Watt-Meyer and Frierson, 2019), leading to the northward movement of the meridional circulation in the SH (Fig. 8a). This meridional circulation shift could result in a weakened (strengthened) meridional circulation in the NH (SH) (Fig. 8b) and, hence, drive δMCD D_M (Fig. 4f) and MMTD (Fig. 6c). Pontes et al. (2021) indicated that the anomalous wind over the southern Pacific is related to the El Niño-Southern Oscillation (ENSO) weakening across models, which could favor a northward shift of the ITCZ. In addition, it should be noted that the northward shift of the ITCZ exists in both the boreal summer (June-July-August) and winter (December-January-February) seasons, accompanied by the northward shift of the meridional circulation in the SH (not shown). Pontes et al. (2021) further suggested that the northward shift of the Pacific ITCZ during austral spring-summer is remarkably related to the ENSO weakening across models, which is associated with the stronger climatological circulation in the SH.

Response in zonal circulation
As mentioned above, δMCD D_Z plays a key role in the changes in PmE over the northern Indian Ocean. As this term is linked to Walker circulation anomalies, we further discuss Walker circulation changes in the mid-Pliocene warm climate.
There is a noticeable diversity in the simulated Pacific Walker circulation (PWC) intensity across the models (Fig. 10a). In addition, previous work has suggested that the PWC intensity is closely tied to the zonal SST and sea level pressure (SLP) gradient during the mid-Piacenzian (Tierney et al., 2019). In this paper, the dSLP and dSST are defined as the difference in SLP and SST across the equatorial Indo-Pacific (160-80 • W, 5 • S-5 • N minus 80-160 • E, 5 • S-5 • N). As expected, the models with an enhanced zonal SST gradient across the equatorial Indo-Pacific tend to produce a weaker zonal SLP gradient and decreased PWC (not shown), with the inter-model correlations of −0.95 and −0.75, re- spectively. Note that the PlioMIP2 models produce a large spread in simulating the changes in dSST (Fig. 10c) and dSLP (Fig. 10d), which is consistent with the results in Fig. 10a. Previous studies have suggested that the eastwest SST gradient was reduced in SST proxies (Tierney et al., 2019). This feature is captured by the CESM2, GISS-E2-1-G, and HadGEM3 models (Fig. 10a, b, c). However, other models (i.e., the CCSM4, CCSM4-Utrecht, EC-Earth3-LR, and NorESM-L models) consistently simulated stronger PWC intensity (Fig. 10a, b, c). That is, the results suggest that the model-simulated changes in the strength of PWC are probably highly model dependent, which might be affected by the different parameterizations (Tierney et al., 2019).
However, the westward shift of the PWC is a robust feature among these models, except the COSMOS and GISS-E2-1-G models (Fig. 10b). To discuss the impact of the PWC shift on atmospheric circulation in the tropics, we further calculate the changes in the zonal mass stream function (ZMS) for the mid-Pliocene with respect to the PI simulation in Fig. 11. As suggested in Fig. 7d, the ZMS in the PI simulation (contours in Fig. 11a) is characterized by ascending motion in the tropical western Pacific and Maritime Continent and descending motion in the western Indian Ocean and eastern Pacific, consistent with previous studies (Kamae et al., 2011;Bayr et al., 2014;Ma and Zhou, 2016;Han et al., 2020). Compared with the PI simulation, the most striking features in the mid-Pliocene simulation are weakened ascending motion over the Maritime Continent and tropical western Pa-cific and strengthened descending motion on the western Indian Ocean, indicating a westward expansion of the PWC (Fig. 11b).
The westward shift of the PWC can also be seen from the potential velocity (Fig. 12). This shows that the center of the anomalous positive values is located in the northern Indian Ocean. In contrast, the center of a negative value exists in the equatorial eastern Pacific and western Atlantic in the low-level troposphere (Fig. 12a). Concurrent, generally opposite anomalies can be seen in the upper-level troposphere (Fig. 12b). Indeed, these features indicate an upward (downward) motion shift from the tropical western Pacific (eastern Pacific) to the west of the Indian Ocean (central Pacific), resulting from the westward expansion of the PWC (Figs. 10b,  11). That is, when divergent/convergent circulations are combined with the climate mean specific humidity (q > 0) in the lower troposphere, they can trigger a negative/positive contribution from the δMCD D_Z term to changes in PmE (Fig. 4g).

Discussion and conclusions
This paper evaluates the changes in the large-scale hydrological cycle during the mid-Pliocene with respect to the PI based on 13 PlioMIP2 simulations. A diagnostic analysis using the moisture budget equation and the Earth's energy budget provides insight into the mechanisms. The main conclusions are summarized in the following.  The PlioMIP2 models show large spatial differences in PmE. The MMM generally depicts a wet regions getting wetter (i.e., the ITCZ, Maritime Continent, and monsoon regions) and dry regions getting drier (i.e., a sinking branch of the Hadley circulation) pattern during the mid-Pliocene warm climate. According to the moisture budget equation, a large part of the changes in PmE at low latitudes are due to the increased specific humidity. However, the thermodynamic component cannot fully explain the changes in PmE. The dynamic effects offset the thermodynamic effects to some extent and even determine a larger contribution to the changes in PmE in the southern tropical Pacific and northern Indian Ocean. We find increased hemispheric asymmetries of the atmospheric energy budget (larger atmospheric energy over NH than SH) during the mid-Pliocene compared with the PI period, which could induce the northward shift of the ITCZ and reorganize atmospheric circulation. These features can result in a weakening meridional circulation in the NH monsoon regions and a strengthening meridional circulation in the SH. In addition, the anomalous meridional circulation can dry the deep tropics but moisten the northern part of the ITCZ. Furthermore, these anomalies dry the SPCZ region and wet its southern part, which is associated with the southward shift of the SPCZ. We also find a robust westward shift in the PWC, which appears to moisten the northern Indian Ocean via the anomalous convergence of zonal circulation.
Our analyses provide a relatively complete understanding of the changes in the large-scale hydrological cycle within the PlioMIP2 ensemble. It is evident that the air could hold more moisture in a warmer climate; thus, the thermodynamic effects amplify the intensity of PmE but do not alter its spatial pattern ( Fig. 2a; Held and Soden, 2006). Note that the hemispheric asymmetries of atmospheric energy could induce regional meridional circulation anomalies and, thus, alter the distribution of PmE anomalies during the mid-Pliocene via the δMCD D_M term at low latitudes. The PlioPMIP2 ensemble simulations suggest that hemispheric asymmetries of atmospheric energy are the key factor altering the spatial pattern of PmE via changes in the local meridional circulation. However, we should note that a noticeable inter-model spread exists in capturing the main features in the past warm climate, particularly for the changes in Walker circulation, such as the large spread in the simulated changes in the intensity of the PWC, dSST, and dSLP in Fig. 10, consistent with previous studies . Further effort to understand the inter-model uncertainty needs to be explored in future work. In addition, previous studies indicate that the storm track (transient eddy component) may play a key role in changes in PmE for middle-high latitudes (Seager et al., 2010;Han et al., 2019a, b). Due to the lack of hourly model data, we mainly discuss the relative contributions from moisture budget components to changes in PmE at low latitudes in this paper. Much more work should be conducted to study the impact of storm tracks on changes in PmE during the mid-Pliocene using hourly data in the future at middlehigh latitudes.
Note that the global temperature during mid-Pliocene is controlled by the combined effects of boundary conditions (e.g., CO 2 level, vegetation, and topography) . Any changes in each boundary condition could induce large-scale hydrological cycling changes. For example, the role of remote biophysical effects in the northern middle-high latitudes is highlighted in driving the variation in monsoon rainfall in low latitudes and the shift of the ITCZ, as needleleaf vegetation expands greatly northward in eastern Eurasia during the mid-Pliocene (Chase et al., 2000;Swann et al., 2014;Mahmood et al., 2014;Zhang and Jiang, 2014;Burls and Fedorov, 2017). Some studies have indicated that uncertainties exist in the boundary conditions of changing South Asian summer monsoon (SASM) hydrological cycling. Sarathchandraprasad et al. (2017) indicated that the tectonically induced reorganization of the Indonesian throughflow can strengthen the SASM during the mid-Pliocene due to the increased cross-equatorial pressure gradient. Recent study by Prescott et al. (2019) highlighted the substantial influence of orbital forcing on the changes in SASM during the mid-Pliocene. The simulations suggest that tectonic uplifts in the South African plateaus can strengthen the SASM as well (Zhang and Liu, 2010). Based on these studies, the boundary conditions applied by the PlioMIP2 models can impact the low-latitude hydrological cycle during the mid-Pliocene. However, the relative impact of boundary conditions on hydrological cycling still remains uncertain. Moreover, not all models carry out the sensitivity experiments designed in PlioMIP2, increasing the difficulty related to exploring their relative contributions to PmE changes. These questions need to be considered further in the future.
Author contributions. QZ and ZH designed the work, ZH wrote the paper under the supervision of QZ. ZH carried out the analyses and programming with help from JC and QW. All of the other coauthors provided the PlioMIP2 model data and commented on the paper.

Competing interests.
The contact author has declared that neither they nor their co-authors have any competing interests.
Disclaimer. Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Special issue statement. This article is part of the special issue "PlioMIP Phase 2: experimental design, implementation and scientific results". It is not associated with a conference.
Acknowledgements. Zixuan Han acknowledges financial support from the National Natural Science Foundation of China (grant no. 42130610), the Fundamental Research Funds for the Central Universities (grant no. B210201009), and the National Key R&D Program of China (grant no. 2017YFC1502303). Jianbo Cheng acknowledges financial support from the National Natural Science Foundation of China (grant no. 42005012) and the Natural Science Foundation of Jiangsu Province (grant no. BK20201058). Qin Wen acknowledges financial support from the National Natural Science Foundation of China (grant no. 42106016) and a project funded by the China Postdoctoral Science Foundation (grant no. 2021M691623). The EC-Earth3 model simulations and the data analysis were performed using the ECMWF computing and archive facilities and the Swedish National Infrastructure for Computing (SNIC) at the National Supercomputer Centre (NSC), which is partially funded by the Swedish Research Council through grant agreement no. 2018-05973. Charles J. R. Williams acknowledges financial support from the UK Natural Environment Research Council within the framework of the SWEET (Super-Warm Early Eocene Temperatures) project (grant no. NE/P01903X/1). Natalie J. Burls acknowledges support from the National Science Foundation (NSF; grant nos. AGS-1844380 and OCN-2002448) and the Alfred P. Sloan Foundation (as a research fellow). Ran Feng acknowledges sponsorship from the U.S. National Science Foundation (grant nos. 1903650 and 1814029). The contributions of Bette L. Otto-Bliesner, Esther C. Brady, and Nan Rosenbloom are based upon work supported by the National Center for Atmospheric Research, which is a major facility sponsored by the NSF under cooperative agreement no. 1852977. The CESM project is primarily supported by the National Science Foundation (NSF). Computing and data storage resources for the CESM and CCSM4 simulations, including the Cheyenne supercomputer (https://doi.org/10.5065/D6RX99HX), were provided by the Computational and Information Systems Laboratory (CISL) at NCAR. Xiangyu Li acknowledges financial support from the National Natural Science Foundation of China (NSFC, grant no. 42005042) and the China Scholarship Council (grant no. 201804910023). The NorESM simulations benefitted from resources provided by UNINETT Sigma2 -the national infrastructure for high-performance computing and data storage in Norway. The work by Anna S. von der Heydt and Michiel L. J. Baatsen was carried out in the framework of the Netherlands Earth System Science Centre (NESSC) program, which is financially supported by the Ministry of Education, Culture and Science (OCW grant no. 024.002.001). Simulations with CCSM4-Utrecht were performed at the SURFsara Dutch national computing facilities and were sponsored by NWO-EW (Netherlands Organisation for Scientific Research, Exact Sciences; project nos. 17189 and 2020.022). Christian Stepanek and Gerrit Lohmann acknowledge computational resources from the Computing and Data Centre of the Alfred Wegener Institute, Helmholtz-Zentrum für Polar-und Meeresforschung. Christian Stepanek and Gerrit Lohmann also acknowledge funding from the Helmholtz Climate Initiative REKLIM and the Alfred Wegener Institute's "Changing Earth-Sustaining our Future" research program. The PRISM4 reconstruction and boundary conditions used in PlioMIP2 were funded by the U.S. Geological Survey Climate and Land Use Change Research and Development Program. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the US Government.
Financial support. This research has been supported by the Swedish Research Council (Vetenskapsrådet; grant nos. 2013-06476 and 2017-04232).
The article processing charges for this open-access publication were covered by Stockholm University.
Review statement. This paper was edited by Laurie Menviel and reviewed by Dipanjan Dey and one anonymous referee.