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 CO₂ 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 mid-Pliocene 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.

1 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₂ rises, including an increase in surface temperature (Xie et al., 2010; 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₂ 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₂ 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₂ 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₂ 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) (Haywood et al., 2013), 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₂ concentrations (Salzmann et al., 2013; Howell et al., 2016) and paleogeography and bathymetry (Otto-Bliesner et al., 2016; 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 (Samakinwa et al., 2020). 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. Zhang et al., 2021), climate sensitivity (Haywood et al., 2020), global monsoons (Q. Zhang et al., 2021), 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 hydrological 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.

2 Data and analytical methods

2.1 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 (Dowsett et al., 2016; Haywood et al., 2016), including soils, lakes, land-ice cover, vegetation, topography, and bathymetry. The CO₂ 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 $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ resolution using bilinear interpolation.

Table 1 PlioMIP2 models used in this study.

2.2 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.,

$$\begin{array}{}\text{(1)}& \mathit{\delta }(\stackrel{\mathrm{‾}}{P}-\stackrel{\mathrm{‾}}{E})\approx \underset{\mathit{\delta }\text{TH}}{\underbrace{-{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{w}}g}}\mathrm{\nabla }\cdot {\int }_{\mathrm{0}}^{p_{\mathrm{s}}}\left({\stackrel{\mathrm{‾}}{V}}_{\mathrm{0}}\mathit{\delta }\stackrel{\mathrm{‾}}{q}\right)\mathrm{d}p}}\underset{\mathit{\delta }\text{MCD}}{\underbrace{-{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{w}}g}}\mathrm{\nabla }\cdot {\int }_{\mathrm{0}}^{p_{\mathrm{s}}}\left(\stackrel{\mathrm{‾}}{q}\mathit{\delta }{\stackrel{\mathrm{‾}}{V}}_{\mathrm{0}}\right)\mathrm{d}p}}+R.\end{array}$$

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:

$$\begin{array}{}\text{(2)}& \left\{\begin{array}{l}{\mathit{V}}_{\mathrm{R}}(\mathit{\lambda },\mathit{\theta },\mathit{\sigma })=u_{\mathrm{R}}(\mathit{\lambda },\mathit{\theta },\mathit{\sigma })\mathit{i}+v_{\mathrm{R}}(\mathit{\lambda },\mathit{\theta },\mathit{\sigma })\mathit{j},\\ {\mathit{V}}_{\mathrm{M}}(\mathit{\lambda },\mathit{\theta },\mathit{\sigma })=v_{\mathrm{M}}(\mathit{\lambda },\mathit{\theta },\mathit{\sigma })\mathit{j}+{\dot{\mathit{\sigma }}}_{\mathrm{M}}(\mathit{\lambda },\mathit{\theta },\mathit{\sigma })\mathit{k},\\ {\mathit{V}}_{\mathrm{Z}}(\mathit{\lambda },\mathit{\theta },\mathit{\sigma })=u_{\mathrm{Z}}(\mathit{\lambda },\mathit{\theta },\mathit{\sigma })\mathit{i}+{\dot{\mathit{\sigma }}}_{\mathrm{Z}}(\mathit{\lambda },\mathit{\theta },\mathit{\sigma })\mathit{k}.\end{array}\right.\end{array}$$

Here, the following continuity equations are satisfied:

$$\begin{array}{}\text{(3)}& \left\{\begin{array}{l}\frac{\mathrm{1}}{\sin \mathit{\theta }}\frac{\partial u_{\mathrm{R}}}{\partial \mathit{\lambda }}+\frac{\mathrm{1}}{\sin \mathit{\theta }}\frac{\partial \left(\sin \mathit{\theta }{v}_{\mathrm{R}}\right)}{\partial \mathit{\theta }}=\mathrm{0},\\ {\displaystyle \frac{\mathrm{1}}{\sin \mathit{\theta }}}{\displaystyle \frac{\partial \left(\sin \mathit{\theta }{v}_{\mathrm{M}}\right)}{\partial \mathit{\theta }}}+{\displaystyle \frac{\partial {\dot{\mathit{\sigma }}}_{\mathrm{M}}}{\partial \mathit{\sigma }}}=\mathrm{0},\\ {\displaystyle \frac{\mathrm{1}}{\sin \mathit{\theta }}}{\displaystyle \frac{\partial u_{\mathrm{Z}}}{\partial \mathit{\lambda }}}+{\displaystyle \frac{\partial {\dot{\mathit{\sigma }}}_{\mathrm{Z}}}{\partial \mathit{\sigma }}}=\mathrm{0}.\end{array}\right.\end{array}$$

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 follows:

$$\begin{array}{}\text{(4)}& \left\{\begin{array}{l}{u}_{\mathrm{R}}=-\frac{\partial R}{\partial \mathit{\theta }},v_{\mathrm{R}}=\frac{\mathrm{1}}{\sin \mathit{\theta }}\frac{\partial R}{\partial \mathit{\lambda }},\\ v_{\mathrm{M}}=-{\displaystyle \frac{\partial H}{\partial \mathit{\sigma }}},{\dot{\mathit{\sigma }}}_{\mathrm{M}}={\displaystyle \frac{\mathrm{1}}{\sin \mathit{\theta }}}{\displaystyle \frac{\partial \left(\sin \mathit{\theta }H\right)}{\partial \mathit{\theta }}},\\ u_{\mathrm{Z}}={\displaystyle \frac{\partial W}{\partial \mathit{\sigma }}},{\dot{\mathit{\sigma }}}_{\mathrm{Z}}=-{\displaystyle \frac{\mathrm{1}}{\sin \mathit{\theta }}}{\displaystyle \frac{\partial W}{\partial \mathit{\lambda }}}.\end{array}\right.\end{array}$$

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,

$$\begin{array}{}\text{(5)}& \mathit{V}={\mathit{V}}_{\mathrm{M}}+{\mathit{V}}_{\mathrm{Z}}+{\mathit{V}}_{\mathrm{R}}\end{array}$$

with the following components:

$$\begin{array}{}\text{(6)}& \left\{\begin{array}{l}u=u_{\mathrm{Z}}+u_{\mathrm{R}}=\frac{\partial W}{\partial \mathit{\sigma }}-\frac{\partial R}{\partial \mathit{\theta }},\\ v=v_{\mathrm{R}}+v_{\mathrm{M}}={\displaystyle \frac{\mathrm{1}}{\sin \mathit{\theta }}}{\displaystyle \frac{\partial R}{\partial \mathit{\lambda }}}-{\displaystyle \frac{\partial H}{\partial \mathit{\sigma }}},\\ \dot{\mathit{\sigma }}={\dot{\mathit{\sigma }}}_{\mathrm{M}}+{\dot{\mathit{\sigma }}}_{\mathrm{Z}}={\displaystyle \frac{\mathrm{1}}{\sin \mathit{\theta }}}{\displaystyle \frac{\partial \left(\sin \mathit{\theta }H\right)}{\partial \mathit{\theta }}}-{\displaystyle \frac{\mathrm{1}}{\sin \mathit{\theta }}}{\displaystyle \frac{\partial W}{\partial \mathit{\lambda }}}.\end{array}\right.\end{array}$$

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 three-pattern 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):

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{\mathrm{1}}{\sin \mathit{\theta }}}{\displaystyle \frac{\partial H}{\partial \mathit{\lambda }}}+{\displaystyle \frac{\mathrm{1}}{\sin \mathit{\theta }}}{\displaystyle \frac{\partial \left(W\sin \mathit{\theta }\right)}{\partial \mathit{\theta }}}+{\displaystyle \frac{\partial R}{\partial \mathit{\sigma }}}=\mathrm{0}.\end{array}$$

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:

$$\begin{array}{}\text{(8)}& \begin{array}{rl}& \mathit{\delta }\text{TH}=\\ & \underset{\mathrm{\Delta }\mathit{\delta }{\text{TH}}_{\mathrm{D}}}{\underbrace{\underset{\mathit{\delta }{\text{TH}}_{\mathrm{D}\mathrm{\_}\mathrm{R}}}{\underbrace{-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }g}}{\int }_{\mathrm{0}}^{p_{\mathrm{s}}}\mathit{\delta }q\mathrm{\nabla }\cdot {\mathit{V}}_{\mathrm{R}\mathrm{0}}\mathrm{d}p}}\underset{\mathit{\delta }{\text{TH}}_{\mathrm{D}\mathrm{\_}\mathrm{Z}}}{\underbrace{-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }g}}{\int }_{\mathrm{0}}^{p_{\mathrm{s}}}\mathit{\delta }q\mathrm{\nabla }\cdot {\mathit{V}}_{\mathrm{Z}\mathrm{0}}\mathrm{d}p}}\underset{\mathit{\delta }{\text{TH}}_{\mathrm{D}\mathrm{\_}\mathrm{M}}}{\underbrace{-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }g}}{\int }_{\mathrm{0}}^{p_{\mathrm{s}}}\mathit{\delta }q\mathrm{\nabla }\cdot {\mathit{V}}_{\mathrm{M}\mathrm{0}}\mathrm{d}p}}}}\\ & \underset{\mathit{\delta }{\text{TH}}_{\mathrm{A}}}{\underbrace{\underset{\mathit{\delta }{\text{TH}}_{\mathrm{A}\mathrm{\_}\mathrm{R}}}{\underbrace{-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }g}}{\int }_{\mathrm{0}}^{p_{\mathrm{s}}}{\mathit{V}}_{\mathrm{R}\mathrm{0}}\cdot \mathrm{\nabla }\mathit{\delta }q\mathrm{d}p}}\underset{\mathit{\delta }{\text{TH}}_{\mathrm{A}\mathrm{\_}\mathrm{Z}}}{\underbrace{-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }g}}{\int }_{\mathrm{0}}^{p_{\mathrm{s}}}{\mathit{V}}_{\mathrm{R}\mathrm{0}}\cdot \mathrm{\nabla }\mathit{\delta }q\mathrm{d}p}}\underset{\mathit{\delta }{\text{TH}}_{\mathrm{A}\mathrm{\_}\mathrm{M}}}{\underbrace{-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }g}}{\int }_{\mathrm{0}}^{p_{\mathrm{s}}}{\mathit{V}}_{\mathrm{M}\mathrm{0}}\cdot \mathrm{\nabla }\mathit{\delta }q\mathrm{d}p}}}},\end{array}\end{array}$$

$$\begin{array}{}\text{(9)}& \begin{array}{rl}& \mathit{\delta }\text{MCD}=\\ & \underset{\mathit{\delta }{\text{MCD}}_{\mathrm{D}}}{\underbrace{\underset{\mathit{\delta }{\text{MCD}}_{\mathrm{D}\mathrm{\_}\mathrm{R}}}{\underbrace{-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }g}}{\int }_{\mathrm{0}}^{p_{\mathrm{s}}}{q}_{\mathrm{0}}\mathrm{\nabla }\cdot \mathit{\delta }{\mathit{V}}_{\mathrm{R}}\mathrm{d}p}}\underset{\mathit{\delta }{\text{MCD}}_{\mathrm{D}\mathrm{\_}\mathrm{Z}}}{\underbrace{-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }g}}{\int }_{\mathrm{0}}^{p_{\mathrm{s}}}{q}_{\mathrm{0}}\mathrm{\nabla }\cdot \mathit{\delta }{\mathit{V}}_{\mathrm{Z}}\mathrm{d}p}}\underset{\mathit{\delta }{\text{MCD}}_{\mathrm{D}\mathrm{\_}\mathrm{M}}}{\underbrace{-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }g}}{\int }_{\mathrm{0}}^{p_{\mathrm{s}}}{q}_{\mathrm{0}}\mathrm{\nabla }\cdot \mathit{\delta }{\mathit{V}}_{\mathrm{M}}\mathrm{d}p}}}}\\ & \underset{\mathit{\delta }{\text{MCD}}_{\mathrm{A}}}{\underbrace{\underset{\mathit{\delta }{\text{MCD}}_{\mathrm{A}\mathrm{\_}\mathrm{R}}}{\underbrace{-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }g}}{\int }_{\mathrm{0}}^{p_{\mathrm{s}}}\mathit{\delta }{\mathit{V}}_{\mathrm{R}}\cdot \mathrm{\nabla }{q}_{\mathrm{0}}\mathrm{d}p}}\underset{\mathit{\delta }{\text{MCD}}_{\mathrm{A}\mathrm{\_}\mathrm{Z}}}{\underbrace{-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }g}}{\int }_{\mathrm{0}}^{p_{\mathrm{s}}}\mathit{\delta }{\mathit{V}}_{\mathrm{Z}}\cdot \mathrm{\nabla }{q}_{\mathrm{0}}\mathrm{d}p}}\underset{\mathit{\delta }{\text{MCD}}_{\mathrm{A}\mathrm{\_}\mathrm{M}}}{\underbrace{-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }g}}{\int }_{\mathrm{0}}^{p_{\mathrm{s}}}\mathit{\delta }{\mathit{V}}_{\mathrm{M}}\cdot \mathrm{\nabla }{q}_{\mathrm{0}}\mathrm{d}p}}}},\end{array}\end{array}$$

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.

3 Changes in hydroclimate during the mid-Pliocene

3.1 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 (Zhang et al., 2019; Li et al., 2020; Feng et al., 2021). Zhang et al. (2016) indicate that the combined influence of SST and CO₂ 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 (Haywood et al., 2020; 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).

Figure 1 The simulated PmE changes in individual PlioMIP2 models. The zonal average of the PmE changes in each model is shown in panel (b) (units: mm d⁻¹).

Figure 2 (a) Changes in the multi-model mean (MMM) PmE for the mid-Pliocene compared with the PI simulation (shading), overlaid by the climatological MMM PmE of the PI simulation (for the contours, a solid line indicates positive values and a dashed line indicates negative values). The red solid curves represent the zero value. (b) The zonal average of the change in PmE, where the shading indicates the interquartile range among models. (Units are mm d⁻¹.)

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₂ 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.

3.2 Previous model–data comparisons of hydrological changes in the PlioMIP2 ensemble

Multi-proxy studies are qualitatively consistent with the results of the PlioMIP2 ensemble (Feng et al., 2021). 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 (Samakinwa et al., 2020).

4 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).

Figure 3 The annual mean changes in the moisture budget components of the mid-Pliocene minus the PI control of the PlioMIP2 multi-model mean, which reflect (a) the thermodynamic term and (b) the estimated change in the thermodynamic term (calculated as the sum of the right-hand side terms in Eq. (8) of the 3P-DGAC decomposition method). The contributions to the change in the thermodynamic effect include the thermodynamic term induced by the climate mean (c) meridional, (d) zonal, and (e) horizontal circulations. The corresponding changes in panels (f–h) and (i–k) are the terms in the components shown in panels (c–e) that are related to the divergent mean flow and change in moisture advection, respectively. Stippling indicates regions where at least 10 of 13 simulations in the model group agree on the sign of the ensemble mean. (Units are mm d⁻¹.)

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 divergence 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 consider 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).

Figure 4 The annual mean changes in the moisture budget components of the mid-Pliocene minus the PI control of the PlioMIP2 multi-model mean, which reflect (a) the dynamic term and (b) the estimated change in the dynamic term (calculated as the sum of the right-hand side terms in Eq. 9 of the 3P-DGAC decomposition method). The contributions to the changes in the dynamic effect include the dynamic term induced by the anomalous (c) meridional, (c) zonal, and (d) horizontal circulations. The corresponding changes in panels (f–h) and (i–k) are the terms in the components shown in panels (c–e) that are related to the change in divergent flow and mean moisture advection, respectively. Stippling indicates regions where at least 10 of 13 simulations in the model group agree on the sign of the ensemble mean. (Units are mm d⁻¹.)

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.

Figure 5 The estimated annual mean changes in the PmE (calculated as the sum of the dynamic term and thermodynamic term) of the mid-Pliocene minus the PI control of the PlioMIP2 multi-model mean (units: mm d⁻¹).

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 humidity 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).

Figure 6 (a) The zonal mean of meridional moisture transport anomalies induced by thermodynamic and dynamic effects for the mid-Pliocene compared with PI simulations. Panels (b) and (c) show the thermodynamic and dynamic terms in panel (a) that are induced by the meridional wind V_M and V_R decomposed from the 3P-DGAC method, respectively. Here, the tropical region is defined as the region between 10^∘ S and 10^∘ N (marked as an orange band), whereas the subtropical region refers to 10–30^∘ N and 10–30^∘ S (marked as a cyan band). The shading indicates 1 standard deviation of individual models departing from the MMM. (Units are g kg⁻¹ m s⁻¹.)

Figure 7 Change in (a) MMM SST (shading; in ^∘C), and (b) specific humidity (shading; in kg kg⁻¹) overlaid by its climate mean for the PI simulation. Panels (c) and (d) show the MMM divergence of the meridional V_M and zonal wind V_Z fields decomposed from the 3P-DGAC method at the 850 hPa level for the PI simulation (in 10⁻⁶ s⁻¹). The circles in panel (a) are the anomalies of reconstructed SST (McClymont et al., 2020) from the alkenone-derived ${\mathrm{U}}_{\mathrm{37}}^{{\mathrm{K}}^{\prime }}$ index (Prahl and Wakeham, 1987) and foraminifera calcite $\mathrm{Mg}/\mathrm{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.

5 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.

5.1 Changes in specific humidity

Figure 7a shows the changes in MMM SST superimposed on the reconstructed SST anomalies (McClymont et al., 2020). 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 (Li et al., 2020). However, Z. Zhang et al. (2021) 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 wind–evaporation–SST feedback (Xie et al., 2010). These SST warming patterns are consistent with current studies (Haywood et al., 2020; 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 $\mathrm{\nabla }\cdot \stackrel{\mathrm{‾}}{\mathit{V}}>\mathrm{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., $-\frac{\mathrm{1}}{\mathit{\rho }g}{\int }_{\mathrm{0}}^{p_{\mathrm{s}}}\mathit{\delta }q\mathrm{\nabla }\cdot {\mathit{V}}_{\mathrm{M}\mathrm{0}}\mathrm{d}p$) to the thermodynamic effect (Fig. 3f) and, hence, changes in PmE.

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., $-\mathrm{\nabla }\cdot {\stackrel{\mathrm{‾}}{\mathit{V}}}_{{\mathrm{Z}}_{\mathrm{0}}}>\mathrm{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., $-\frac{\mathrm{1}}{\mathit{\rho }g}{\int }_{\mathrm{0}}^{p_{\mathrm{s}}}\mathit{\delta }q\mathrm{\nabla }\cdot {\mathit{V}}_{\mathrm{Z}\mathrm{0}}\mathrm{d}p)$ 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).

Figure 8 (a) Changes in the annual mean mass stream function (MSF; shading; in 10¹⁰ kg s⁻¹) of meridional circulation for mid-Pliocene with respect to the PI simulation, overlaid by the climate mean MSF for the PI simulation (contours). The meridional wind V_M is decomposed from the 3P-DGAC method. Solid curves indicate positive values, and dashed curves indicate negative values. Stippling indicates regions where at least 10 of 13 simulations in the model group agree on the sign of the ensemble mean. (b) Changes in the annual mean intensities (in 10¹⁰ kg s⁻¹) of meridional circulation in the NH and SH. (c) The latitudes of the center of annual mean precipitation (Pcent) between 20^∘ S and 20^∘ N (in ^∘). (d) Hemispheric asymmetry (NH minus SH) of the energy flux into the atmosphere (in W m⁻²). (e) Changes in the integrated atmospheric meridional heat transport across the Equator (in PW).

5.2 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 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 (Feng et al., 2020).

Figure 9 The changes in meridional wind V_M decomposed from the 3P-DGAC method at the 850 hPa level (vectors; in m s⁻¹), overlaid by its divergent circulations (shading; in 10⁻⁶ s⁻¹). Only vectors where at least 10 of 13 simulations in the model group agree on the sign of the ensemble mean are shown here.

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 (Zhang et al., 2013; 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⁻² more net radiation than the SH (Fig. 8d), which could induce an increased cross-equatorial 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.

Figure 10 Changes in the annual mean of the (a) intensities of the Pacific Walker circulation (PWC; in 10¹⁰ kg s⁻¹), (b) longitude of the PWC cell center (in ^∘), (c) dSST (in K), and (d) dSLP (in hPa) for mid-Pliocene simulations compared with the PI simulation. Here, the PWC intensity is defined as the vertically integrated zonal mass stream function (ZMS; Bayr et al., 2014; Schwendike et al., 2014) averaged in the equatorial Pacific (140^∘ E–120^∘ W), and the location of the PWC cell center is the longitude of the maximum ZMS. 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).

5.3 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, respectively. 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 east–west 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 Pacific and strengthened descending motion on the western Indian Ocean, indicating a westward expansion of the PWC (Fig. 11b).

Figure 11 (a) Changes in ZMS (shading; in 10¹⁰ kg s⁻¹) averaged between 10^∘ S and 10^∘ N for the mid-Pliocene with respect to the PI simulation, overlaid by the climate mean ZMS for the PI simulation (contours). The zonal wind V_Z is used to calculate ZMS, which is decomposed from the 3P-DGAC method. The contours represent the climate mean ZMS for the PI simulation. Solid curves indicate a positive value, and dashed curves show a negative value. Stippling indicates regions where at least 10 of 13 simulations in the model group agree on the sign of the ensemble mean. Panel (b) shows the vertically integrated ZMS from panel (a). The gray and pink shading indicates 1 standard deviation of individual models' departure from the MMM mean of the MSF for the PI and mid-Pliocene simulations, respectively.

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 ($\stackrel{\mathrm{‾}}{q}>\mathrm{0}$) in the lower troposphere, they can trigger a negative/positive contribution from the δMCD_{D_Z} term to changes in PmE (Fig. 4g).

Figure 12 Changes in the potential function of zonal wind V_Z at (a) 850 hPa and (b) 250 hPa (shading; in kg s⁻¹), corresponding to the divergent mode of the wind field (vectors; in kg m⁻¹ s⁻¹). V_Z is decomposed from the 3P-DGAC method. The contours represent the climate mean of the potential function for the PI simulation. Solid curves indicate positive values, and dashed curves indicate negative values. Red solid curves represent zero values. The vectors and stippled regions are where at least 10 of 13 simulations in the model group agree on the ensemble mean.

6 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 (Oldeman et al., 2021). 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 middle–high latitudes.

Note that the global temperature during mid-Pliocene is controlled by the combined effects of boundary conditions (e.g., CO₂ level, vegetation, and topography) (Haywood et al., 2016). 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.