Lacustrine plains, lake terraces, and lake shorelines are the landform relics of lake evolution. The lake terrace relics and corresponding lacustrine deposits, which represent the stable stage of lakes during their expansion period, are direct evidence of the evolutionary history of lakes and for reconstructing high lake levels. The lakes on the Ordos Plateau have been retreating since the late Quaternary. Hongjian Nuur and Bojianghaizi developed three lake terraces (Figure 2). T2 and T3 developed before the late Pleistocene, and T1 developed in the early and middle Holocene (Yuan et al., 1987; Shi, 1991). Holocene lake terrace and lacustrine strata in the MUD are generally no more than 10 m higher than the modern lake surface. The lacustrine sediments are gray-white or gray-green silty sand and fine sand with horizontal bedding. The sedimentary texture is loose with poor cementation. The lake terraces are stepped terrain with a flat top and a steep front edge. Lacustrine strata remain partially on the terrace. Based on remote sensing image interpretation, lake basins of Hongjian Nuur (HJN), Dalengturu Nuur (DLTR), Baozhai Nuur (BZN), Sharibuduyin Nuur (SRBDY), and Bojianghaizi (BJHZ) were selected for field investigation and sampling (Figure 1B, No.1–5).

HJN is the largest lake in MUD. Shorelines and lake terraces are clearly visible in south HJN. The shorelines formed in recent decades are roughly parallel to the modern shoreline (Figures 2A,B). The top of the lacustrine strata in T1 is about 9 m above the modern lake. The sediment is gray-white very fine sand. The terrace of DLTR is in the south of the lake. The terrace is about 5.5 m higher than the modern lake, slightly inclined toward the lake (Figure 2C). The lacustrine strata is gray-green fine sand with clay, about 17 cm in thickness (Figure 2D). The Holocene terrace of BZN is located on the southeast beach of the lake about 9 m above the modern lake. The lacustrine layer is gray-white silt sand with well-developed horizontal bedding (Figure 2E). The lake terrace of SRBDY, located in the east of the lake, is about 8 m higher than the modern lake. The lacustrine strata at the top is about 40 cm, with gray-white fine sand and visible and well-developed horizontal bedding (Figure 2F). The T1 in the north of Bojianghaizi is about 7 m higher than the modern lake with gray-green and gray-black lacustrine clay (Figure 2G).

Two quartz Optical Stimulated Luminescence (OSL) dating samples were collected in the DLTR section. Pure quartz grains 38–63 μm size were extracted from the middle part of the sample tube using the procedure employed by the OSL Laboratory of the State Key Laboratory of Earth Surface Processes and Resource Ecology, Beijing Normal University. The SAR-SGC method was used to determine the equivalent dose (De) (Lai and Ou, 2013).The measurement equipment was a Risø DA-20 TL/OSL reader equipped with blue diodes (k = 470 ± 20 nm) and a 90Sr/90Y radioactive beta source. The concentrations of Uranium (U), Thorium (Th), and Potassium (K) were measured by inductively coupled plasma mass spectrometry (ICP-MS). The water content was measured by weighing the samples before and after drying. Finally, we obtained the boundary ages of the lacustrine deposit in the DLTR section.

The lake area is governed by lake inputs and outputs. The lake level remains stable when the water inputs equal outputs. Lake inputs are lake surface precipitation and drainage runoffs, whereas the outputs are water surface evaporation and groundwater discharge. Based on the reconstructed lake area, the lake hydrologic model can be used to calculate the precipitation of the corresponding period. The model used in this study followed Goldsmith et al. (2017): P × A C × F R + P × A L = E × A L + G (1) where P is mean annual precipitation(m); A C is the catchment area (m²) where the area of the lake is subtracted from the area of the catchment; F R represents the fraction of the watershed precipitation that reaches the lake (including surface and groundwater runoff); A L is the area of the lake (m²); E represents the mean annual evaporation of the lake surface (m); and G represents groundwater discharge from the lake per year (m³).

The catchment area ( A C ) was calculated from topographic maps with contour lines of 5 and 10 m. Considering the inter-annual fluctuations of the lake level (Xu et al., 2019), the modern area of the lake ( A L ) was averaged from data of two periods: a topographic map drawn in the 1970s and a geomorphic map of the MUD (Wu et al., 2017). The precipitation was interpolated from the mean annual precipitation of meteorological stations in the MUD for 56 years (1961–2016). The highstand area of five lakes was calculated from topographic maps according to the elevation of the lacustrine strata. The watershed area of the lake highstand was the same as the modern catchment area because the elevation of each watershed was much higher than the modern lake and Holocene highstand. All the catchments are closed basins with no hydraulic connection between them. Thus, groundwater discharge from the lake ( G ) can be neglected for the hydrologic model. The meteorological data of each basin was obtained by interpolation.

The evaporation of the lake surface was estimated using a modified Penman equation (Penman, 1948; Sene et al., 1991): E d = Δ H / L + γ E a Δ + γ (2) where E d is evaporation per day (m/day); ∆ is the slope of the saturation vapor pressure curve at air temperature t (hpa/°C); H is net radiation, J/(m²∙day); L represents latent heat of vaporization (J/kg); γ is psychrometric constant (hpa/°C); and E a is atmospheric conductance (m/day).

The climatological data (radiation, temperature, air pressure, vapor pressure, wind velocity, and sunshine duration) were obtained from the Chinese Meteorological Agency for the period 1961–2016 for meteorological stations located in the MUD. The mean annual evaporation of each meteorological station was calculated by summing evaporation per day ( E d ). The evaporation of each lake surface ( E ) was mathematically interpolated from these stations. The runoff fraction ( F R ) can be calculated using the closed-basin hydrological model. In addition, it can be estimated by the Budyko relationship (Budyko, 1974; Koster et al., 2006): F R = 1 − E c / P = 1 − [ D ( tanh ⁡ 1 / D ) ( 1 − cosh ⁡ D + sinh ⁡ D ) ] 0.5 (3) where E c is catchment evaporation calculated as D = R n e t /P L , where D is the dryness index, R n e t is mean annual net radiation, P is annual precipitation, L is latent heat of vaporization, and R n e t / L is potential evaporation. Potential evaporation can be measured by the radiative balance or by pan evaporation (Koster et al., 2006; Goldsmith et al., 2017). It also can be approximately estimated using the Penman equation (Yang et al., 2012; Li et al., 2020).In our calculation, the Penman equation was used.

In the modern lake hydrologic model, the mean annual evaporation of meteorological stations in the MUD was calculated using the Penman equation. Then, the evaporation ( E ) of the five lake basins was calculated by mathematical interpolation. The modern runoff fraction ( F R ) of each basin was calculated using the lake hydrologic model and Budyko relationship, respectively. The results show that the two runoff fractions ( F R 1 and F R 2 ) are similar, indicating that both the modern lake hydrological model and Budyko relationship produced suitable estimates.

In this paper, first, the model was applied to HJN, DLTR, BZN, SRBDY, and BJHZ to reconstruct precipitation during the lake highstand. Second, the paleolake areas of other drainages (No. 6–24 in Figure 1B) were estimated using this model under the reconstructed precipitation pattern. In the paleolake hydrologic model, the evaporation was also calculated using the Penman equation. Considering the influence of wind speed and temperature on evaporation, we estimated that the evaporation during the lake’s highstand was at most 15% higher than the modern value. For the paleolakes, precipitation (or paleolake area) and runoff fraction could not be constrained independently. Therefore, we used both the closed-basin lake hydrologic model and the Budyko relationship to assess these parameters.

A series of precipitation data was extracted from published articles to reconstruct the paleoprecipitation pattern in MUD combining with reconstructed precipitation by the model in section 3.2. The published precipitation data were extracted according to the following methods:1) the paleoprecipitation was extracted immediately from articles with exact data; 2) the paleoprecipitation was calculated using modern precipitation (1961–2016) multiplied by change percentage given in the articles; 3) with regard to data in precipitation series, data from the chosen period were digitalized and averaged. With the method above, forty precipitation data was extracted in the MUD and its surrounding region in north China.

Based on field investigations, the highstand of five lakes was reconstructed using a topographic map and a digital elevation model (Figure 3). The area of the paleolakes is 1.65–3.56 times that of modern lakes, with the smallest ratio in HJN. The total area of paleolakes in the five catchments is 158.2 km², which is 2.05 times that of modern lakes. There was only one lake during the highstand period in the DLTR and SRBDY drainage. The drainage basins in HJN, BJHZ, and BZN contain more than one lake. In each basin, the elevation of the watershed between lakes is the upper limit of the upstream paleolake level. We also inferred that paleolakes existed in some depressions during the highstand period. For example, the Shejiahaizi (haizi means lake in the dialect of study area), located in depression northeast of HJN, was a lake that regressed decades ago according to the topographic map and remote sensing image. During the Holocene highstand period, Shejiahaizi was probably an independent lake that was ∼3 m higher than the HJN based on lacustrine landform. In the BZN drainage, lake Baozhai Nuur and the adjacent Baga Nuur were connected to a unified lake due to rising lake levels. In the west of the drainage, there are two depressions where paleolakes may develop. In the BJHZ drainage, two adjacent lakes were connected during the high lake level period.

The quartz optically stimulated luminescence (OSL) ages of the lacustrine strata in DLTR are listed in Table 1. The quartz OSL ages confirmed the lake highstand occurred during the early to middle Holocene.

The modern precipitation ( P ), modern and highstand evaporation ( E ), catchment area ( A C ), and runoff fraction ( F R ) required in the hydrologic model are listed in Table 2. The interpolated modern mean annual evaporation ranged from 995 to 1,065 mm with a rough trend of increasing from east to west. The evaporation of the paleolakes was calculated using modern values multiplied by 1.15. The runoff fraction ( F R ) of the paleolake drainages during lake highstand was estimated by Budyko relationship. They ranged from 5% to 8%, higher than at present in each catchment. The runoff fraction ( F R ) in both periods shows a decreasing trend from southeast to northwest.

The precipitation of five lakes during the highstand was estimated using the lake hydrologic model (Table 2). The reconstructed precipitation of HJN in the eastern MUD was 500 mm. In the middle MUD, the precipitation of BZN and BJHZ was 471 and 472 mm, respectively. The precipitation of DLTR and SRBDY in the west was 443 and 444 mm, respectively. The reconstructed precipitation was 120–165 mm higher than the modern value. The increase in precipitation was 32–60%, and the increase in the west was significantly higher than that in the east.

Precipitation patterns in the MUD during the early to middle Holocene were reconstructed by synthesizing precipitation data in section 4.2 and published data. A total of forty-five precipitation data from the early to middle Holocene in the MUD and surrounding areas were used (Supplementary Table S1). The isohyet of the MUD was reconstructed by interpolation (Figure 4). The results show that the precipitation during the early to middle Holocene in the MUD decreased from about 550 mm in the southeast to about 350 mm in the northwest, which is consistent with present levels. The 400 mm isohyet moved 130–170 km to the northwest, roughly consistent with the modern 250 mm isohyet. The moving distance in the southwest is slightly smaller than that in the northeast. The amplification of precipitation gradually increased from southeast to northwest. Due to the difference in amplification, the precipitation gradient in much of the MUD is lower than at present. The large gradient is distributed near the 400 mm isohyet in the west MUD.

The lake highstand of five drainage areas was reconstructed based on lacustrine landform and strata in section 4.1. Based on the reconstructed precipitation pattern in the MUD, the paleolake area of each drainage (No. 6–24) during the early to middle Holocene was estimated. First, the runoff fraction (F_R2) of each drainage basin was estimated using the Budyko relationship. Second, the lake area of each basin was calculated using the lake hydrologic model (Table 2). The total area of reconstructed paleolakes during the early to middle Holocene was 896.1 km², which is 2.96 times that of modern lakes. The lake area of each drainage is two to five times that of modern lakes with the largest increase in the southwest MUD. The area of the paleolake in the Wuhudong Nuur, Gouchi, and Beidachi basins (No. 22–24 in Table 2) are 4.57, 4.84, and 4.97 times that of modern lakes, respectively.

Based on the reconstructed lake area in each drainage (No. 6–24), we also attempted to reconstruct the paleolake pattern. For the catchment comprising one single lake, the paleolake expanded by topographic contour around the modern lake based on the calculated area. In the basin comprising a series of catenulated lakes, the elevation of the watershed between lakes was regarded as the altitude of the paleolake upstream. The upstream paleolake area was calculated using closed contour and the corresponding altitude. The downstream lake area was calculated as the total area minus the upstream lake area. Due to the rise in lake levels, adjacent lakes may have become unified lakes. Therefore, the modern 108 lakes decreased to eighty-four during the early to middle Holocene. In addition, there are small, closed depressions in certain catchments. These depressions were inferred as paleolakes during the Holocene highstand. A total of forty-three inferred lakes with an area more than 0.1 km² were identified. According to the method above, there were 127 paleolakes in the closed flow area of MUD during the early to middle Holocene. The paleolake pattern is shown in Figure 4.

As one of the required parameters in the lake hydrologic model, open water evaporation is typically estimated by observing different pans. However, the observation data of a 20 m² evaporation pan, which is approximately equal to open water evaporation, is rare in the MUD. The observation data of D20 and 601B pans are significantly higher than open water evaporation in arid and semi-arid areas because of the pan wall. In certain regions, open water evaporation is estimated using observation data of pans multiplied by a conversion coefficient (Linacre, 1994). However, the conversion coefficient in the MUD was difficult to constrain. Therefore, a modified Penman equation was used to estimate the evaporation of the lake surface (Penman, 1948; Sene et al., 1991). The Penman equation has been used to calculate evaporation in arid and semi-arid China (Li et al., 2001; Yang et al., 2010). The difference between calculated values and measured data in 20 m² evaporation pans was no more than 10% according to experiments conducted in the arid regions of Xinjiang, in western China (Zhang and Zhou, 1992).

The observed evaporation of two pans and calculated evaporation using the Penman equation are listed in Table 3. The ratios E D 20 / E and E 601 B / E are 0.4–0.52 and 0.81–0.97, respectively. Both ratios decreased from the semi-arid southeast to the arid northeast. The average conversion coefficient between D20 pan and open water evaporation in arid northwest China and Inner Mongolia Plateau was estimated at 0.44 and 0.52, respectively (Xie and Wang, 2007). The E D 20 / E ratio in the MUD was roughly between these two coefficients. The conversion coefficient between 601B pan and 20 m² pan in North China ranged from 0.83 to 0.93, which is close to the E 601 B / E ratio in the MUD. The decreasing trend from semi-humid Henan to semi-arid Inner Mongolia was also consistent with the MUD (Qian et al., 1998). This indicates the evaporation calculated using the Penman equation was approximately equal to the open water evaporation.

In the paleolake hydrologic model, the evaporation of lakes was also estimated using the Penman equation. The temperature during the early to middle Holocene was about 1°C higher than modern times (Fang et al., 2011). Based on the sensitivity of evaporation to temperature change, a 1°C rise would produce evaporation rate changes of about 5–6% in the MUD. We also examined the sensitivity of evaporation to wind speed change. A 20% change in wind speed can lead to about a 10% change in evaporation, which is nearly consistent with the examination in Dali Lake in eastern Inner Mongolia (Goldsmith et al., 2017). Thus, considering the possible change in temperature and wind speed, the evaporation was assumed at 15% higher than the modern value at most.

The lake hydrologic model is one of the common methods for reconstructing precipitation. The precipitation reconstructed in HJN was 500 mm. The result is approximate to the value 510 ± 25 mm estimated from pollen data from the Hejialiang section about 15 km east of the lake (Chen et al., 1993). The reconstructed precipitation in BZN was 471 mm. The rainfall calculated based on pollen assemblage in Taobao and Wushen profiles was 470 ± 25 and 461.8 mm, respectively (Shi, 1991; Chen et al., 1993). These two profiles are about 7 and 15 km to the southeast of BZN, respectively. This indicated that the lake hydrologic model and pollen assemblage produced suitable estimates. The mean annual precipitation of the HJN lake was ∼380 mm. The precipitation during the early to middle Holocene was ∼32% higher than present. The result is consistent with reconstructed results near the 400 mm isohyet in the MUD and southeast Inner Mongolia Plateau (Jiang et al., 2006; Xu et al., 2010). However, the result was much smaller than the estimated precipitation of Qinghai Lake (Shi et al., 1992). This may reflect that the migration distance of the monsoon rain belt varied in different regions of the monsoon edge area.

In the EASM area of China, the relative increase in precipitation during the early to middle Holocene varied from ∼10 to 15% in the humid southeastern coastal area to more than 70% in the arid northwest area (Yang and Williams, 2003; Stebich et al., 2015; Liu and Li, 2017; Li et al., 2018; Wu et al., 2018). The reconstructed relative increase in the MUD rose from ∼32% in the southeast to ∼60% in the northwest. The rising trend is consistent with the whole EASM area in China. However, the extent of relative increase ranges significantly from southeast to the northwest, indicating that the monsoon edge area was extremely sensitive to climatic change. During the early to middle Holocene, the precipitation gradient decreased significantly in the MUD due to the difference in precipitation increase. This probably reflects the monsoon rain belt and the monsoon edge area migration to the northwest. In consequence, the MUD was stably dominated by the EASM. The spatio-temporal pattern of climate showed a relatively stable, warm, and humid climate for millennia. It can also be confirmed from the peat and paleosol development process in the MUD. The development of peat and lacustrine in the MUD since the last deglaciation indicates the regional climate fluctuated considerably until 8.5 ka. Thereafter there were humid intervals until 6ka, during which lacustrine sediments were deposited (Liu et al., 2018b). The paleosol was generally regarded as the indicator of warm and humid climate. In the aeolian sand/paleosol sequence of the MUD, the duration of paleosol is much longer than the aeolian sand during the early to middle Holocene (Gao et al., 1993; Ding et al., 2021).

The lakes in northern China generally reached their highstand due to the increase in precipitation in the early and middle Holocene. The lake landforms and lacustrine strata of HJN and BJHZ are direct evidence of the development of high lakes in the study area. The decrease in the grain size of lacustrine deposits in the lake cores of the Baahar Nuur and Qigai Nuur in north MUD indicates an increase in lake level during the early and middle Holocene (Guo et al., 2007; Sun and Feng, 2013). The clay deposits in Chagan Nuur in the western MUD indicate a stable deep-water sedimentary environment of ∼9–5 ka (Yang, 1997). Clay sediments in Hetongchahagan Nuur in the western MUD indicate a stable deep-water sedimentary environment at∼9—6 ka (Bernasconi et al., 1997). Grain size analysis of the DSGL profile in the southern MUD showed that the Salawusu Lake began to form in the early Holocene and reached its maximum in the middle Holocene (Liu et al., 2018b). These lake records above indicate that the lake highstands during the early to middle Holocene in the MUD were almost simultaneous and stable.

The total area of lakes in the closed drainage area is 2.96 times that of modern lakes. The relative height between the highstand and the modern lake surface is about 5–9 m. The increase of height and area was smaller than the lakes in Badain Jaran Desert and Tengger Desert (Yang and Williams, 2003; Zhang et al., 2004). Most of the lakes in the study area locate in bedrock valleys or aeolian bedrock depressions. The topography of the lake basins is mainly controlled by bedrock gullies. The modern aeolian activity has little impact on topography. The stable terrain is conducive to the reconstruction of paleolake pattern. According to the reconstruction results, the number of paleolakes increased along with the rise of the lake level and the expansion of the lake area. The adjacent lakes in certain basins connected and formed unified lakes due to the rise in water level. In addition, closed depressions were extracted from topographic maps and remote sensing images. Geomorphic relics, such as lake shorelines, can be distinguished around the depressions. Salt crystals on the surface of depressions indicated paleolakes may have existed (Figure 5). We inferred these closed depressions formed lakes due to the rise of precipitation during the early to middle Holocene. In each catchment, the altitude of the lake highstands was lower than the watershed between drainages. The pattern of closed drainage area is consistent with the present. There is no evidence from the lake landform, lacustrine strata, and climatic conditions supporting the formation of a unified megalake in the closed drainage area during the early to middle Holocene.