A dipole pattern of orbital-scale precipitation oxygen isotope variation in North African monsoon region and the driving mechanism

On orbital timescales, the North African (NAF) monsoon variability is featured by dramatic fluctuations between wet and dry periods, which have played a significant role in early human migration and the development of agricultural civilizations. However, the spatial patterns of hydroclimate response, particularly changes in rainfall and precipitation oxygen isotopes (δ¹⁸O_p) remain poorly constrained due to the scarcity of proxy records. Here, we use the isotope-enabled Community Earth System Model (iCESM) to investigate the spatial-temporal variations of both rainfall and δ¹⁸O_p across the NAF region (15° W–35° E, 8° N–25° N) on orbital timescales. Our analysis shows that both δ¹⁸O_p and rainfall exhibit a clear precessional signal. Enhanced Northern Hemisphere summer insolation (NHSI) intensifies land-sea thermal contrast, thereby strengthening monsoon circulation and leading to widespread increases in monsoon rainfall. In contrast, δ¹⁸O_p presents a spatially dipole pattern, with depletion in the southern NAF (15° W–35° E, 8° N–17° N) and enrichment in the northern part (15° W–35° E, 17° N–25° N). Tagging experiments further reveal that the depletion in the south is primarily driven by en route depletion resulting from rainout process along air mass trajectories over the African continent, whereas enrichment in the north results from shifts in moisture sources, with a reduced contribution from distant sources and an increased influence of local sources. This work advances our understanding of past hydroclimate variability in the NAF region, while also highlighting the challenges associated with reconstructing past variations in the distribution of δ¹⁸O_p based on discrete site data.

1 Introduction

The North African (NAF) monsoon region (15° W–35° E, 8° N–25° N) is one of the most vulnerable hotspots under climate change (Lézine et al., 2011). Monsoon variability in this region exerts substantial environmental and socioeconomic effects, particularly across Sahelian countries where economies rely heavily on rainfed agriculture (Sultan et al., 2005; Sultan and Gaetani, 2016). On orbital timescales, the NAF region is highly sensitive to orbital forcing and has experienced dramatic fluctuations between wet and dry periods (deMenocal, 2004; deMenocal and Tierney, 2012; Ehrmann et al., 2017; Kutzbach et al., 2020). These changes not only drive the early human migration to the Nile Basin and the rise of agricultural civilization (McDougall et al., 2005; Vaks et al., 2007; Kutzbach et al., 2020), but also have profound impacts on Asian monsoon, ENSO activity, and even the polar ice sheets through atmospheric teleconnections (Muschitiello et al., 2015; Pausata et al., 2020). Given the important role of NAF monsoon in the global climate change, a clear understanding of its past evolution is critical for improving predictions of its future behavior.

Numerous modeling studies have been conducted extensively to investigate the NAF monsoon variability on orbital timescales. For example, Kutzbach and Liu (1997) conducted Mid-Holocene (6 ka BP) simulations with a general circulation model that asynchronously couples the atmosphere and the ocean. Their simulations show that increased insolation induces an intensified land-sea thermal contrast and enhanced summer monsoon rainfall. Similarly, Tuenter et al. (2003) used an Earth system model of intermediate complexity (ECBilt) to analyze orbital-scale signals over Africa. They demonstrated that a stronger and more northward monsoon over NAF occurs at times of minimum precession or maximum obliquity, when boreal summer insolation is high. This rainfall enhancement was attributed to a stronger land-sea thermal contrast, intensified low-level monsoonal circulation, and increased moisture transport. Similar results were also found in Bosmans et al. (2015) using the EC-Earth model. To further explore the evolution of the NAF monsoon, Kutzbach et al. (2020) performed transient simulations covering the past 140 000 years using a dynamic atmosphere-ocean general circulation model with higher resolution. Their simulations showed that orbital precession causes high seasonality in Northern Hemisphere (NH) insolation, with stronger and northward extended summer monsoon rainfall and increased winter rains in the Mediterranean Basin. These combined effects increased vegetation and narrowed the width of the Saharan-Arabian desert. Overall, the NAF monsoon rainfall has been extensively studied.

In contrast, the precipitation oxygen isotope (δ¹⁸O_p) has received comparatively less attention, largely because earlier models lacked coupled water isotope modules. Given that most proxy reconstructions are made using measurements of isotopic ratios in natural archives such as ice cores, speleothems and corals (Wang et al., 2001; Bar-Matthews et al., 2003; Cruz et al., 2005; Cheng et al., 2013), a more comprehensive understanding of the controls on isotope ratios is critical in paleoclimate research. In the context of North Africa, past hydrological variations have been inferred from δ¹⁸O signal of speleothem calcite (δ¹⁸O_c) from caves in Israel (Bar-Matthews et al., 2003). The δ¹⁸O_c in Israel serves as an indicator of past NAF monsoon variability, as wet or dry conditions in the NAF region are propagated to the Mediterranean Sea through Nile River and subsequently influence the isotopic composition of cave carbonates in Israel. This δ¹⁸O_c shows a strong precessional signal, with the enhanced Northern Hemisphere summer insolation (NHSI) corresponding to a depleted δ¹⁸O value, indicating an enhanced NAF rainfall (Bar-Matthews et al., 2003; Brahim et al., 2023). This is consistent with other moisture-related proxies from lakes (Kohfeld and Harrison, 2000; Armitage et al., 2015), terrigenous sediment (deMenocal et al., 2000; Revel et al., 2010; Ehrmann et al., 2016; Skonieczny et al., 2019; Blanchet et al., 2021), and pollen (Lézine et al., 2005; Fersi et al., 2016), showing that a higher NHSI corresponds to stronger NAF rainfall. Therefore, a major feature of the isotopic composition of precipitation in NAF region is the anticorrelation between the amount of precipitation and the proportion of heavier isotopes in the precipitation, called the “amount effect” (Dansgaard, 1964).

However, it is still unclear whether the NAF δ¹⁸O_c signal responds primarily to changes in local rainfall or is controlled by other factors. This question arises because a growing number of studies have shown that many factors other than local rainfall amount can affect δ¹⁸O_p, including changes in the relative contribution of moisture sources (Cole et al., 1999; Vuille et al., 2003), isotope composition during convective activity such as sub-cloud evaporation and diffusive exchanges between raindrops and the surrounding vapor (Lee and Fung, 2007; Risi et al., 2008; Kurita, 2013; Moore et al., 2014), and upstream rainfall effects (Lee et al., 2009; Pausata et al., 2011; Liu et al., 2014; Shi et al., 2025). These factors may contribute to isotope variations, thereby complicating the interpretation of isotope data in paleoclimate.

Over the past few decades, isotope-enabled models have evolved as valuable and well-established tools for improving our understanding of the relationship between water isotopes and climate variables. However, there are still controversies regarding the spatial patterns of δ¹⁸O_p and its forcing mechanisms over the NAF region. Herold and Lohmann (2009), using an isotope-enabled General Circulation Model (ECHAM4), identified a dipole pattern in δ¹⁸O_p across NAF, characterized by enrichment in the west and depletion in the east. They attributed western enrichment to reduced upstream rainfall depletion from Atlantic, and eastern depletion to increased local rainfall and enhanced moisture transport from the Atlantic. A similar dipole structure was later reported by Battisti et al. (2014) using the same model. Nonetheless, this west-east dipole structure has been challenged by subsequent studies showing a spatially coherent depletion pattern across the NAF region. For example, Cauquoin et al. (2019) have demonstrated a clear “amount effect” across the NAF monsoon region, with significant δ¹⁸O_p depletion coinciding with increased monsoon rainfall. This finding is further supported by Shi et al. (2023). More recently, Shi et al. (2025) examined both spatial and temporal variations in rainfall and δ¹⁸O_p based on 24 time-slice experiments spanning a full precession cycle. Their results reveal widespread δ¹⁸O_p depletion throughout North Africa, but also identify a positive temporal δ¹⁸O_p-rainfall relationship in certain inland areas, contrast to the classical “amount effect”. In these regions, the authors propose that δ¹⁸O_p variability is influenced mainly by rainfall changes in upstream areas rather than local rainfall. Therefore, despite growing efforts to simulate δ¹⁸O_p, spatial distribution and interpretations remain inconsistent across different modeling studies.

In this study, we use the isotope-enabled fully coupled Community Earth System Model (iCESM; Brady et al., 2019) developed by the National Center for Atmospheric Research (NCAR) to perform transient simulations for the past 150 000 years. In addition, we perform moisture tagging experiments to track moisture sources for precipitation and to investigate the corresponding δ¹⁸O_p response. We aim to answer the following central questions: How do monsoon rainfall and δ¹⁸O_p respond to insolation on orbital timescales in NAF region? Do they exhibit spatially coherent patterns? And what is the climatic significance of δ¹⁸O_p in NAF region? The paper is arranged as follows. Section 2 describes the Data and methods; Sect. 3 presents the Data-model comparison; Sects. 4 and 5 discuss the mechanisms underlying rainfall and δ¹⁸O_p changes, respectively; and Sect. 6 provides the Conclusions.

2 Data and methods

2.1 The 150 000-year simulation

In this study, the isotope-enabled Community Earth System Model (iCESM) is employed, which has active atmosphere, land, ocean, river transport, and sea ice component models linked through a coupler (Brady et al., 2019). We use the version 1.3 of iCESM, with a resolution of f19_g16. The atmosphere component model is the CAM5.3, which has a horizontal resolution of 1.9° × 2.5° (latitude × longitude) and 30 hybrid vertical levels (Neale et al., 2010). The land component is the CLM4, which shares the same horizontal grid as the atmosphere model (Oleson et al., 2010). The ocean and sea ice components are the POP2 and CICE4. POP2 and CICE4 here use a common grid of 320 × 384, a displaced-pole grid with poles in Greenland and Antarctica, and a nominal 1° resolution with enhancement near the equator and in the North Atlantic. The POP2 has 60 vertical layers (Smith et al., 2010; Hunke, 2010). Water isotope ratios, and the associated fluxes and isotopic fractionations, are tracked in all of the components of the hydrologic cycle: atmospheric water vapor and clouds, soil moisture and other land surface water pools, oceans, and sea ice. The iCESM has been successfully applied in multi-scale climate studies to investigate mechanisms linking the hydrological cycle and isotopic signals (Tabor et al., 2018; Hu et al., 2019; He et al., 2021; Bao et al., 2023; Wen et al., 2024). It has been shown to adequately capture key features of δ¹⁸O_p, albeit with a weaker amplitude compared to observations (Brady et al., 2019; Nusbaumer et al., 2017; Wen et al., 2024).

We conduct a 150 000-year transient simulation that is driven solely by variations in Earth's orbital parameters (i.e., precession, obliquity, and eccentricity; Berger, 1978). All other boundary conditions, such as greenhouse gas concentrations, ice sheet extent, and vegetation distribution, are held constant at pre-industrial levels (Wen et al., 2024). Although the ignorance of vegetation and dust feedbacks over the NAF region tends to suppress the magnitude of changes in rainfall and δ¹⁸O_p values (Waldmann et al., 2010; Pausata et al., 2016; Tierney et al., 2017a; Messori et al., 2019; Tabor et al., 2020), the orbital-forcing experiment still captures the dominant precessional signal and the overall phase of regional climate change (Pokras and Mix, 1987; Patricola and Cook, 2007; Weber and Tuenter, 2011; Roe et al., 2016; Cheng et al., 2020). This is sufficient for the purpose of this study.

The experiment starts 170 000 years ago with the last 150 000 years used for analysis. Due to the limitation of computational resources, the experiment is accelerated by 100-time to improve efficiency. That is, at the end of each year's simulation, we advance the orbital parameters by 100 years, so the actual length of our simulation is 1500 years. This acceleration scheme is suitable for the monsoon analysis, given that the response time of the atmosphere-upper ocean system is much faster than the shortest orbital cycle (Kutzbach et al., 2008; Wen et al., 2024). The previous analysis shows that the acceleration method leads to a delayed response only in the deep-ocean temperature, but has little impact on the near-surface quantities like temperature and precipitation (Lorenz and Lohmann, 2004; Timm and Timmermann, 2007; Yin and Berger, 2015). Since this study focuses on surface ocean-atmospheric conditions, the 100-factor acceleration scheme would not affect the conclusions regarding the NAF region.

2.2 Calculation of water isotopes from iCESM simulation

Since the observed speleothem isotope records reflect the combined influence of moisture from all seasons, the annual δ¹⁸O_p are calculated using the precipitation-weighted δ¹⁸O from 12 months following the previous works (Tabor et al., 2018; He et al., 2021):

$$\begin{array}{}\text{(1)}& {\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{p}}={\sum }_{m=\mathrm{1}}^{\mathrm{12}}{\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{m}}\cdot {\displaystyle \frac{P_{\mathrm{m}}}{P}},\end{array}$$

where m denotes the calendar month, P_m and δ¹⁸O_m represent the precipitation and δ¹⁸O in precipitation from the mth month output by the model, respectively. P is the annual total precipitation.

In order to directly compare our simulations with the observed calcite speleothem records of δ¹⁸O_c, we need to get simulated δ¹⁸O_c from the simulated δ¹⁸O_p. First, we need to convert the simulated δ¹⁸O_p from the V-SMOW scale (Vienna-Standard Mean Ocean Water; δ¹⁸O_p−SMOW) to the PDB scale (Pee Dee Belemnite; δ¹⁸O_p−PDB) following Coplen et al. (1983):

$$\begin{array}{}\text{(2)}& {\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{p}-\mathrm{PDB}}=\mathrm{0.97002}\times {\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{p}-\mathrm{SMOW}}-\mathrm{29.98}.\end{array}$$

Then, we get the simulated δ¹⁸O_c following O'Neil et al. (1969):

$$\begin{array}{}\text{(3)}& {\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{c}}={\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{p}-\mathrm{PDB}}+\mathrm{2.70}\times {\mathrm{10}}^{\mathrm{6}}/T^{\mathrm{2}}-\mathrm{3.29},\end{array}$$

where T represents the absolute cave temperature. Since this temperature approximately equals the annual mean surface air temperature at the cave site (Fairchild et al., 2012), the model annual mean surface air temperature from the corresponding model grid is used for calculation.

2.3 Tagging experiments

In order to quantify the changes in moisture sources of δ¹⁸O_p in the NAF monsoon region between high and low NHSI periods, two water tagging experiments are carried out at the times of 127 ka (high NHSI periods) and 116 ka (low NHSI periods) using the atmospheric model iCAM5.3. These two time slices correspond to the period of greatest changes in solar insolation across the past 150 000 years. The external forcing in both tagging experiments is based on the orbital parameters specific to their respective periods. The boundary conditions and sea surface isotope ratios are derived from 1000-year climatological mean states corresponding to each period, from the coupled iCESM simulation. The experiments can track water vapor evaporating from the source regions until they follow the hydrological processes in the model to the region where they rain out (sink region). Following Wen et al. (2024), the global source regions have been divided into 25 subregions, with 13 covering the ocean and the rest covering the land (see Fig. A1 in the Appendix). Each tagging experiment is run for 40 years, with the last 20 years used for analysis.

Here, the original 25 source regions are grouped into five broader regions for simplicity (Fig. A1): the African continent (AFR), the South Atlantic Ocean (EQA+SSA), the Indian Ocean (EQI+SSI), the North Atlantic Ocean (NNA+SNA), and the rest of the globe. This regrouping is based on their proximity to the NAF region and their relevance as primary contributors to precipitation patterns over NAF region. For the mean climate state, the African continent, the South Atlantic, the Indian Ocean, and the North Atlantic contribute approximately 39.5 %, 22.5 %, 16 %, and 15 % to NAF rainfall, respectively (Table A1). Collectively, these four regions account for 93 % of the total NAF rainfall. Thus, the rainfall and δ¹⁸O_p at a grid point are the sums from all five source regions:

$$\begin{array}{}\text{(4)}& P={\sum }_{i=\mathrm{1}}^{\mathrm{5}}{P}_i,{\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{p}}={\sum }_{i=\mathrm{1}}^{\mathrm{5}}{\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_i\cdot {\displaystyle \frac{P_i}{P}},\end{array}$$

where P_i and δ¹⁸O_i represent the rainfall and δ¹⁸O_p from the ith tagging source region, respectively.

2.4 Decomposition of δ¹⁸O_p

The response of δ¹⁸O_p between high and low NHSI periods, Δδ¹⁸O_p=Δδ¹⁸O_p_high − Δδ¹⁸O_p_low, can be decomposed in two steps. In the first step, Δδ¹⁸O_p is decomposed into two terms: the response induced by the changes in precipitation weight $\mathrm{\Delta }\left(\frac{P_i}{P}\right)$ and the response induced by the changes in isotope ratio Δδ¹⁸O_pi:

$$\begin{array}{}\text{(5)}& \begin{array}{rl}\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{p}}& ={\sum }_{i=\mathrm{1}}^{\mathrm{5}}{\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{p}i}\cdot \mathrm{\Delta }\left({\displaystyle \frac{P_i}{P}}\right)\\ & +{\sum }_{i=\mathrm{1}}^{\mathrm{5}}{\displaystyle \frac{P_i}{P}}\cdot \mathrm{\Delta }{\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{p}i}.\end{array}\end{array}$$

The term in left-hand side is plotted in Fig. 7b, and two terms in right-hand side of Eq. (5) are plotted in Fig. 8a and d.

In the second step, the Δδ¹⁸O_pi can be further decomposed into three terms following Tabor et al. (2018) in Eq. (6). The first term denotes the changes in vapor δ¹⁸O_i in the source region, which is produced by the net effect of local rainfall and evaporation at the source region. The second term denotes the changes in en route depletion of vapor δ¹⁸O_i due to rainout along its trajectory from the source region to the sink region. The third term denotes the changes in local condensation enrichment due to the transition from vapor to rainfall in the sink region, along with any post-condensation processes such as rain evaporation. These terms are written as follows:

$$\begin{array}{}\text{(6)}& \begin{array}{rl}\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{p}i}=& \mathrm{\Delta }{\left({\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{v},\mathrm{source}}\right)}_i\\ & +\mathrm{\Delta }{\left({\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{v},\mathrm{sink}}-{\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{v},\mathrm{source}}\right)}_i\\ & +\mathrm{\Delta }{\left({\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{p},\mathrm{sink}}-{\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{v},\mathrm{sink}}\right)}_i.\end{array}\end{array}$$

The three processes are displayed in Fig. 8e–g.

2.5 Moist static energy budget

The column-integrated moist static energy (MSE) budget is used to accurately calculate the energy changes in the atmosphere column and to better understand the mechanisms of rainfall changes in the monsoon region. The MSE balance equation is referenced from Hill et al. (2017) and Wen et al. (2022).

$$\begin{array}{}\text{(7)}& \left\{\stackrel{\mathrm{‾}}{\mathit{\omega }}{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{h}}{\partial p}}\right\}\approx {\stackrel{\mathrm{‾}}{F}}_{\mathrm{net}}-\left\{\stackrel{\mathrm{‾}}{v}\cdot {\mathrm{\nabla }}_{\mathrm{p}}\stackrel{\mathrm{‾}}{h}\right\}-\left\{{\mathrm{\nabla }}_{\mathrm{p}}\cdot \stackrel{\mathrm{‾}}{\left(v{}^{\prime }h{}^{\prime }\right)}\right\}-\left\{{\partial }_t\stackrel{\mathrm{‾}}{h}\right\},\end{array}$$

where $\left\{\cdot \right\}={\int }_{\mathrm{0}}^{P_{\mathrm{s}}}\cdot \left(\frac{\mathrm{d}p}{g}\right)$ denotes the column integration from the surface P_s to the top of the atmosphere. h represents the MSE: $h=c_{\mathrm{p}}T+gz+L_{v}q$. F_net is the net energy flux into the atmosphere column, which is calculated as the net radiation flux at the top of atmosphere (TOA) minus the net radiation flux at the surface. The net radiation flux at TOA = net shortwave − net longwave. The net radiation flux at the surface = net shortwave − net longwave − sensible heat − latent heat. Overbars are the monthly mean and primes are the transient eddy. Based on this equation, the vertical advection of MSE ($\left\{\stackrel{\mathrm{‾}}{\mathit{\omega }}{\partial }_{\mathrm{p}}\stackrel{\mathrm{‾}}{h}\right\}$) can be inferred from the sum of the net energy flux into the atmosphere column (F_net), horizontal advection ($\left\{\stackrel{\mathrm{‾}}{v}\cdot {\mathrm{\nabla }}_{\mathrm{p}}\stackrel{\mathrm{‾}}{h}\right\})$, and transient eddy activity ($\left\{{\mathrm{\nabla }}_{\mathrm{p}}\cdot \stackrel{\mathrm{‾}}{\left(v{}^{\prime }h{}^{\prime }\right)}\right\}$). Note that the $\left\{{\partial }_t\stackrel{\mathrm{‾}}{h}\right\}$ is very small in the equilibrium state, so we neglect it in the following analysis. The transient eddy term is calculated as a residual from the atmospheric energy budget. The vertical integral of MSE stratification in the troposphere (${\partial }_{\mathrm{p}}\stackrel{\mathrm{‾}}{h}$) is mostly negative in pressure coordinate. Therefore, the inferred positive vertical MSE advection ($\left\{\stackrel{\mathrm{‾}}{\mathit{\omega }}{\partial }_{\mathrm{p}}\stackrel{\mathrm{‾}}{h}\right\}>\mathrm{0}$) corresponds to ascending vertical motion ($\stackrel{\mathrm{‾}}{\mathit{\omega }}<\mathrm{0}$) and thus the intense precipitation, and vice versa (Chen and Bordoni, 2014).

3 Data-model comparison

The iCESM effectively captures the key features of modern NAF climate and δ¹⁸O_p (Fig. 1). Specifically, a low-level southwesterly flow transports humid air from the equatorial Atlantic to the African continent and converges with dry air from the Sahara, leading to heavy rainfall in the Sahel and its southern areas (Fig. 1a) (Sultan and Janicot, 2003; Rose et al., 2016; Selami et al., 2021; Datti et al., 2025). These features generally align with modern observations (Fig. 1b) and other model simulations (Cook and Vizy, 2019; Shi et al., 2024). The simulated NAF δ¹⁸O_p shows a more enriched values in the Sahel regions due to their proximity to the moisture source region over the tropical Atlantic Ocean, and more depleted values farther inland, as the air masses undergo repeated condensation and rainfall cycles during their movement toward the continental interior (shading in Fig. 1c) (Shi et al., 2023). This simulated pattern generally agrees with observations from the Global Network of Isotopes in Precipitation (GNIP) (circles in Fig. 1c) and other climate models (Cauquoin et al., 2019; Shi et al., 2023, 2025), although the model exhibits more negative values compared to the GNIP data (Fig. 1d).

Figure 1 Modern climatology of the North Africa (NAF) summer monsoon and δ¹⁸O_p. (a) Simulated boreal summer (JJA) rainfall (shading; mm d⁻¹) and 900 hPa wind (vector; m s⁻¹). (b) Same as (a) but for observation during 1940–2024 from ERA5 (https://cds.climate.copernicus.eu/, last access: October 2025). (c) Simulated annual precipitation-weighted δ¹⁸O_p (shading; ‰) and observed δ¹⁸O_p from GNIP (https://nucleus.iaea.org/wiser/explore/, last access: October 2025) (circle; ‰). (d) Point-to-point scatter plot of observed GNIP δ¹⁸O_p versus model-simulated δ¹⁸O_p. The simulated data is averaged over the last 1000 years (1 ka-0) of the transient simulation. In (a)–(c), the black rectangle region (15° W–35° E, 8° N–25° N) is the study area of the NAF monsoon and is used for subsequent regional averaging. The blue curve plots the African Nile River. The red triangles mark the locations of the proxy sites in Fig. 2.

Moreover, our simulations can well reproduce the NAF monsoon evolution during the past 150 000 years (Fig. 2). In this work, the NAF monsoon region is defined as 8° N to 25° N and 15° W to 35° E, which is denoted as black box in Fig. 1 following Wen et al. (2022). The simulated δ¹⁸O_c varies in phase with NHSI, exhibiting strong precessional signal (orange line in Fig. 2b). Specifically, higher NHSI corresponds to more depleted δ¹⁸O_c values. This is in agreement with observed δ¹⁸O_c records from Soreq and Peqiin caves (black and grey lines in Fig. 2b), although the simulated amplitude is smaller than that in the reconstructions. It should be noted that although the Soreq and Peqiin caves are located outside the defined NAF monsoon region, their δ¹⁸O_c records primarily reflect hydroclimate changes over the NAF monsoon region. This is because the hydroclimate signal from the NAF monsoon region can be transmited to the Eastern Mediterranean Sea through Nile River discharge, where it is recorded in both seawater δ¹⁸O and speleothem δ¹⁸O_c (Ayalon et al., 2002; Bar-Matthews et al., 2003; Rohling et al., 2015).

Figure 2 Data-model comparison. (a) Northern Hemisphere summer insolation (NHSI) at 30° N (red; W m⁻²) and Precession parameter (grey; Berger, 1978). (b) Simulated-δ¹⁸O_c derived from simulated-δ¹⁸O_p (see Sect. 2.2 for details) in the NAF region (orange; ‰) and observed δ¹⁸O_c records from Soreq Cave (black; ‰) and Peqiin Cave (grey; ‰), Israel (Bar-Matthews et al., 2003). (c) Simulated NAF rainfall (green; mm d⁻¹) and color reflectance (black; %) at 540 nm from ODP Site 968, Mediterranean (Ziegler et al., 2010). The correlation coefficient calculated for simulated-δ¹⁸O_c and simulated rainfall is −0.82.

Additionally, our simulations show an increased rainfall when NHSI is high (green line in Fig. 2c), which is consistent with rainfall records from ODP Site 968 (organic-rich layers, i.e., 540 nm reflectance in Fig. 2d; Ziegler et al., 2010). The sapropel layers in this sediment core reflect deep-water anoxia and increased primary productivity, which is mainly related to the increased run-off inputs into the eastern Mediterranean driven by the NAF monsoon-induced continental rainfall (Ziegler et al., 2010; deMenocal and Tierney, 2012; Grant et al., 2016, 2017). The appearance of the sapropel layers corresponds to the intensification of the NAF monsoon rainfall.

Therefore, the high NHSI corresponds to regional averaged lower δ¹⁸O_c values and higher rainfall, exhibiting a strong negative correlation between the two variables ($r=-\mathrm{0.82}$; Fig. 2). This relationship is consistent with the “amount effect” in terms of the statistical relationship. Notably, a strong precession cycle is evident in summer insolation at 30° N (Berger and Pestiaux, 1984; Berger et al., 1993), as well as in both rainfall and δ¹⁸O_c. These findings point out that the δ¹⁸O_c and rainfall variations are highly sensitive to precession-driven insolation changes, emphasizing the importance of precessional forcing in shaping low-latitude climate dynamics.

We further examine the coherent spatial-temporal response of rainfall and water isotope by conducting Multivariate Empirical Orthogonal Function (MVEOF) analysis. Previous studies have shown that the δ¹⁸O_p plays a dominant role in determining the speleothem δ¹⁸O_c (Cheng et al., 2012, 2016; Liu and Battisti, 2015). Therefore, for simplicity, we subsequently use δ¹⁸O_p for further analyses. The results show that the explained variance of the first MVEOF mode accounts for 35 % of the total variance, with the time coefficient following NHSI and exhibiting a strong precessional signal (Fig. 3a).

Figure 3 The general pattern of δ¹⁸O_p and NAF summer climate change during the past 150 000 years. (a) Normalized principal component time series (black) of the first MVEOF mode (PC1; explained 35.6 % of the total variance) between summer rainfall and δ¹⁸O_p and NHSI (red; W m⁻²). The red circles represent the high and low NHSI periods used in the following composite analysis in Fig. 4 and 5. (b) Surface temperature (shading; K) and sea level pressure (contour; Pa) regressed on the normalized time coefficient of PC1 [black in (a)]. The solid (dashed) contours denote positive (negative) values. (c) As in (b) but for rainfall (shading; mm d⁻¹) and 900 hPa wind (vector; m s⁻¹). (d) As in (b) but for surface moist static energy (MSE) (shading; J m⁻²) and atmospheric precipitable water (contour; kg m⁻²). (e) As in (b) but for precipitation-weighted δ¹⁸O_p (shading; ‰). The black rectangle is the NAF monsoon region.

Spatial analysis reveals that NHSI induces a spatially homogeneous rainfall response across the NAF region (Fig. 3c). Physically, when the NHSI is high, there is a dramatic warming over the Eurasian continent, which triggers the anomalous low-pressure over the continent and enhances the land-sea thermal contrast (Fig. 3b). This strengthens southwesterly winds that transport moisture from the equatorial Atlantic into Africa, resulting in pronounced vapor convergence over the Sahara-Sahel region (Fig. 3c, d). Concurrently, increased surface MSE destabilizes the atmospheric column, promoting upward vertical motion (Patricola and Cook, 2007; Roe et al., 2016). Together, these processes enhance deep convection and significantly increase rainfall across the NAF monsoon region (Fig. 3c). Additionally, the increased cloud formation and rainfall cause a pronounced cooling between 10° N and 20° N, consistent with previous findings (Patricola and Cook, 2007; Herold and Lohmann, 2009; Marzocchi et al., 2015; Wen et al., 2022). Despite this cooling, Patricola and Cook (2007) suggest that elevated low-level moisture content dominates the MSE change, further destabilizing the atmospheric column and reinforcing convection. Thus, although surface temperatures decrease in this latitudinal band (Fig. 3b), MSE continues to rise due to increased humidity, which is sufficient to sustain deep convection (Fig. 3d). A strengthened NAF monsoon also weakens the northeasterly winds over the Sahara (Fig. 3c). This is consistent with the dust records observed in Atlantic sediment cores located offshore at the NAF continent (Skonieczny et al., 2019; O'Mara et al., 2022; Crocker et al., 2022). Furthermore, orbitally forced land warming amplifies the interhemispheric temperature gradient, shifting the ITCZ northward (Schneider et al., 2014; Fig. 4). This shift extends monsoon rainfall farther north and contracts the Saharan arid belt (Fig. 3c).

Figure 4 Changes in model rainfall and ITCZ position during boreal summer between high and low NHSI periods. (a) Changes in summer rainfall (shading; mm d⁻¹). The contours denote the ITCZ position, which is defined as the latitude of the precipitation centroid between 30° N and 30° S following Voigt et al. (2016). The red contour is for high NHSI periods, and the blue contour is for low NHSI periods. (b) Seasonal latitudinal shifts of the area-averaged ITCZ position over the longitude range of 15° W–35° E. Compared to the low NHSI periods, the ITCZ shifts northward by an average of 2.3° during the high NHSI periods.

Since the observed δ¹⁸O_c records evolves coherently with the NAF monsoon rainfall as in “amount effect” (Fig. 2b, c), one may expect a spatially uniformly distributed δ¹⁸O_p response to insolation forcing similar to rainfall response (Fig. 3c). However, the simulated δ¹⁸O_p exhibits a dipole response with increasing in the northern NAF but decreasing in the southern part (Fig. 3e). This pattern shows a positive correlation with rainfall in the north but a negative correlation in the south. This implies the complex regional response of δ¹⁸O_p to external forcing.

4 Mechanism of rainfall change

Here, we employ the MSE budget to elucidate the physical mechanisms behind monsoon rainfall responses. In the climatology of low NHSI periods, the mean state exhibits intense rainfall over the Sahel and its southern regions (Fig. 1a, b), accompanied by positive vertical MSE advection (Fig. 5a). This rainfall is initiated by strong energy input into the atmospheric column (Fig. 5b), which drives vertical MSE advection in the tropics (Fig. 5a) and, consequently, sustains deep convective rainfall in these areas (Fig. 1a). These findings align with previous studies indicating that the vertical MSE advection is primarily balanced by the net energy flux (Chen and Bordoni, 2014; Hill et al., 2017). To the north of the rainfall band, pronounced subsidence is observed, manifested by negative vertical MSE advection (Fig. 5a), consistent with the low rainfall amounts and the presence of the Sahara Desert in this region (Fig. 1a). This suppressed convection is caused by negative horizontal MSE advection (Fig. 5c), in line with anticyclonic circulation and divergent sinking flow (Fig. 1a). The transient eddy term, which reflects northward moisture transport associated with African easterly waves (Hill et al., 2017), exhibits a meridional dipolar structure (Fig. 5d) that is opposite to the rainfall distribution.

Figure 5 Composite analysis of summer MSE budget. (a–d) Vertical MSE advection term ($\left\{\stackrel{\mathrm{‾}}{\mathit{\omega }}{\partial }_{\mathrm{p}}\stackrel{\mathrm{‾}}{h}\right\}$), net radiation term (${\stackrel{\mathrm{‾}}{F}}_{\mathrm{net}}$), horizontal MSE advection term ($-\left\{\stackrel{\mathrm{‾}}{v}\cdot {}_{\mathrm{p}}\stackrel{\mathrm{‾}}{h}\right\}$), and transient eddy term ($-\left\{{\mathrm{\nabla }}_{\mathrm{p}}\cdot \stackrel{\mathrm{‾}}{\left(v{}^{\prime }h{}^{\prime }\right)}\right\})$ during low NHSI periods. (e–h) Same as (a–d), but for the difference between high NHSI and low NHSI periods. The black rectangle is the NAF monsoon region.

Comparing high with low NHSI periods, the Sahara-Sahel region experiences increased energy input into the atmospheric column (Fig. 5f). This enhanced net energy forcing primarily stems from two factors (Fig. 6): First, greater shortwave radiation absorption in the atmosphere; and second, increased surface latent heat flux emission to the atmosphere resulting from higher rainfall and evaporation. The increased net energy in the Sahara-Sahel region directly intensifies upward motion (Fig. 5e), which intensifies the NAF monsoon and leads to increased rainfall (Fig. 3c). Concurrently, strengthened southwesterly winds contribute to positive horizontal MSE advection (Fig. 5g), indicating energy transport from adjacent tropical high-MSE regions into the Sahara-Sahel. This horizontal energy flux further sustains the increased rainfall and expands the influence of the monsoon system. The transient eddy term exhibits a north-south dipole pattern (Fig. 5h), but it contributes little to the area-averaged vertical MSE advection.

Figure 6 Changes in energy entering the atmosphere during boreal summer between high and low NHSI periods. (a) Difference in downward shortwave radiation between the top of atmosphere and the surface, denote the net shortwave radiation into atmosphere. (b) Same as (a) but for the net longwave radiation change. (c) Surface sensible heat flux. (d) Surface latent heat flux.

Additionally, we employ moisture tagging experiments to quantify the contributions of different moisture sources to the NAF monsoon rainfall response (Fig. 7a). Although the mean state rainfall in the NAF region is primarily supplied by the African continent (about 46 %), the rainfall increase from low NHSI periods to high NHSI periods is additionally contributed from the Atlantic Ocean (Fig. 7a). Quantitatively, the North Atlantic Ocean, South Atlantic Ocean, and African continent contributes 18 %, 38 %, and 28 % of rainfall increase. These moisture contributions align with the dynamical framework in which intensified southwesterly monsoon winds enhance moisture transport from the Atlantic Ocean.

Figure 7 Tracking ¹⁶O and δ¹⁸O_p from source regions (high-low NHSI in tagging experiments). Summer rainfall difference (mm d⁻¹) between high and low NHSI periods originating from moisture source regions (a1) North Atlantic Ocean (NAt), (a2) South Atlantic Ocean (SAt), (a3) Indian Ocean (Ind), (a4) African continent (AFR), (a5) other regions and (a6) the global total of 25 sub regions. (b1)–(b6) Same as (a1)–(a6) but for the difference of δ¹⁸O_p between high and low NHSI periods. The black rectangles mark the NAF monsoon region.

5 Mechanism of δ¹⁸O_p change

As illustrated above, δ¹⁸O_p exhibits a dipole response with enrichment in the north but depletion in the south in response to high NHSI (Fig. 3e). Given the overall increase in rainfall across the NAF region, the change in δ¹⁸O_p indicates a positive rainfall-δ¹⁸O_p relationship in the north and negative relationship in the south. This contrast highlights the complexity of regional δ¹⁸O_p responses to external forcing and raises a critical question: Which hydrological processes are truly reflected by these δ¹⁸O_p variations? To explore the mechanisms behind this pattern, we analyze the results from the aforementioned moisture tagging experiments.

Moisture tagging experiments perfectly reproduce meridional dipole pattern in δ¹⁸O_p between high and low NHSI periods (Fig. 7b6 vs. Fig. 3e). In light of this dipole response, we divide the monsoon region into two subregions for further analysis: a northern region (15° W–35° E, 17° N–25° N) and a southern region (15° W–35° E, 8° N–17° N). The tagging experiments indicate that the dipole response of δ¹⁸O_p in the NAF region are primarily influenced by moisture originating from the Atlantic Ocean and the African continent (Fig. 7b), which is largely consistent with the moisture contribution to rainfall changes discussed earlier. Specifically, the δ¹⁸O_p enrichment in the northern region is mainly attributed to moisture from the North Atlantic, whereas depletion in the southern region is predominantly contributed by the African continent and South Atlantic.

To further understand the mechanism behind the dipole response of δ¹⁸O_p, we first separate the changes in δ¹⁸O_p into two terms: the contribution from changes in precipitation weight δ¹⁸O${}_{\mathrm{p}i}\cdot \mathrm{\Delta }\left(\frac{P_i}{P}\right)$ (Fig. 8a) and that from changes in the isotope ratio $\frac{P_i}{P}\cdot \mathrm{\Delta }{\mathit{\delta }}^{\mathrm{18}}$O_pi (Fig. 8d) based on the Eq. (5). The decomposition suggests that dipole δ¹⁸O_p responses in the northern and southern NAF regions are governed by distinct hydrological processes. The δ¹⁸O_p enrichment in the northern NAF region is almost entirely attributable to changes in precipitation weight from North Atlantic (Fig. 8a), whereas the δ¹⁸O_p depletion in the southern NAF region is mainly controlled by change in the value of isotope ratio from the African continent, as well as South Atlantic and Indian Ocean (Fig. 8d). In the following sections, we discuss the mechanisms driving δ¹⁸O_p variability in the northern and southern NAF regions separately.

Figure 8 Decomposition of δ¹⁸O_p response. (a1–a6) δ¹⁸O_p response due to changes in precipitation weight ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{p}i}\cdot \mathrm{\Delta }\left(\frac{P_i}{P}\right)$. (b1–b6) Changes in precipitation weight $\mathrm{\Delta }\left(\frac{P_i}{P}\right)$ from each source region. (c1–c6) Climatological δ¹⁸O value from each source region. (d1–d6) δ¹⁸O_p response due to changes in δ¹⁸O_p value $\frac{P_i}{P}\cdot \mathrm{\Delta }{\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{p}i}$, which are further decomposed into the change due to (e1–e6) source δ¹⁸O_v, (f1–f6) en route depletion and (g1–g6) condensation enrichment. In addition to marking the NAF monsoon region with a black rectangle, a north-south division of the NAF region is also marked in the sixth column of subplots.

In the northern NAF region, as we illustrated above, the δ¹⁸O_p enrichment is almost entirely attributable to contributions related to precipitation weight from North Atlantic (Fig. 8a). To further explore this, we plot the changes in precipitation weight $\mathrm{\Delta }\left(\frac{P_i}{P}\right)$ along with climatological mean isotopic value δ¹⁸O in Fig. 8b and c. It shows that in the northern NAF region, the δ¹⁸O_p enrichment is almost entirely attributable to changes in precipitation weight from North Atlantic. Here, we examine how changes in precipitation weight leads to δ¹⁸O_p enrichment in the region. Climatologically, rainfall in the northern NAF region is primarily supplied by two moisture sources: the remote North Atlantic (∼ 28 %) and the local African continent (∼ 41 %) (Table A2). From low to high NHSI periods, the precipitation weight from the remote North Atlantic source decreases (Fig. 8b1), whereas that from the local African source increases (Fig. 8b4). As a result, the reduction in the North Atlantic contribution (∼ −11 %) is nearly offset by the increase from the African source (∼ 13 %) (Table A2). However, the δ¹⁸O values of these sources differ markedly: moisture originating from the North Atlantic is significantly more depleted (δ¹⁸O ≈ ∼ −8 ‰) compared to that from the African continent (δ¹⁸O ≈ ∼ −2 ‰) (Fig. 8c1 vs, c4). Therefore, despite the comparable magnitude of changes in precipitation weight, the contrast in isotopic values between the two sources results in a net enrichment of δ¹⁸O_p. In other words, from low to high NHSI periods, the reduced contribution of distal, isotopically depleted moisture source and an increased contribution of nearby, isotopically enriched moisture source, ultimately driving the positive δ¹⁸O_p anomaly in the northern NAF region.

In the southern NAF region, the δ¹⁸O_p depletion is mainly controlled by change in the value of isotope ratio from the African continent. Our decomposition analysis (based on Eq. 6) shows that rainfall depletion from African continent (Fig. 8f4) determines the δ¹⁸O_p depletion in the southern NAF region. This occurs because the southwesterly winds are enhanced over the NAF regions, which migrates the ITCZ northward by 2.3° (Fig. 4) and extending rainfall eastward into the Horn of Africa (Fig. 3c). This notable expansion of the monsoon system is accompanied by the transport of air parcels away from their source locations, during which they experience rainout processes, leading to progressive isotopic depletion. In addition, the en route depletion from South Atlantic and Indian Ocean sources is also responsible for δ¹⁸O_p decrease in west and east of southern NAF region (Fig. 8f2 and f3), respectively, due to rainout process. In contrast, the other two terms associated with source effect (Fig. 8e) and local condensation (Fig. 8g) have minimal influence on δ¹⁸O_p changes in the region. Therefore, although the negative correlation between rainfall and δ¹⁸O_p observed in the southern NAF region appears consistent with the isotopic “amount effect”, our study indicates that changes in NAF δ¹⁸O_p cannot be simply attributed to local rainfall. Instead, they are significantly influenced by en route depletion processes. This finding suggests that future interpretations of oxygen isotope records in this region should be approached with caution.

6 Conclusions

In this study, we investigate the spatial-temporal evolutions of both rainfall and δ¹⁸O_p across the NAF region using transient simulations spanning the past 150 000 years. We find that the NAF hydroclimate is controlled by low-latitude insolation and exhibits a strong precessional signal, with regional average rainfall and δ¹⁸O_p negatively correlated, seemingly consistent with the sense of “amount effect”. However, rainfall and δ¹⁸O_p display distinct spatial patterns. During high NHSI periods, rainfall shows a regionally uniform increase, whereas δ¹⁸O_p exhibits a north-south dipole pattern, characterized by enrichment in the northern part and depletion in the southern part. The widespread increase in rainfall across the NAF region is driven by enhanced solar radiation, which destabilizes the atmospheric column and intensifies monsoon circulation, further sustained by horizontal advection of MSE. In contrast, the north-south dipole response of δ¹⁸O_p is governed by different hydrological processes. Tagging experiments suggest that δ¹⁸O_p enrichment in the northern NAF region results from shifts in moisture sources, with a reduced contribution from the distant Atlantic Ocean and an increased contribution from local continental sources. On the other hand, δ¹⁸O_p depletion in the southern NAF region is primarily attributed to en route rainout processes, as intensified monsoon circulation transports air parcels away from their source regions, during which they undergo progressive isotopic depletion.

The results of our simulations are broadly consistent with previous studies. The simulated rainfall expands poleward, in agreement with lake-level reconstructions (Tierney et al., 2011) and marine sediment records from the western Sahara (Tierney et al., 2017a). In addition, the rainfall extends eastward into the Horn of Africa, driven by strong southwesterly winds, which is supported by marine sediment evidence from the Gulf of Aden (Tierney et al., 2017b). The West African precipitation anomaly penetrates farther north than its East African counterpart, consistent with earlier modeling results (Schurgers et al., 2007). The simulated δ¹⁸O_p also aligns well with speleothem and seawater δ¹⁸O records from the eastern Mediterranean Sea, which reflect hydroclimate conditions over the NAF region via Nile River discharge (Ziegler et al., 2010). Furthermore, the spatial pattern of δ¹⁸O_p resembles that reported by Shi et al. (2023) with different climate models, with depleted values in the southern NAF region and enriched values in the northwestern part. However, in their simulation, the area of δ¹⁸O_p enrichment is more spatially confined.

Nonetheless, some discrepancies remain between our simulated δ¹⁸O_c and previous studies. For instance, the simulated δ¹⁸O_c exhibits a much narrower range of variability, from approximately −7 ‰ to −10 ‰ compared to the observed range of −2 ‰ to −8 ‰ . This mismatch may arise from either model biases or uncertainties in the proxy data. Relative to pollen-based reconstructions (Tierney et al., 2017a; Braconnot et al., 2012), our simulated precipitation anomalies are less pronounced (Bartlein et al., 2011; Shi et al., 2022), which likely results in weaker isotopic signals. In addition, our simulation considers only insolation forcing, whereas in reality, other factors such as lakes, vegetation, and ecosystem dynamics also play important roles. Model resolution is another key factor that can influence the results (Werner et al., 2011). On the other hand, uncertainties also exist in the proxy data themselves. For example, δ¹⁸O_p signals preserved in speleothems may reflect a complex combination of local environmental processes rather than large-scale climate signals alone (Lachniet, 2009).

Direct spatial comparison between simulated δ¹⁸O and proxy observations on orbital timescales is not currently feasible, as no direct δ¹⁸O reconstructions are available from the African continent. The north–south dipole pattern of δ¹⁸O_p simulated in this study differs from the west–east dipole pattern previously using the ECHAM4 model (Herold and Lohmann, 2009). This discrepancy likely stems from differences in model-generated climate responses. In Herold and Lohmann (2009), simulated rainfall decreases over the western Sahara, which contrasts with sediment records indicating increased rainfall in that region (Tierney et al., 2017a). This negative rainfall bias may have led to an underestimation of isotopic depletion processes, resulting in δ¹⁸O_p enrichment in their simulations. Our simulated north–south dipole pattern also differs from the widespread δ¹⁸O_p depletion over North Africa simulated by MPI-ESM-wiso (Cauquoin et al., 2019). In that study, the simulated rainfall increase does not extend to the Horn of Africa, contrary to sediment records from the Gulf of Aden (Tierney et al., 2017b). This regional limitation likely underestimates the rainfall increase and, consequently, the moisture contribution from the African continent. The resulting discrepancy in moisture sources may explain the absence of δ¹⁸O_p enrichment in the northern NAF region in their simulations compared with ours.

Although direct spatial validation of δ¹⁸O_p against orbital-scale observations remains limited due to the lack of continent-based δ¹⁸O proxies, the strong agreement between our simulated rainfall patterns and multiple proxy reconstructions lends confidence to our simulated δ¹⁸O_p distributions. Furthermore, our tagging experiments enable a quantitative assessment of the processes driving δ¹⁸O_p variations. In the southern NAF, upstream en route depletion of vapor δ¹⁸O accounts for 81 % of the total depletion, while the contribution of the “amount effect” is negligible. This quantitative approach advances beyond previous studies that relied primarily on statistical relationships between δ¹⁸O_p and rainfall to infer underlying mechanisms (Cauquoin et al., 2019). The dipole pattern of δ¹⁸O_p simulated here provides a distinct isotopic fingerprint that can be tested against future proxy evidence.

Our results reveal that the forcing mechanisms governing δ¹⁸O_p across North Africa are complex. While the “amount effect” is widely invoked to interpret δ¹⁸O_p in tropical monsoon regions, its applicability varies across different subregions of North Africa. In the southern NAF, although δ¹⁸O_p and rainfall are negatively correlated, that is consistent with the “amount effect” in a statistical sense, the primary control on δ¹⁸O_p is en route depletion processes rather than local precipitation amount. The interpretation of δ¹⁸O_p on orbital timescales also differs from that on seasonal timescales. At the seasonal scale, Risi et al. (2008) demonstrated that convective activity is the main control on the isotopic composition of precipitation in North Africa, supporting the “amount effect.” This contrast highlights the complex, scale-dependent relationship between δ¹⁸O_p and rainfall. On orbital timescales, it is essential to consider changes in rainfall along upstream vapor transport trajectories when interpreting isotopic proxies of past climate change. Although δ¹⁸O_p cannot be used to directly infer past monsoon rainfall intensity, it remains a robust indicator of monsoon strength. Overall, depleted δ¹⁸O_p values across North Africa are associated with an intensified NAF monsoon, characterized by enhanced southwesterly winds and increased monsoon rainfall.

There are still caveats in our simulation. We only consider insolation forcing. Other factors such as ice sheet dynamics and vegetation feedbacks may also influence the amplitude of the simulated results. For instance, the simulated northward extension of rainfall in the NAF region is more limited compared to reconstructions, which suggest that the African Humid Period extended as far north as ∼ 31° N (Tierney et al., 2017a). This underestimation of rainfall increase may be due to the absence of dynamic vegetation and dust emission feedbacks in our model (Pausata et al., 2016), and could also constrain the simulated range of δ¹⁸O_c. Furthermore, recent studies indicate that polar ice sheets can influence tropical monsoon systems by altering the interhemispheric pressure gradient (An et al., 2024). The lack of ice sheet forcing in our simulation may therefore contribute to an underestimation of rainfall over North Africa, potentially affecting the amplitude and variability of the simulated δ¹⁸O_c as well.

Appendix A: Additional material

Figure A1 Source region delineation in the tagging experiments and regrouping for analysis. (a) Regional delineation of the global ocean. The black dashed boxes mark the ocean sub regions used for the experiments. The blue polygons mark the regrouped ocean regions used for the analysis. (b) Similar to (a) but for regional delineation of the global land.

Table A1 Regional average of annual mean rainfall (Rainfall; mm d⁻¹), precipitation weight (Wgt; proportion of total precipitation) and annual mean δ¹⁸O in precipitation (‰) in the NAF monsoon region (15° W–35° E, 8° N–25° N).

Table A2 Same as Table A1, but for the northern NAF monsoon region (15° W-35° E, 17° N-25° N).

Table A3 Same as Table A1, but for the southern NAF monsoon region (15° W-35° E, 8° N-17° N).