On the importance of moisture conveyor belts from the tropical eastern Pacific for wetter conditions in the Atacama Desert during the mid-Pliocene

Geomorphic and sedimentologic data indicate that the climate of today's hyper-arid Atacama Desert (northern Chile) was more humid during the mid-Pliocene to Late Pliocene. The processes, however, leading to increased rainfall in this period are largely unknown. To uncover these processes we use both global and regional kilometre-scale model experiments for the mid-Pliocene (3.2 Ma). We found that the PMIP4–CMIP6 (Paleoclimate Modelling Intercomparison Project–Coupled Model Intercomparison Project) model CESM2 (Community Earth System Model 2) and the regional model WRF (Weather Research and Forecasting) used in our study simulate more rainfall in the Atacama Desert for the mid-Pliocene in accordance with proxy data, mainly due to stronger extreme rainfall events in winter. Case studies reveal that these extreme winter rainfall events during the mid-Pliocene are associated with strong moisture conveyor belts (MCBs) originating in the tropical eastern Pacific. For present-day conditions, in contrast, our simulations suggest that the moisture fluxes rather arise from the subtropical Pacific region and are much weaker. A clustering approach reveals systematic differences between the moisture fluxes in the present-day and mid-Pliocene climates, both in strength and origins. The two mid-Pliocene clusters representing tropical MCBs and occurring less than 1 d annually on average produce more rainfall in the hyper-arid core of the Atacama Desert south of 20^∘ S than what is simulated for the entire present-day period. We thus conclude that MCBs are mainly responsible for enhanced rainfall during the mid-Pliocene. There is also a strong sea-surface temperature (SST) increase in the tropical eastern Pacific and along the Atacama coast for the mid-Pliocene. It suggests that a warmer ocean in combination with stronger mid-tropospheric troughs is beneficial for the development of MCBs leading to more extreme rainfall in a +3 ^∘C warmer world like in the mid-Pliocene.

1 Introduction

The Atacama Desert in northern Chile (19–26^∘ S) is considered to be one of the driest deserts on Earth, with a mean annual rainfall of only a few millimetres in the hyper-arid core and along the Coastal Cordillera (e.g. Houston, 2006; Reyers et al., 2021). Due to its special geographic position, the aridity in the Atacama Desert is long-lasting, with its onset likely dating back to the Oligocene (ca. 33.9–23 Ma), thus making it also one of the oldest deserts (e.g. Dunai et al., 2005; Evenstar et al., 2017). Paleoclimate records and geomorphic studies from various sites of the Atacama Desert, however, indicate variations in the overall aridity over time, which was repeatedly punctuated by more humid phases on millennial and orbital timescales (Ritter et al., 2019; Diederich et al., 2020). Based on geomorphic and sedimentologic data, Amundson et al. (2012) and Hartley and Chong (2002) postulated that nowadays dominantly hyper-arid conditions commenced during the Late Pliocene (ca. 3–2 Ma). This shift to overall hyper-aridity postdates a period of more humid or even semi-arid conditions in the Atacama Desert during the mid-Pliocene to Late Pliocene, which is illustrated in the widespread occurrence of lake, salar and fluvial deposits in many parts of the desert (e.g. Gaupp et al., 1999; Kirk-Lawlor et al., 2013; Jordan et al., 2014; Evenstar et al., 2016; Ritter et al., 2018; Vásquez et al., 2018). While such proxy records provide a reliable picture of the timing and magnitude of past pluvial phases, the mechanisms controlling these climate shifts and in particular the involved atmospheric processes are hard to derive from them and are thus not well understood. Intervals of increased rainfall in the southern Atacama Desert are mostly attributed to a northward displacement of mid-latitudinal westerlies and accompanied extratropical winter cyclones (Vuille and Ammann, 1997; Stuut and Lamy, 2004; Bartz et al., 2020), which suggest a southwestern moisture source. Past rainfall variations in the northern Atacama Desert and the Andes have been linked to latitudinal shifts of the extratropical westerlies in the Southern Hemisphere (Amidon et al., 2017). For the warmer mid-Pliocene climate, the multi-model mean of the PlioMIP2 (Pliocene Model Intercomparison Project) models for instance indicates that the Hadley cell was shifted northward and the Walker circulation was shifted westward (e.g. Han et al., 2021). Also, cutoff lows as seen in March 2015 have been proposed as a possible mechanism for wetter conditions in the past (Jordan et al., 2019). Wetter conditions during the warm Pliocene are furthermore attributed to a mean state more similar to today's El Niño conditions, and according to this the onset of the hyper-aridity concurrently occurred with the change to the present-day El Niño–Southern Oscillation (ENSO) variability (Amundson et al., 2012). Hartley and Chong (2002) concluded that the shift from semi-arid to hyper-arid conditions in the Atacama Desert was controlled by global cooling rather than by local factors, like the cold Humboldt Current.

Another perspective on potential drivers for wetter paleoclimate episodes can be inferred from the mechanisms involved in the March 2015 extreme rainfall event that caused severe damage in the Atacama Desert. The event was associated with a mid-tropospheric cutoff low and anomalously warm tropical sea-surface temperatures (SSTs) along the Atacama coast, paired with transport of large amounts of water vapour from the tropical eastern Pacific to the Atacama Desert at the foreside of the low-pressure system (Bozkurt et al., 2016). Based on the characteristic isotopic composition of rainwater (e.g. Herrera and Custodio, 2014) from the March 2015 event, Jordan et al. (2019) proposed that the processes involved in this event might also play an important role in increased paleoclimate rainfall in the Atacama Desert.

Böhm et al. (2021, hereafter B2021) systematically investigated moisture conveyor belts (MCBs) in the Atacama Desert for the present-day climate that also occurred during the March 2015 event. MCBs are elongated bands of strong poleward water vapour fluxes in the free troposphere (e.g. Knippertz and Martin, 2007), which often originate in tropical latitudes and occur in connection with troughs or cutoff lows. Based on different datasets for the recent past, B2021 indicate that moisture transport associated with MCBs is decoupled from the maritime boundary layer. Despite their rareness of about four regional MCBs per year which make landfall in the Atacama Desert, MCBs produce 40 %–80 % of the total rainfall in the hyper-arid core and along the Coastal Cordillera for modern climate conditions (B2021).

The aim of the present study is to assess to what extent MCBs act as drivers of rainfall activity in the Atacama Desert on geological timescales, specifically during the relatively wet and warm mid-Pliocene period. If MCBs played an important role in that period, this would imply that in addition to the previously suggested regions southwest or east of the Atacama Desert (Stuut and Lamy, 2017; Bartz et al., 2020; Amidon et al., 2017) also the tropical southeastern Pacific northwest of the desert could be a potential moisture source for increased humidity in the mid-Pliocene, like assessments of the regional rainfall under present-day climate suggest (Bozkurt et al., 2016; Jordan et al., 2019; Böhm et al., 2021). However, the past and present constellations of the global atmospheric and oceanic circulations are substantially different. The Atacama Desert is characterized by a complex topography, including a steep coastal cliff and the western slopes of the Andes, as shown in Fig. 1. As a consequence, strong rainfall gradients occur, in particular in the west-to-east direction (Houston, 2006; Reyers et al., 2021). Due to the rather coarse horizontal resolution of global climate models, the kilometre-scale changes in the orography and rainfall patterns cannot be adequately simulated by design (e.g. Ludwig et al., 2019; Fiedler et al., 2020). We therefore perform kilometre-scale simulations for a limited-area domain over the Atacama Desert for the mid-Pliocene at around 3.2 Ma. Section 2 gives an overview of the datasets and methods used in our study. Section 3 consists of two subsections, where we present large-scale and regional paleoclimate changes as obtained by global and downscaled PMIP4–CMIP6 (Paleoclimate Modelling Intercomparison Project–Coupled Model Intercomparison Project) models (Sect. 3.1) and a potential key driver for the detected rainfall changes (Sect. 3.2). A summary and discussion close this paper (Sect. 4).

Figure 1 Map showing the outer (50 km×50 km resolution) and inner WRF model domain (red box, 10 km×10 km resolution) and the topographic heights (shading, in m).

2 Data and methods

2.1 Model simulations

We use global climate model (GCM) output from the fourth phase of the Paleoclimate Modelling Intercomparison Project, which is endorsed by the sixth phase of the Coupled Model Intercomparison Project (PMIP4–CMIP6; Kageyama et al., 2018). Here, we examine the output of a mid-Pliocene warm-period experiment (3.2 Ma), as part of the Pliocene Model Intercomparison Project Phase 2 (PlioMIP2; Haywood et al., 2016). In the PlioMIP2 simulations, the CO₂ concentration is set to 400 ppm, while other greenhouse-gas concentrations and the orbital parameters are as for the pre-industrial period (1850) to be consistent with the setup of PlioMIP2 experiments. Paleoenvironmental reconstructions from the Pliocene Research, Interpretation and Synoptic Mapping project Phase 4 (PRISM4) are used as paleogeography boundary conditions, e.g. topography, bathymetry, ice sheets and vegetation cover. A significant feature of these reconstructions is the strong reduction of the Greenland and Antarctic ice sheets compared to today. Dowsett et al. (2016) describes PRISM4 in detail.

From PlioMIP2 we select simulations from the Community Earth System Model 2 (CESM2) of the National Center for Atmospheric Research (NCAR) (Danabasoglu et al., 2020) for our study. This choice has been made based on the availability of boundary data and the validation results of GCMs participating in PlioMIP2 assessed previously (e.g. Feng et al., 2020; Pontes et al., 2020; Han et al., 2021). In addition to CESM2, two other PlioMIP2–CMIP6 models are analysed, namely IPSL-CM6A-LR (Boucher et al., 2018) from the Institut Pierre-Simon Laplace and GISS-E2.1-G (Bauer et al., 2020) from NASA. Only for these models are the historical and the mid-Pliocene experiment available. Our PlioMIP2 intercomparison for the Atacama Desert revealed the following:

• i.

  CESM2 reproduces the increased rainfall for the Atacama region compared to the present day in accordance with proxy data (Sect. 3.1).

• ii.

  CESM2 shows an expected SST increase for the mid-Pliocene, which qualitatively agrees with reconstructions (Sect. 3.1).

This makes the CESM2 output the most useful amongst the three GCMs from PlioMIP2 to perform the high-resolution simulations for the Atacama Desert.

For the atmosphere CESM2 has a horizontal resolution of $\mathrm{1.25}{}^{\circ }\times \mathrm{0.9}{}^{\circ }$ in the longitudinal and latitudinal direction, respectively. The horizontal resolution of the ocean is constant in the zonal direction (1.125^∘) and varies from 0.27 to 0.64^∘ in the meridional direction. We select a 30-year period in the mid-Pliocene experiment of CESM2 that includes a strong ENSO and rainfall variability (model years 231–260) to represent internal variability. To quantify changes relative to the present day, we further use a 30-year time slice of the CMIP6 historical experiment (years 1985–2014) of CESM2.

The output of CESM2 is used as initialization and boundary data for simulations with the Weather Research and Forecasting model version 3.9 (WRF; Skamarock et al., 2008). The lateral boundary fields are 6-hourly updated. We yield a horizontal resolution of 10 km×10 km via a double one-way nesting, with the first nesting from $\mathrm{1.25}{}^{\circ }\times \mathrm{0.9}{}^{\circ }$ to 50 km×50 km and the second nesting from 50 km×50 km to 10 km×10 km. The outer and inner model domains for our WRF simulations are depicted in Fig. 1. For our analysis we use the output of both nesting steps (i.e. WRF simulations with 50 and 10 km, respectively).

Reyers et al. (2021) demonstrate that the WRF model reproduces present-day rainfall in the Atacama Desert when driven with reanalysis data, despite the complexity of the terrain. We therefore use the same model setup in terms of the physical parameterizations for this study on the mid-Pliocene climate. The WRF simulation uses prescribed time-varying SSTs from CESM2 that are updated every 24 h. The same 30-year mid-Pliocene (mP) and historical (hist) periods as used for the analysis of the global CESM2 model are simulated with WRF (hereafter WRF_mP and WRF_hist). For the WRF_mP simulation the greenhouse-gas concentrations and the orbital parameters are prescribed as in the global PlioMIP2 simulations (see above). Paleorecords reveal that the land cover strongly changed during the Pliocene on global scales, but a reconstruction for South America is difficult due to the sparse data coverage (Dowsett et al., 2016). As a consequence, reconstructed land cover for the Atacama Desert is implausible. To avoid uncertainties arising from artefacts in reconstructed land-cover changes and for a better comparability of our WRF simulations, we therefore prescribe the same present-day vegetation mask in all WRF simulations. For an additional evaluation of the WRF simulations for the Atacama Desert, rainfall of WRF_hist is compared to the present-day validation simulation with WRF from Reyers et al. (2021), which was driven by ERA-Interim reanalysis data (WRF_era).

2.2 Analysis strategy

For the assessment of the dynamical processes involved in rainfall in the Atacama Desert we first quantify atmospheric moisture fluxes from the WRF simulations. This is done by computing the integrated water vapour flux (IWVF), which consists of a zonal (IWVF_u) and a meridional (IWVF_v) component:

$$\begin{array}{}\text{(1)}& {\displaystyle }{\text{IWVF}}_u=\underset{z\mathrm{1}}{\overset{z\mathrm{2}}{\int }}uq\mathrm{d}z,\text{(2)}& {\displaystyle }{\text{IWVF}}_v=\underset{z\mathrm{1}}{\overset{z\mathrm{2}}{\int }}vq\mathrm{d}z,\end{array}$$

where q is the specific humidity, z is the height, and u and v are the zonal and the meridional wind components, respectively. B2021 demonstrate that MCBs are associated with elevated moisture transport decoupled from the maritime boundary layer. To take this into account we compute IWVF for low levels between 0 and 2000 m above ground level (a.g.l.) and for upper levels above 2000 m a.g.l. separately.

For the classification of MCBs in our study area, we apply machine learning techniques to the IWVF over the offshore region of western central South America. Due to the rareness of MCBs in the Atacama Desert for the present-day climate (B2021), we choose a method that can detect outliers. This is a combination of self-organizing maps (SOMs, Kohonen, 2001) and k-means clustering (e.g. Hartigan and Wong, 1979). This combined method has already been applied successfully in climate research (e.g. Ohba et al., 2016). SOMs are a machine learning technique often used for pattern recognition. In this technique, high-dimensional data are projected onto a visually interpretable two-dimensional grid by mapping similar input vectors close to each other. As input vectors we use 12-hourly IWVF patterns as diagnosed from our WRF simulations. IWVF values below a certain threshold are set to zero before the SOM clustering. This threshold-based methodology is inspired by approaches which are used to identify low-level atmospheric rivers (e.g. Rutz et al., 2014) and ensures that aside from the strength of the MCBs also information about their shape is provided for the clustering. For WRF_hist an IWVF threshold of 250 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$ is chosen. Due to the fact that during the warmer mid-Pliocene the atmosphere is generally moister, the threshold for WRF_mP is set to a higher value of 350 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$. Sensitivity tests indicate that slight differences in the thresholds yield similar results from the clustering, although the choice of these threshold values is somewhat arbitrary. The input vectors are projected on a 7×7 grid SOM, thus resulting in 49 IWVF classes. To further reduce these data, the 49 IWVF classes are processed via k-means clustering, which is an iterative algorithm that obtains an optimal partition of the input data (here, the 49 IWVF classes from the SOM) into k clusters. The algorithm computes so-called cluster centroids, and an optimal solution is found when the squared distances between the cluster members and their respective cluster centroids are minimized. The resulting cluster centroids represent the final clusters. For a better comparability, we aim for the same number of final clusters for WRF_hist and WRF_mP. Nine final clusters are optimal, since the squared distances between the centroids converge towards a minimum value for more than nine clusters for both WRF_hist and WRF_mP. Note that the order of clusters is arbitrary.

For the identification of mid-tropospheric troughs off the Atacama coast we adopt the automated identification of Knippertz (2004). A grid point P is defined as a trough point, where the west-to-east gradient in the 500 hPa geopotential height exceeds a threshold of 25 gpm (geopotential metres) per 10^∘. Here, we apply the detection algorithm to the CESM2 data instead of the regionally limited WRF output since a larger spatial coverage is needed for identifying synoptic-scale troughs. The 500 hPa geopotential height is interpolated to a regular $\mathrm{2}{}^{\circ }\times \mathrm{2}{}^{\circ }$ horizontal grid first. The interpolation is done for consistency with the coarse reanalysis data used in the original development of the detection algorithm.

3 Results

3.1 Mid-Pliocene against present-day climate

The Pliocene is considered to be the most recent period in Earth history with globally sustained temperatures higher by several degrees Celsius and thus a potential analogue for the conditions projected by CMIP6 scenarios for the future. During the stable mid-Pliocene period between 3.264 and 3.025 Ma the global annual mean temperature is, based on climate models, assumed to be up to more than 3 ^∘C warmer than the present day (Haywood et al., 2013; Haywood and Valdes, 2014). The model results are supported by proxy data indicating a global SST anomaly for the mid-Pliocene vs. pre-industrial levels of 2.3 and 3.2–3.4 ^∘C based on foraminifera $\mathrm{Mg}/\mathrm{Ca}$ and alkenones or alkenones only, respectively (McClymont et al., 2020). Specifically in the upwelling regions at the Peruvian margin, Dekens et al. (2007) reconstructed a Pliocene–modern SST change of 2.9 ^∘C. Figure 2a shows the mean annual near-surface air temperature changes for the mid-Pliocene against the present day (historical experiment) from CESM2. A strong polar amplification is visible over Antarctica, with mid-Pliocene temperatures being partly more than 10 ^∘C warmer than for present-day conditions. A particular strong warming is also revealed in the tropical eastern Pacific and along the western coast of South America, which is associated with changes in the SST in the tropical eastern Pacific and in the upwelling zone off the Chilean coast. Here, the SST is simulated to be up to 5 ^∘C warmer during the mid-Pliocene (Fig. 2b). This SST difference pattern qualitatively agrees with the PRISM reconstructions for this area (Dowsett et al., 2013; see their Fig. 5). A strong precipitation increase is found for the Altiplano region, where mean annual rainfall in the mid-Pliocene is up to 400 mm higher than in the historical experiment (Fig. 2c). More rainfall is also simulated for the Atacama Desert. These results for more rainfall are broadly consistent with proxy records for the wetter conditions in the mid-Pliocene compared to present-day conditions in the region, as listed in Table 1.

Figure 2 Changes (mid-Pliocene minus historical) for mean annual (a) near-surface air temperature (at 2 m above ground in ^∘C), (b) sea-surface temperature (in ^∘C) and (c) rainfall (in mm yr⁻¹) as simulated by the global CESM2 model.

Table 1 Proxy data for wetter conditions in the region of the Atacama Desert that fall into the mid-Pliocene.

In contrast, IPSL-CM6A-LR and GISS-E2.1-G simulate drier conditions in the Atacama Desert for the mid-Pliocene compared to the present day (see Figs. S1c and S2c in the Supplement). Both models also exhibit a strong polar amplification like CESM2, but the warming in the Southern Hemisphere mid-latitudes and tropics is much weaker when compared to CESM2 (Figs. S1a and S2a). In particular the atmospheric and surface temperature increase in the tropical eastern Pacific and off the Atacama coast is less pronounced in these two models (Figs. S1a and b and S2a and b) compared to CESM2. Hence, amongst the three available PlioMIP2 models the CESM2 model has the best agreement with proxy data regarding both SST and rainfall in the study area.

We evaluated the rainfall from WRF_hist against a WRF simulation that downscales the ERA5 reanalysis for the same domain and spatial resolution (WRF_era). There are quantitative differences in rainfall, but the aridity is overall satisfyingly reproduced by the WRF simulation that used data from the historical simulation of CESM2 at the lateral boundaries (WRF_hist). Specifically, the spatial patterns and the seasonal cycle of rainfall are qualitatively captured by WRF_hist (compare Fig. 3f–j with Fig. 3a–e). Both WRF_era and WRF_hist show the present-day hyper-aridity in the Atacama Desert, with mean annual rainfall of less than 2 mm yr⁻¹ north of 22^∘ S. Rainfall slightly increases towards the south, but a particular strong west-to-east gradient related to the topographic characteristics of the region is visible in the northeastern part of the model domain. The rainfall at the western slopes of the Andes is strongest in summer (DJF, December–January–February), whereas strongest rainfall in the hyper-arid core of the Atacama Desert is simulated for winter (JJA, June–July–August). Annual and seasonal rainfall amounts tend to be regionally overestimated by WRF_hist against WRF_era, but the hyper-aridity with only a few millimetres of rainfall per year is well simulated (Fig. 3f–j). We therefore conclude that the WRF simulations using CESM2 as boundary conditions are suitable for our research interest.

Figure 3 Present-day mean annual and seasonal rainfall as simulated by (a–e) WRF_era and (f–j) WRF_hist and (k–o) mid-Pliocene changes (WRF_mP minus WRF_hist) in the mean annual and seasonal rainfall. Values for (a, f, k) the whole year, (b, g, l) DJF (December–January–February), (c, h, m) MAM (March–April–May), (d, i, n) JJA (June–July–August) and (f, j, o) SON (September–October–November). All results are shown for the simulations with 10 km resolution, and values are given in millimetres per year and millimetres per season, respectively.

For the mid-Pliocene (WRF_mP, Fig. 3k), our simulations yield more humid conditions, with nearly a doubling in the annual rainfall for most parts of the desert south of 20^∘ S compared to WRF_hist (compare Fig. 3k and f). This mid-Pliocene rainfall increase in the hyper-arid core of the Atacama Desert is mainly due to more winter rainfall (June–August, Fig. 3n). We therefore focus on this season in the remainder of our study and first test the two following hypotheses:

• i.

  The intensity of rain days is similar in both periods, but there are more rain days in WRF_mP compared to WRF_hist.

• ii.

  Individual rain days are more intensive in WRF_mP than in WRF_hist.

The seasonal cycle of rainfall shifts from winter-dominated in the south to summer-dominated in the north of the Atacama Desert (Houston, 2006; Reyers et al., 2021). In order to account for these differences, we split the hyper-arid core into a northern and southern region (see black boxes in Fig. 3n), which has also been done similarly in previous studies (e.g. Böhm et al., 2020). To quantify the cause of the winter rainfall increase (see Fig. 3n), we compute the spatial means of daily rainfall over these two regions. Then we rank the spatially averaged daily rainfall amounts according to their magnitude and display the results in percentile–percentile plots for the mid-Pliocene against the present-day simulation (Fig. 4). For both the northern and the southern boxes, there is no increase in the number of rainfall events in WRF_mP when compared to WRF_hist, and appreciable rainfall amounts are rare in both simulations. Instead, rain days at the upper end of the distribution are stronger during the mid-Pliocene. For example, for the northern (southern) box daily rainfall is approx. 12 mm (55 mm) per grid point on average for the strongest event in WRF_mP but only 5 mm (21 mm) for the strongest event in WRF_hist. Hence, our simulations suggest that the more humid conditions during the mid-Pliocene is due to rain days with more intense rainfall events by a factor of about 2, not due to more rain days.

Figure 4 Percentile–percentile plots for daily winter rainfall (in mm d⁻¹) in WRF_mP (y axis; midPlio) against WRF_hist (x axis; HIST) as spatial means over (a) the northern box and (b) the southern box, as marked in Fig. 3n.

3.2 Potential drivers for stronger rainfall events in the mid-Pliocene

We assess physical processes, potentially driving the stronger winter rainfall events in the mid-Pliocene. To that end, we conduct several case studies of days with strong winter rainfall events in the northern and southern Atacama Desert (see black boxes in Fig. 3n). Figure 5 shows the rainfall patterns of exemplary top winter rainfall events as simulated by WRF_mP. It turned out that 5 out of the top 10 rain days in the northern box are also among the top days in the southern box. As an example, two of these events are displayed in Fig. 5a and b. The events are associated with elongated bands of strong upper-level IWVF (above 2000 m a.g.l.), which have the typical properties of MCBs, as identified in B2021 for present-day conditions. Further, the MCBs originate in the tropical eastern Pacific region and transport moisture along the coastal offshore region towards the Atacama Desert. Strong rainfall events which are restricted to the northern Atacama Desert are in some cases associated with relatively weaker MCBs but with a similar spatial pattern with a northwestern-to-southeastern moisture flux (Fig. 5c and d). A rather zonally oriented moisture flux is found for one top rainfall events which is restricted to the southern Atacama Desert (Fig. 5e). The majority of the southern events, however, also coincide with strong northwesterly MCBs with an origin in the tropical eastern Pacific, e.g. shown in Fig. 5f. These results indicate that for both the southern and northern Atacama Desert, MCBs and thus the same physical processes are responsible for the wetter conditions in the mid-Pliocene. In contrast, upper-level IWVFs during or prior to the top 10 events in WRF_hist are weaker (Fig. 6). Additionally, the IWVF patterns associated with the WRF_hist events more often show a rather zonal orientation and thus originate from the subtropical Pacific (Fig. 6b, d, e and f).

Figure 5 Daily rainfall over land (green–blue shading, in mm d⁻¹) during top winter events in WRF_mP (obtained for the inner model domain with 10 km resolution) and the maximum upper-level IWVF (above 2000 m a.g.l.) over ocean (arrows and red–blue shading, in $\mathrm{kg}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$) on the day of or the day prior to the rainfall event (obtained from the outer model domain with 50 km resolution). (a, b) Events which are among the top events of both the northern and southern Atacama Desert. (c, d) Events which are only among the top events of the northern Atacama Desert. (e, f) Events which are only among the top events of the southern Atacama Desert.

Figure 6 As Fig. 5 but for the top winter events in WRF_hist.

Interestingly, low-level IWVF (below 2000 m a.g.l.) is negligible during all northern and southern top rainfall events in both WRF_mP and WRF_hist (Figs. S3 and S4 in the Supplement), consistent with present-day cases (B2021). Upper-level moisture transport via MCBs into the study area is therefore largely decoupled from the maritime boundary layer in both present-day and mid-Pliocene conditions.

The MCBs identified in the case studies for WRF_mP often have a northwestern-to-southeastern orientation (Fig. 5a, b, d and f) and are therefore associated with a southward moisture transport from the tropics. Figure 7 shows the probability density function of 12-hourly southward upper-level IWV mass fluxes integrated between 80 and 70^∘ W along the 20^∘ S transect (see yellow line in Fig. 5a) in WRF_mP and WRF_hist. The extremely strong southward mass fluxes are more frequent in WRF_mP than in WRF_hist (inset in Fig. 7). For example, moisture mass fluxes of approx. 400×10⁶ kg s⁻¹ occur more than 15 times in 30 Southern Hemisphere winters of WRF_mP but are rare in WRF_hist. Additionally, fluxes of more than 600×10⁶ kg s⁻¹ lacking in WRF_hist regularly occur in WRF_mP. These findings indicate that there are cases in WRF_mP when the large moisture reservoir over the tropical eastern Pacific is tapped via MCBs, advecting moist air to the Atacama Desert.

Figure 7 Probability density function (y axis, in total number of dates) of 12-hourly southward upper-level (above 2000 m a.g.l.) IWV mass fluxes (x axis, in 10⁶ kg s⁻¹) integrated between 80 and 70^∘ W at 20^∘ S (see yellow line in Fig. 5a) in WRF_mP and WRF_hist (outer model domain with 50 km resolution) for Southern Hemisphere winter. The inset zooms into the results above 200×10⁶ kg s⁻¹.

The case studies (Figs. 5 and 6) illustrate the key role of MCBs in increased rainfall in the mid-Pliocene. The question arises as to what extent the characteristics of the MCBs during the mid-Pliocene were systematically different from those under present-day conditions. We perform an objective clustering of upper-level wintertime IWVF (Sect. 2.2) to quantify potential differences. The nine final IWVF clusters for WRF_mP are displayed in Fig. 8. The cluster with the highest frequency of occurrence is cluster 8, which represents cases without MCBs and occurs for 79 d per winter on average. Instead, this cluster includes days with easterly IWVF, which in some cases is associated with deep convection over the Altiplano plateau (Reyers et al., 2021). Consequently, a high rainfall amount in the northeastern part of the study area is found for this cluster. Cluster 1, 2, 5, 7 and 9 include more zonally oriented IWVF with landfall in central or southern Chile with origins further towards the south of the Pacific, e.g. when compared to cluster 4. Also cluster 3 is associated with moisture advection from regions south of 20^∘ S but includes also some transport from further north with an axis more tilted towards meridional directions compared to cluster 1, 2, 5, 7 and 9. These clusters with dominant moisture advection from regions south of 20^∘ S play a minor role in winter rainfall in the Atacama Desert (see red–blue shading in Fig. 8). In contrast, IWVF cluster 4 and 6, which occur approximately every second winter in WRF_mP, produce the largest rain amount during winter in the hyper-arid core of the Atacama Desert. These IWVF clusters represent strong MCBs originating in the tropical eastern Pacific region, which make landfall at the coastal region of the Atacama Desert. In WRF_hist, the moisture fluxes in the IWVF cluster are generally weaker than in WRF_mP (Fig. 9). As for WRF_mP, the more zonally oriented IWVF clusters with landfall south of 25^∘ S are associated with only a little rainfall (cluster 1, 4, 6, 7, 8 in Fig. 9). The most crucial IWVF clusters for winter rainfall in WRF_hist are cluster 9 and in particular cluster 3, which represent northwesterly MCBs. This is in agreement with the findings of B2021, suggesting that the moistest MCBs in the recent climate have their origin in the northwest (see Fig. 7 in B2021). However, compared to the MCBs of WRF_mP (Fig. 8), the MCBs in WRF_hist (in Fig. 9) have much lower IWVF and can have different transport characteristics, e.g. indicated by wind maxima that are shifted towards the south and east. Take, for instance, the lower IWVF and different wind pattern in cluster 4 (6) for WRF_mP in Fig. 8 against cluster 9 (3) for WRF_hist in Fig. 9. The strengthening of the northwesterly MCBs in WRF_mP might suggest that the wind speeds associated with these MCBs are higher in the mid-Pliocene leading to rainfall on land, but this is not visible in the output. We computed the composite of mean wind speeds for each of the MCB clusters associated with rain in the desert. The resulting composite of wind speeds in the clusters is marked for a typical moisture transport height of 4000 m a.s.l. in Figs. 8 and 9. Interestingly, the mean wind speeds in the two clusters in WRF_mP have a similar magnitude like in WRF_hist. However, the region with peak winds is shifted to the northwest in WRF_mP. We therefore conclude that the MCBs during the mid-Pliocene that lead to rain in the Atacama Desert are not only stronger but also in some cases have origins and characteristics that are different from those that occur under present-day conditions.

Figure 8 Final mid-Pliocene wintertime IWVF clusters (yellow–green shading, in $\mathrm{kg}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$) as obtained by the combined SOM and k-means clustering for WRF_mP (outer model domain with 50 km resolution). The numbers in the red boxes are the frequencies of occurrence of the individual cluster (in days per Southern Hemisphere winter). Red–blue shading shows the fraction (in %) of the rainfall associated with the individual clusters to the total winter rainfall (inner model domain with 10 km resolution). The black contours in cluster 4 and 6 show the composite mean wind speeds at 4000 m a.s.l. (15, 17 and 19 m s⁻¹). Note that we show IWVF exceeding 350 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$ in accordance with the definition of MCBs for the mid-Pliocene (see “Data and methods”).

Figure 9 As Fig. 8 but for WRF_hist. The black contours in cluster 3 and 9 show the composite mean wind speed at 4000 m a.s.l. (15, 17 and 19 m s⁻¹).

Despite their rare occurrence, MCB cluster 4 and 6 of WRF_mP (Fig. 8) are associated with a relatively large amount of rainfall. In most parts of the Atacama Desert more than 60 % of the total winter rainfall in WRF_mP is associated with only these two MCB clusters, and for some regions north of 22^∘ S the fraction is even close to 100 % (Fig. 10a). Further, in the hyper-arid core south of 20^∘ S the mid-Pliocene winter rainfall produced by these two clusters is more than 60 % of the total annual rainfall (Fig. 10b), and it even clearly exceeds the total annual rainfall in WRF_hist (Fig. 10c). This indicates that MCBs originating in the tropical eastern Pacific are indeed a key driver for the more humid conditions in the Atacama Desert during the mid-Pliocene.

Figure 10 Fraction (in %) of the rainfall associated with MCB cluster 4 and 6 of WRF_mP to (a) total winter rainfall in WRF_mP, (b) total annual rainfall in WRF_mP and (c) total annual rainfall in WRF_hist (inner model domain with 10 km resolution).

The strength of MCBs is determined by both the winds and the atmospheric moisture content. The former is controlled by pressure gradients. We see that MCBs in the mid-Pliocene appear at the foreside of mid-tropospheric troughs over the subtropical southeastern Pacific. This is also the case for the regional MCBs for present-day conditions (B2021). Interestingly, the frequency of occurrence of troughs in the region is less during the mid-Pliocene compared to present day (Fig. 11a). Thermodynamic and dynamical changes in PlioMIP2 models lead to wetter tropics, particularly in the tropical eastern Pacific (Han et al., 2021). Furthermore, the subtropical anticyclones over the southern Pacific intensified and expanded further to the west in PlioMIP2 (Pontes et al., 2020). The change in the synoptic-scale circulation is a potential reason for the lower frequency of troughs reaching the subtropical Pacific offshore the Atacama Desert (Fig. 11a). Fewer troughs in the warmer past are consistent with findings for future global warming; e.g. Priestley and Catto (2022) found that the Southern Hemisphere extratropical storm track density decreases under future warming based on CMIP6. However, the troughs that occur extraordinarily close to the Equator (at or north of 20^∘ S) are stronger in the mid-Pliocene, measured by the mean pressure gradients (Fig. 11b), which implies an intensification and northward shift of the flow at their northern flank. This finding is consistent with wind composites we show for the northwesterly MCB clusters (Figs. 8 and 9).

Figure 11 CESM2 mid-Pliocene minus CESM2 historical for (a) number of trough points (in number per decade) and (b) mean trough point strength (in gpm).

The atmospheric moisture content is influenced by the SST. In a sensitivity study, Bozkurt et al. (2016) show that a hypothetical reduction in the SST in the tropical eastern Pacific during the March 2015 rainfall event, which occurred during anomalously high SST, significantly decreases the precipitable water along the Peruvian and Chilean coast and the rainfall amount in the Atacama Desert. We therefore assess the role of the SST in the development of MCBs originating over the tropical eastern Pacific. To this end, we computed SST composites for all days within MCB cluster 4 and 6 of WRF_mP. For those days, the SST composite is up to 1.2 ^∘C warmer than the median of the winter SST for the entire mid-Pliocene period (Fig. 12a). Over the equatorial and tropical eastern Pacific, the SST composite even exceeds the 70th winter SST percentile (Fig. 12b). As shown in Fig. 2b, this corresponds to the region with the strongest SST increase in the mid-Pliocene. Hence, this mid-Pliocene SST increase in combination with the stronger upper-level troughs (Fig. 11b) is a plausible scenario to explain why exceptionally strong MCBs, as identified in cluster 4 and 6 in WRF_mP, occur in the mid-Pliocene but are not seen during present-day conditions (WRF_hist). Our results support the finding that higher SSTs lead to stronger rainfall in the Atacama Desert, broadly consistent with the March 2015 case studied by Bozkurt et al. (2016).

Figure 12 SST composite mean over all days with MCB cluster 4 and 6 minus the (a) 50th percentile of winter SST and (b) 70th percentile of winter SST in the mid-Pliocene experiment of CESM2.

These results indicate that stronger troughs and warmer SSTs during the mid-Pliocene favour the formation of strong MCBs. It is not clear, however, how MCBs in the Atacama region develop in detail. For the southwestern United States, Knippertz and Martin (2007) attributed the formation of an MCB to a baroclinic zone and a strong subtropical jet at the foreside of a cutoff, as well as to strong convection in the Intertropical Convergence Zone (ITCZ). Although plausible, this explanation may not necessarily be applicable for MCBs in the southeastern Pacific. For the Atacama Desert under present-day conditions, B2021 traced the air in the strongest MCBs back to the Amazon basin. Their backward trajectories hardly indicate any moisture enrichment when passing the Pacific (see Fig. 4d of B2021), and consequently they conclude that the associated moisture most likely originates from the continental Amazon basin. Interestingly, the isotopic analysis of the collected rainwater in coastal regions from the March 2015 event indicates a tropical Pacific origin for the water (Jordan et al., 2019). While B2021 did not specifically address this particular event, they show tropical and subtropical Pacific origins for cases that are weaker in terms of humidity along the identified trajectories. The uncertainty in identifying the most representative target location and time for the trajectory calculation may cause stronger events to appear weaker.

4 Conclusion and outlook

We analysed global PMIP4–CMIP6 model simulations and performed regional kilometre-scale experiments to assess drivers for the increased rainfall amounts in the Atacama Desert during the warm mid-Pliocene as reported from geological evidence. Our model results suggest that the enhanced rainfall during the mid-Pliocene can be primarily explained by increased winter rainfall. Focussing on this season, we find evidence that the main drivers of the increased rainfall are moisture conveyor belts (MCBs) advecting moisture from the tropical eastern Pacific to the Atacama Desert. This conclusion is based on the following findings:

• i.

  More humid conditions in the mid-Pliocene experiment are caused by stronger rainfall events that can be twice as strong compared to the present day. The frequency of occurrence of rainfall events remains nearly unchanged.

• ii.

  The top mid-Pliocene rainfall events are often associated with strong elongated upper-level moisture fluxes. These share the typical characteristics of MCBs with an origin in the tropical eastern Pacific advecting moist air to the Atacama Desert.

• iii.

  A clustering reveals that MCBs in the mid-Pliocene are stronger and show distinct differences in the exact origin and pathways compared to those for the present day. The two mid-Pliocene clusters for MCBs suggest that the tropical eastern Pacific is the moisture source, in contrast to present-day simulations (historical experiment) for which these northwesterly MCBs are much weaker and their origin is further to the southwest.

• iv.

  The rainfall associated with the northwesterly MCBs can be up to 100 % of the regional mid-Pliocene winter rainfall and more than 60 % of the annual rainfall in the mid-Pliocene. Further, these rain amounts exceed the total annual rainfall in the present-day climate in some regions south of 20^∘ S.

We find that only the PMIP4–CMIP6 model CESM2 simulates increased rainfall in the Atacama region for the mid-Pliocene experiment, while the other models (IPSL-CM6A-LR, GISS-E2.1-G) simulate decreased rainfall. Simultaneously, CESM2 is the only model that realistically captures the sea-surface warming in the tropical and eastern Pacific. In the other GCMs the temperature increase in the mid-Pliocene is underestimated when compared to proxy data. However, it remains unclear whether the discrepancies in the simulated rainfall change patterns are caused by differences in the sea-surface temperature (SST). There is no clear improvement in tropical rainfall in tandem with reduced biases in SSTs in CMIP models (Fiedler et al., 2020). We cannot rule out that the representation of atmospheric circulation and clouds (Bony et al., 2015) contributes to differences in rainfall changes in the CMIP6 models.

Our regional evaluation is interesting in the context of the relatively high climate sensitivity of CESM2 (Gettelman et al., 2019), which might be seen as an outlier in a larger ensemble of CMIP6 simulations for other time periods (Burls and Sagoo, 2022). It was proposed to use paleosimulations as a test bed for climate model performance to constrain climate sensitivity (Burls and Sagoo, 2022; Zhu et al., 2022). Our results suggest that paleosimulations paired with regional downscaling to the kilometre scale might also be useful for better understanding and predicting regional climate changes with global warming, e.g. for the hydrological cycle that remains an outstanding challenge for global models with parameterized convection. If our mid-Pliocene simulation is a useful out-of-sample test, the fact that CESM2 outperforms other models with lower climate sensitivity for the mid-Pliocene climate in the region of the Atacama Desert would support high climate sensitivity. It would be valuable to have data from more global model simulations for the Pliocene or other warm climates for similar downscaling experiments in future research, especially from CMIP6 models with high climate sensitivity. This endeavour also requires further development of proxy data for paleoclimates, of which there are still a limited number for the Pliocene.

To the best of our knowledge this is the first time that a PMIP4 mid-Pliocene experiment has been downscaled with a regional climate model. We think that the unique datasets give new insights into the processes that are crucial for the mid-Pliocene climate in the Atacama Desert as an analogy for a future warmer world and are valuable for the interpretation of paleoclimate proxies in the region, e.g. for understanding physical processes leading to different MCBs. One interesting aspect for future research could stem from a detailed analysis of processes involved in some of the top mid-Pliocene rainfall events to gain a better understanding of the development of MCBs in the 3 ^∘C warmer world of the mid-Pliocene and to put it into the context of the El Niño–Southern Oscillation (ENSO). Solman and Menendez (2001) state that El Niño events result in an equatorward shift of the winter storm track over the subtropical Pacific. Our results suggest that anomalously high SSTs might be beneficial for the formation of MCBs from the tropical Pacific. This implies that a Pliocene mean state more similar to today's El Niño conditions, like that seen in some studies (e.g. Wara et al., 2005), may indeed contribute to wetter conditions in a warmer world via the processes identified in our study. However, the El Niño-like state of the Pliocene is subject to discussion (e.g. Fedorov et al., 2006; Watanabe et al., 2011; Oldeman et al., 2021), and hence the role of ENSO is not yet fully understood. Atmospheric rivers are not unlike MCBs and are more strongly changing in the northern Pacific than in the southern Pacific between the mid-Pliocene and pre-industrial times (Menemenlis et al., 2021), which suggests that a dynamical assessment with high spatial resolutions for the mid-Pliocene would also be interesting in other regions.

Furthermore, our findings have implications for rainfall extremes in a potential future 3 ^∘C warmer world, since the mid-Pliocene experiment is regarded as a potential analogue to what might be expected for the future. As such, our work suggests that regional rainfall extremes exceeding historic observations might occur. Such events involve atmospheric dynamics that are similar but not identical to what we know from the recent past.