Global River Discharge and Floods in the Warmer Climate of the Last Interglacial

We investigate hydrology during a past climate slightly warmer than the present: the last interglacial (LIG). With daily output of preindustrial and LIG simulations from eight new climate models we force hydrological model PCR‐GLOBWB and in turn hydrodynamic model CaMa‐Flood. Compared to preindustrial, annual mean LIG runoff, discharge, and 100‐yr flood volume are considerably larger in the Northern Hemisphere, by 14%, 25%, and 82%, respectively. Anomalies are negative in the Southern Hemisphere. In some boreal regions, LIG runoff and discharge are lower despite higher precipitation, due to the higher temperatures and evaporation. LIG discharge is much higher for the Niger, Congo, Nile, Ganges, Irrawaddy, and Pearl and lower for the Mississippi, Saint Lawrence, Amazon, Paraná, Orange, Zambesi, Danube, and Ob. Discharge is seasonally postponed in tropical rivers affected by monsoon changes. Results agree with published proxies on the sign of discharge anomaly in 15 of 23 sites where comparison is possible.


Introduction
With ongoing climate change, the need is pressing to anticipate the response of the hydrological cycle to increasing global temperatures. This can be achieved through modeling frameworks that include emission scenarios and simulation of climate, hydrology, and river hydrodynamics (Arnell & Gosling, 2016;Hirabayashi et al., 2013;van Vliet et al., 2013;Ward et al., 2017). One opportunity to increase our grasp of the Earth's hydroclimate is offered by studying past responses of hydrological systems to different climate states that occurred during the Earth's history, on which geology grants evidence (Harrison et al., 2015;Lunt, Elderfield, et al., 2013). As climate swung between glacials and interglacials of the Quaternary, global, or at least hemispheric, temperatures sometimes exceeded those of the present day (PAGES, 2016). The most recent of those intervals is the last interglacial (LIG; 129 to 116 kyr ago). The LIG holds special relevance to the future, because of the unique combination of two factors: its proximity in geological time, which means that the conditions of most climate system components are comparable to the present, and its partial temperature analogy with possible futures. Atmospheric surface temperatures were higher than the preindustrial (PI) period and probably slightly higher than the present, with largest anomalies in the Northern Hemisphere and particularly at its high latitudes (CAPE Members, 2006;Turney & Jones, 2010); however, global quantifications are still inadequate (Capron et al., 2017). Global estimates of LIG sea surface temperature anomalies compared to the present range between 0°C and +0.9°C (Hoffman et al., 2017;Turney et al., 2020) and mean ocean temperature anomaly is estimated to be approximately +1°C (Shackleton et al., 2020). The warmth of the LIG was due to the different latitudinal distribution of seasonal insolation and the resulting internal climate feedbacks, principally related to Arctic amplification of warming. The LIG has been abundantly studied in terms of temperature, using climate models (e.g., Lunt, Abe-Ouchi, et al., 2013;Otto-Bliesner et al., 2013 and with proxies covering some Northern Hemisphere continents and many ocean subbasins (Capron et al., 2017;Hoffman et al., 2017;Turney et al., 2020), and recently also in terms of precipitation (Nikolova et al., 2013;Pedersen et al., 2017). In particular, the higher boreal LIG precipitation likely resulted from high-latitude warming and associated reduction in boreal latitudinal temperature gradients, ultimately linked to differences in insolation (Scussolini et al., 2019). While reconstructions of LIG runoff and river discharge are available for some locations (supporting information Table S1), for example, for the Nile (Wu et al., 2018) and Congo rivers (Gingele et al., 1998), the hydrological implications of the different atmospheric and land climate of the LIG remain unexplored to this day. The application of hydrological modeling to study large-scale changes before the twentieth century is rare. Exceptions are the modeling of Holocene discharge for a set of large river basins (Aerts et al., 2006;Ward et al., 2007).
In this study we (1) model global anomalies between LIG and PI surface runoff, river discharge, and floods; (2) validate those anomalies by comparison with the available proxies; and (3) quantify modeled anomalies of hydrological extremes, in discharge and floods. To these ends, we use results from the latest generation of climate models to force state-of-the-art global hydrological and hydrodynamic models. Our discharge and flood results form the first global picture of hydrological conditions for the LIG.

Materials and Methods
To simulate the hydrology and hydrodynamics of the LIG and the PI, we use data from eight global climate models of the Paleoclimate Modelling Intercomparison Project Phase 4: CESM1.2, CESM2, EC-EARTH3.2, HadGEM3-GC3.1, IPSL-CM6-LR, MPI-ESM 1.2.01p1-LR, NorESM1-F, and NUIST-CSM. Depending on the model, the atmosphere's resolution ranges from 1°to 2°, and simulations last between 100 and 300 yr. Simulations are conducted at equilibrium conditions, with forcing reflecting 127 ka for the LIG and CE 1850 for the PI. From the climate models, daily total precipitation and near-surface air temperature are inputted in the global hydrological model PCR-GLOBWB Version 2 (Sutanudjaja et al., 2018). We test the climate model data against climate reanalysis (Text S2 and Figures S1 and S2) and determine that bias correction of climate variables prior to input in the hydrological model is not warranted. Input data are bilinearly interpolated to the 0.5°resolution of the hydrological simulations. The time step is daily. Potential reference evapotranspiration is calculated from air temperature using the Hamon method (Hamon, 1963). Simulations reflect natural conditions as much as possible, and human factors are excluded. Both LIG and PI simulations use present-day maps for land cover, river networks, and sea levels. In hydrological simulations, we let the groundwater storage adjust to the different climatology and then spin up the model for 50 yr. From PCR-GLOBWB we analyze surface runoff. We then interpolate the daily surface runoff from PCR-GLOBWB for input into the global hydrodynamic model CaMa-Flood (Yamazaki et al., 2011), which solves the local inertial equation and runs at 0.25°resolution. From CaMa-Flood's output we analyze routed river discharge, flood (inundated) area, and flood volume (i.e., flood area × flood depth) per grid cell. All variables are averaged across the whole model time series, unless specified. Further, we investigate changes at the river basin scale, for 27 major rivers selected to approximately cover all climatic regions ( Figure S10). We report anomalies in seasonal discharge at the river mouth in the form of average annual hydrographs, along the boreal water year (October to September). To correctly compare the two simulations, we account for the seasonal implications of a different Earth's orbit during the LIG by applying the angular definition of calendar (Bartlein & Shafer, 2019). See Text S1 for additional details on each model and simulation and Text S2 and Figures S3-S5 for a test of the sensitivity of our results with respect to bias correction of the climate data sets.
For the vast majority of regions, anomalies of LIG runoff and precipitation coincide in sign ( Figure 2). In departure from this rule, some areas present negative runoff anomaly where LIG precipitation is positive, mostly over eastern and northern Asia, northern continental Canada, northern South America, and Iberian Peninsula. In most of these cases, the LIG runoff coefficient is lower ( Figure S7), indicating that a smaller portion of precipitation converts into runoff. This is mainly a result of higher evaporation ( Figure S9) due to higher air temperature ( Figure S6). Comparing anomalies in discharge (Figures 3 and S11) and precipitation, areas where negative LIG discharge corresponds to positive precipitation are found mostly over northeast Asia (Amur and Lena river basins), northern Canada, and the east coast of North America; positive discharge corresponds to negative precipitation over smaller areas: mainly northeast Brazil, northwestern North America, and Madagascar. Last, comparing runoff and discharge, areas of positive discharge anomaly and negative runoff emerge mostly over northeast Asia, northern South America, Iberian Peninsula, and central Asia, whereas small areas of negative discharge and positive runoff are found over East and central Asia, northern North America, and southern Australia.  Ensemble anomaly in annual mean precipitation from the climate models (a) and in surface runoff from hydrological model PCR-GLOBWB (b). Areas for which more than two models (out of eight) disagree with the sign of the anomaly in precipitation and runoff are hatched. Also shown (c) are areas where the sign of anomalies in precipitation and runoff coincides or diverges, in a discrete color scale; white indicates that either percentage anomaly is comprised between −2% and +2%; gray indicates no runoff values.

Annual and Seasonal Discharge of Main Rivers
Annual hydrographs show much higher discharge at the river mouth in the LIG than PI for basins in Africa: Nile, Niger, and Jubba-Shebelle ( Figure 4). LIG discharge is clearly higher almost year-round for the Congo, Irrawaddy, Lerma-Santiago, and Douro. It is clearly lower for rivers Flinders, Murray-Darling, Zambesi, Orange, Paraná, Amazon, Mississippi, Saint Lawrence, Danube, and Ob and slightly lower for the Yenisei, Rhine, and Mackenzie. A common pattern is that several rivers with higher LIG yearly discharge display somewhat negative anomalies in spring/summer: the Ganges-Brahmaputra-Meghna, Pearl, Lena, and Anadyr. The opposite happens for the Rhine and Mackenzie, where positive summer anomalies emerge amidst annually negative anomalies. Intermodel spread in the percent anomaly of seasonal and annual discharge varies across rivers and tends to be larger in regions with higher LIG precipitation, runoff, and discharge: in the North African monsoon and adjacent region (Nile, Niger, and Jubba-Shebelle), in the Indian monsoon (Ganges-Brahmaputra-Meghna, Irrawaddy), in the South American monsoon (Orinoco), and in northeast Asia (Lena, Anadyr). An exception is the large intermodel spread for the Murray-Darling, where precipitation and runoff anomalies also diverge across models (Figure 2).
For most basins, the seasonal pattern of discharge is similar across LIG and PI simulations, except from some rivers in the tropics, where strong changes in the timing of monsoonal precipitation and in the domains of monsoons (Scussolini et al., 2019) result in a different shape of the hydrograph. Seasonality changes notably for the Congo, where a bimodal LIG hydrograph substitutes a near-unimodal PI one. For several basins, the seasonal contrast in discharge is stronger in the LIG, for example, for the Orinoco, Douro, Rhine, Ganges-Brahmaputra-Meghna, Irrawaddy, Yangtze, and Pearl. Furthermore, in many basins the peak of annual discharge is shifted forward, by days up to weeks: in the Amazon, Orange, Zambesi, Mississippi, Saint Lawrence, Douro, Danube, Rhine, and Yangtze.
To represent changes in the annual timing of discharge, we inspect changes in the timing of the annual center of mass of flow (CMF) (Stewart et al., 2005) (bottom lines in Figure 4). The CMF of the LIG is anticipated by as much as 42 days in the case of the Niger and postponed by as much as 36 days in the Irrawaddy. In many rivers of the Northern Hemisphere the CMF occurs later in the LIG (Mackenzie, Mississippi, Saint  Table S1. Discharge from hydrological model PCR-GLOBWB is shown in Figure S8 for reference.

SCUSSOLINI ET AL.
Lawrence, Douro, Rhine, Danube, Ob, Irrawaddy, Yenisei, and Pearl), in many cases reflecting enhanced precipitation in summer and early autumn. The Niger and the Yangtze are the only boreal rivers for which LIG CMF occurs substantially earlier. It must be noted that the Niger is also anomalous in that its LIG discharge is +200% larger. In the Southern Hemisphere, the CMF of LIG rivers shifts substantially only in the Orange and Zambesi, where it occurs later, and in the Mekong and Murray-Darling, where it occurs earlier.
Although the seasonal pattern changes widely in the Congo, its CMF is virtually the same. In other basins with much higher LIG discharge (Nile and Jubba-Shebelle), the CMF shifts by only 2 days.

Extreme Discharge and Floods
Global area-weighted average anomalies in annual maxima of flood metrics are more accentuated than anomalies in annual average precipitation, runoff, and discharge ( Figure 1). Anomalies in annual maximum of flood volume are larger than flood area. Globally, LIG simulations display +11% (with +5% to +21% intermodel spread) flood area and +21% (+10% to +43%) flood volume. Especially in the Northern Hemisphere, the positive anomalies are very strong, with +39% (+24% to +57%) flood area and +82% flood volume, with very large intermodel spread for the latter: +45% to +127%. For the Southern Hemisphere, anomalies are −17% (−20% to −15%) for flood area and −27% (−34% to −24%) for flood volume. The much larger anomalies of floods metrics cannot be attributed to changes in extreme precipitation (not shown), as these are even smaller than changes in mean precipitation. Geographic patterns of anomalies in flood area and flood volume at the 100-yr return period ( Figure S13) diverge in many regions from anomalies in runoff ( Figure 2) and discharge (Figure 3), mainly over northern Asia, North America, and India. Also, they are more spatially heterogeneous, pointing to the relevance of subbasin processes to modeled floods. To assess the link between precipitation intensity and hydrological extremes, we calculate annual maximum precipitation in consecutive intervals of 5, 15, and 30 days ( Figure S14). In the Northern Hemisphere, geographical coincidence in the sign of anomaly is higher between flood area (or volume) and 5-day maximum precipitation than between flood area (or volume) and average precipitation ( Figure S13). This indicates that, in the  Figure S5), for a set of return periods from 2 to 100 yr. Red and blue lines represent last interglacial and preindustrial climates, respectively. Individual models are in thin lines, and ensemble means in thick lines. Discharge is calculated for each basin in proximity of the river mouth; flood area is summed for the whole basin.

Geophysical Research Letters
model, extremes in boreal flooding are driven by intense precipitation over multiple days rather than by average precipitation.
For discharge and flood area in the same river basins analyzed above, the ensemble-averaged extreme values increase from return periods of 2 to 100 yr in a similar fashion ( Figure 5; flood volume is omitted because the shape of its curves is comparable to flood area). For most basins, patterns of extremes of intense precipitation are very similar for 5-, 15-, and 30-day intervals and mostly also resemble patterns of extreme discharge and flood area (Figure S.15). Notable deviations from this similarity emerge: (1) in the Mississippi, where LIG discharge and flood are lower than PI, whereas extremes of intense precipitation are virtually the same across simulations; and (2) in the Saint Lawrence and Flinders, where LIG discharge extremes are slightly lower, whereas LIG flood area extremes are much lower.

Comparison to Proxies of Discharge/Runoff
We revise proxy studies that report patterns of LIG runoff and discharge, to compare to the modeling results. We select mostly records of marine sediments under the influence of either discharge downstream from main river mouths or broadly of runoff from the adjacent land. Often, the interpretation of such records is dubious, because of the presence of multiple mechanisms that may concur to explain the proxy signal, such as strong sea level variations of the glacial-interglacial cycles (e.g.,  and potential shifts in the location of river mouths that are poorly constrained topographically. We find published proxies for 32 sites, of which 23 are suitable for comparison with models (Table S1). The modeled discharge agrees with the proxies on the sign of LIG anomaly at 15 sites: the west coast of North America; the north and west coasts of South America; at several sites on the northern, eastern, and southern coasts of Africa; at the southeast of India; and at western Australia (Figure 3). Disagreement occurs over the Mississippi basin and downstream its river mouth, at the northern coast of South America (where model resolution is likely inadequate to capture small-scale gradients), the western coast of Northern Africa, the eastern Mediterranean region, over the Limpopo basin, and at a site in Indonesia.

Discussion and Conclusions
Because reconstructing climatic and hydrological signals from proxies is hard, it is tempting to extend the interpretation of proxies of, for example, precipitation, to also represent hydrological variables as runoff, discharge, and flood, or conversely. Although our model results indicate that for most areas the sign of anomaly in those variables coincides, the discrepancies that emerge for nonnegligible land portions (Figures 2, S11, and S13) provide a cautionary message for the interpretation of proxy records. On the other hand, areas where modeled variables diverge and areas of stronger anomaly could inform site selection for proxy reconstructions that help shed light into mechanisms of LIG hydroclimate.
In the only modeling studies that offer points of paleoclimatic comparison for our results, Aerts et al. (2006) and Ward et al. (2007) simulated discharge for major rivers during the mid-Holocene, when insolation and warming anomalies were geographically comparable, though smaller in magnitude, to those of the LIG. Their annually integrated anomalies agree in sign with our LIG ones for the Ganges, Congo, Lena, Murray-Darling, Amazon, Zambesi, and Nile. They slightly disagree for the Danube, Rhine, and Mekong and completely disagree for the Mississippi. Attributing discrepancies with those results is unattainable because the effect of different climate forcing and very different modeling framework cannot be separated. By virtue of the possible similarities in the geographic pattern of temperature anomalies, another term of comparison for LIG discharge are the ensemble mean results of simulations for the end of this century (2070-2100) under high-emission scenarios, although key differences in forcing between these climates (Bony et al., 2013) limit the merit of this comparison. Our LIG positive discharge anomalies agree with future projections from CMIP3 models (van Vliet et al., 2013) over northeast Asia, Alaska, central Africa, and the Arabian Peninsula. LIG negative discharge anomalies agree over southern Africa and parts of North America and of Europe. Glaring disagreement emerges over Northern Africa and the Mediterranean, tropical South America, Southeast Asia, and particularly the Himalaya region. Compared to CMIP5-based projections, anomalies in LIG discharge hydrographs are similar for Rhine, Ganges, and Yangtze and different for Amazon, Niger, Nile, Mississippi, Mackenzie, and Lena (Krysanova et al., 2017). LIG floods are lower over many parts of Europe, in agreement with CMIP3-based future projections (Arnell & Gosling, 2016) but in disagreement with projections based on downscaled CMIP5 results (Alfieri et al., 2015).

Geophysical Research Letters
Precipitation is an uncertain product of global climate models (e.g., Huang et al., 2020; see also Text S1). Nevertheless, precipitation from an almost identical ensemble of climate models has been validated by comparison to proxies (Scussolini et al., 2019). That study reports that climate models seem to reasonably reproduce the main features of LIG precipitation, especially in the Northern Hemisphere. Besides, comparison of proxies with a similar ensemble of models shows broad agreement in LIG Arctic amplification of warmer temperatures (Otto-Bliesner et al., 2020). Our proxy-model comparison on discharge and runoff supplies additional confidence in the models' capabilities of reproducing past climate change, suggesting that even for a cascade of climate-hydrological-hydrodynamic models there is more agreement than disagreement with proxies. It must be noted that there are inherent limits of this model-proxy comparison. Beyond those raised by Scussolini et al. (2019), an additional limit derives from the fact that LIG simulations are carried out with current maps of river network, land cover, permafrost, and sea levels. Although these were likely different in the LIG (e.g., the top of the permafrost may have been deeper in warmer boreal climate), improvements in this direction are postponed until global data sets for this period will emerge. Further, although we found more proxies than were found for Holocene discharge (Aerts et al., 2006), these are fewer than for LIG precipitation (Scussolini et al., 2019), and the number of suitable comparisons is too small to draw conclusions about the geographic pattern of model-proxy agreement. Nevertheless, with growing interest and with accelerated development in the fields of paleohydrology and paleo floods (Baker, 2013;Wilhelm et al., 2017Wilhelm et al., , 2018, it is realistic that the coverage and resolution of proxies for runoff, discharge, and flooding will increase, enabling more extensive comparison with the models presented here. Additionally, the above mentioned forms of validation with proxies enable selecting models that best reproduce observations for the LIG and that can therefore be assumed to better capture the internal feedbacks of the climate system in a warmer climate. With ad hoc methodologies, this could in turn enable differential weighting of models in the context of ensemble studies (Eyring et al., 2019).
Lastly, our hydrological results inform the study of the LIG history of ecosystems and humans at the continental scale. For example, humid LIG conditions have been inferred for sectors of North Africa that were crucial to prehistoric human migrations (Larrasoana et al., 2013). The timing of human dispersal from Africa is not well constrained, with current hypotheses based on genomics and archeology suggesting migrations in the time window 130-50 ka, via the Bab el Mandeb strait or the Sinai Peninsula (López et al., 2015), and pioneer hydrodynamic studies indicating possible LIG migration to the Mediterranean along a western route (Coulthard et al., 2013). For 127 ka, increased precipitation emerges from proxies and models over the Sahel, Sahara, Arabian Peninsula, Levant, and Middle East, with further-reaching and longer monsoon seasons (Scussolini et al., 2019). We here simulate higher LIG annual runoff and discharge over North Africa and the southwestern Arabian Peninsula. All this points to the existence of a river network and to greening of key regions that may have sustained the path of human tribes migrating across and out of the continent.

Data Availability Statement
The results of the hydrological and hydrodynamic models that form the basis for the analyses here presented are available as netCDF files at public repository (https://doi.org/10.5281/zenodo.3677737).