A first investigation of hydrogeology and hydrogeophysics of the Maqu catchment in the Yellow River source region

The Tibetan Plateau is the source of most of Asia's major rivers and has been called the Asian Water Tower. Detailed knowledge of its hydrogeology is paramount to enable the understanding of groundwater dynamics, which plays a vital role in headwater areas like the Tibetan Plateau. Nevertheless, due to its remoteness and the harsh environment, there is a lack of field survey data to investigate its hydrogeology. In this study, borehole core lithology analysis, soil thickness measurement, an altitude survey, hydrogeological surveys, and hydrogeophysical surveys (e.g. magnetic resonance sounding – MRS, electrical resistivity tomography – ERT, and transient electromagnetic – TEM) were conducted in the Maqu catchment within the Yellow River source region (YRSR). The hydrogeological surveys reveal that groundwater flows from the west to the east, recharging the Yellow River. The hydraulic conductivity ranges from 0.2 to 12.4 m d⁻¹. The MRS sounding results, i.e. water content and hydraulic conductivity, confirmed the presence of an unconfined aquifer in the flat eastern area. Based on TEM results, the depth of the Yellow River deposits was derived at several places in the flat eastern area, ranging from 50 to 208 m. The soil thickness measurements were done in the western mountainous area of the catchment, where hydrogeophysical and hydrogeological surveys were difficult to be carried out. The results indicate that most soil thicknesses, except on the valley floor, are within 1.2 m in the western mountainous area of the catchment, and the soil thickness decreases as the slope increases. These survey data and results can contribute to integrated hydrological modelling and water cycle analysis to improve a full-picture understanding of the water cycle at the Maqu catchment in the YRSR. The raw dataset is freely available at https://doi.org/10.17026/dans-z6t-zpn7 (Li et al., 2020a), and the dataset containing the processed ERT, MRS, and TEM data is also available at the National Tibetan Plateau Data Center with the link https://doi.org/10.11888/Hydro.tpdc.271221 (Li et al., 2020b).

1 Introduction

With a huge amount of water storage, the Tibetan Plateau (TP) acts as the “Water Tower of Asia” (Qu et al., 2019; Wang et al., 2017), recharging many major Asian rivers including the Salween, Mekong, Brahmaputra, Irrawaddy, Indus, Ganges, Yellow, and Yangtze rivers (Immerzeel et al., 2009), feeding more than 1.4 billion people (Immerzeel et al., 2010) and promoting regional social and economic development (Xiang et al., 2016). Due to climate change, the TP has experienced accelerated temperature rise over the past decades (Huang et al., 2017). Since the 1950s, the warming rate over the TP ranges between 0.16–0.36 ^∘C per decade and rises to 0.50–0.67 ^∘C per decade from the 1980s (Kuang and Jiao, 2016). The retreating glaciers and snow cover, decreasing wetland area, and rising snow lines indicate that the hydrological system on the TP is undergoing profound changes (Kang et al., 2010; Li et al., 2019; Xu et al., 2016; Yao et al., 2013).

So far, the groundwater-related studies on the TP are mainly satellite-based, focusing on using the Gravity Recovery and Climate Experiment (GRACE) to estimate terrestrial water storage, which consists of surface and subsurface water (Haile, 2011; Jiao et al., 2015; Zhong et al., 2009). Among those studies, Xiang et al. (2016) separated the groundwater storage from terrestrial water storage observed by GRACE using four land surface models: the Community Land Model (CLM), Mosaic, Noah, and the Variable Infiltration Capacity (VIC) model of the Global Land Data Assimilation System (GLDAS) (Rodell et al., 2004), Climate Prediction Center (CPC) soil moisture data (Fan and Van Den Dool, 2004), and a glacial isostatic adjustment model.

An integrated surface-groundwater model is essential for improving understanding of different processes quantitatively (Graham and Butts, 2005). To set up an integrated surface-groundwater model, different kinds of data are needed for the parameterization of land surface and subsurface but also for atmospheric forcing and state variables that are required for model calibration and validation. Land surface data such as topography, land cover, and soil parameters can be obtained from digital elevation models (DEMs) and regional or global soil databases (Su et al., 2011; Zhao et al., 2018). Atmospheric forcing data, including precipitation, air temperature, wind velocity, and other variables, are available from regional or global meteorological datasets (Li et al., 2017; Su et al., 2013, 2020; Yang, 2017).

Subsurface data, like hydrogeological information (e.g. lithology, water table depth, hydrogeological parameters) and state variables (e.g. hydraulic heads and soil moisture content at root zone or of deeper layers), usually require in situ measurements. These hydrogeology-related data are the most difficult ones to acquire. Efforts have been made to develop the global map of permeability (Gleeson et al., 2014, 2011), hydraulic conductivity (Gupta et al., 2021; Montzka et al., 2017), groundwater table depth (Fan et al., 2013), and groundwater volume and distribution (Gleeson et al., 2016). However, global subsurface databases suffer large uncertainty, mainly because of spatial heterogeneity and scale effects that can be resolved only by field measurements. The TP lacks such measurements due to its remoteness and harsh environment (Yao et al., 2019). Lately, the Chinese Academy of Science (CAS) launched the CAS Earth Poles project (X. Li et al., 2020, 2021), which aims to address the bottleneck of polar (high-mountain cold region) data curation, integration, and sharing, which is expected to overcome some of the aforementioned difficulties.

The conventional way to acquire hydrogeological information in an unknown area is by augering or drilling boreholes and carrying out hydraulic tests, for example, pumping tests (Vouillamoz et al., 2012). However, due to the harsh environment of the TP and the high costs and time-consuming nature of the traditional hydrogeological survey methods, little hydrogeological work has been done on the TP.

The hydrogeophysical methods are up and coming in hydrogeological studies (Chirindja et al., 2016). They have been applied in various conditions, for example in wetlands (Chambers et al., 2014), rivers (Steelman et al., 2015), proglacial moraines (McClymont et al., 2011), karst regions (McCormack et al., 2017), and volcanic systems (Di Napoli et al., 2016; Fikos et al., 2012). Compared to other hydrogeophysical methods, such as seismic, gravity, and resistivity methods, magnetic resonance sounding (MRS) is the only method that is able to detect the free water in the subsurface directly (Lubczynski and Roy, 2003, 2004) and quantify hydrogeological parameters and water storage (Lachassagne et al., 2005; Legchenko et al., 2002, 2018; Lubczynski and Roy, 2007).

The MRS method is based on the nuclear magnetic resonance (NMR) principle, widely implemented in medical magnetic resonance imaging (MRI). In contrast to MRI, the MRS excitation is done at the Larmor frequency corresponding to the local magnetic field of the Earth, which propagates and gets attenuated through conductive subsurface media. Therefore, at a given position of the Earth, MRS signal response depends on the strength of the Earth's magnetic field, and on the strength of the natural and artificial noise (Lubczynski and Roy, 2004) but also on the subsurface resistivity. The electrical resistivity measurement is suggested to be jointly used with MRS (Braun and Yaramanci, 2008; Descloitres et al., 2007; Vouillamoz et al., 2002). Electrical resistivity tomography (ERT) is one of the predominantly employed hydrogeophysical methods to estimate the subsurface electrical resistivity (Herckenrath et al., 2012; Jiang et al., 2018). It has been widely applied together with MRS to explore regional hydrogeology (Vouillamoz et al., 2003; Descloitres et al., 2008; Pérez-Bielsa et al., 2012). Furthermore, the transient electromagnetic survey (TEM), referred to as the time domain electromagnetic (TDEM) method, also provides subsurface resistivity but is able to achieve deeper penetration than ERT. The ERT has already been used on the TP by Gao et al. (2019) and You et al. (2013) to investigate permafrost. Nevertheless, we have not seen the combined use of MRS, ERT, and TEM for hydrogeophysical surveys over the TP.

Investigations on various fields, such as geomorphology, climate change, glacier, and permafrost, have been done on the TP based on different DEMs. Zhang et al. (2006) analysed the geomorphic characteristics of the Minjiang drainage basin with SRTM (Shuttle Radar Topography Mission) data. Wei and Fang (2013) assessed the trends of climate change and the spatiotemporal differences over the TP from 1961–2010, with a generalized temperature zone–elevation model and SRTM. Ye et al. (2015) calculated the glacier elevation change in the Rongbuk catchment from 1974 to 2006 based on topographic maps and ALOS (Advanced Land Observing Satellite). Niu et al. (2018) mapped permafrost distribution throughout the Qinghai–Tibet Engineering Corridor based on ASTER (the Terra Advanced Spaceborne Thermal Emission and Reflection Radiometer) global DEM. However, different DEMs used in different studies may lead to potential inconsistencies for understanding relevant physical processes (Zeng et al., 2019, 2015). For the Maqu catchment, located within the TP and being the focus of this study, it is crucial to evaluate the accuracy of different DEMs to select the best one and eventually calibrate it using ground measurements since altitude is used to calculate hydraulic heads, which determine and control groundwater flow. Therefore, it is important to evaluate the accuracy of satellite DEM products with a ground-based real-time kinematic global positioning system (GPS-RTK); such works have not been common over the TP.

This study jointly uses hydrogeological and hydrogeophysical methods, including aquifer tests, MRS, ERT, TEM, and other necessary approaches at the Maqu catchment in the Yellow River source region (YRSR) on the TP. The paper focuses on field hydrogeological and hydrogeophysical surveys and corresponding datasets, aiming to fill the scientific data gap on the TP. In what follows, the study area is introduced in Sect. 2. Borehole core lithology analysis, soil thickness measurement, altitude surveys, hydrogeological surveys, and hydrogeophysical surveys are presented in Sect. 3. The results are documented and uncertainties are discussed in Sect. 4. Data availability is given in Sect. 5. Conclusions are made in Sect. 6.

2 Study area

The study area is a catchment (33^∘43^′–33^∘58^′ N, 101^∘51^′–102^∘16^′ E) in Maqu County, China (Fig. 1). It is located at the north-eastern edge of the TP, the first major bend of the Yellow River. Maqu County is regarded as the “reservoir” of the YRSR. The part of the Yellow River in contact with the Maqu catchment is about 36 km long, while the whole length passing Maqu County is 433.3 km. When the Yellow River flows through Maqu County, the annual run-off increases by 10.8×10⁹ m³, accounting for 58.7 % of the total run-off of 18.4×10⁹ m³ of the Yellow River in the YRSR (Wang, 2008). The Maqu catchment is characterized by a cold climate with dry winter and warm summer, assigned as Dwb in the updated Köppen–Geiger climate classification (Peel et al., 2007). The annual mean temperature is about 1.8 ^∘C, and the precipitation is around 620 mm annually. The catchment is covered by short grasses used for grazing by yaks and sheep. The elevation ranges between 3367 and 4017 m a.s.l. according to ALOS PALSAR (Phase Array type L-band Synthetic Aperture Radar) RT1 (high terrain correction resolution). There is a reservoir in the catchment (Fig. 1c) with functions of grassland irrigation and flood control.

Figure 1 The geographical location of Maqu catchment, Maqu County, and Ruoergai basin on the TP and geomorphologic map: (a) the geographical location and boundary of the TP (Zhang et al., 2014a, b); (b) the geographical location of the Maqu catchment, Maqu County, and Ruoergai basin (Li and Gao, 2019); (c) the geomorphologic map of Maqu catchment.

Based on the field survey of geomorphology and geology (Compton, 1962; Dackombe and Gardiner, 2020), the catchment can be divided into two parts, the flat eastern area and the western mountainous area. The western mountains are feldspathic quartzose sandstone and sandy slate with soil covered at the top, while the eastern, relatively flat part is composed mainly of sediments, such as alluvial deposits with intercalated aeolian units. The eastern part, together with its extension outside of the study area, is called the Ruoergai basin (Fig. 1b). Surface processes cause erosion, mixing, unmixing, and redistribution of alluvial materials within the thick alluvia accumulation in the eastern part. Geomorphological characterization was carried out in the Maqu catchment in 2018, and the three terraces were identified (Fig. 1c).

Some previous works have been done in or around the catchment. Su et al. (2011) monitored the soil moisture and soil temperature from 5 to 80 cm below the ground surface. Dente et al. (2012) assessed the reliability of AMSR-E (Advanced Microwave Scanning Radiometer for Earth Observing System) and ASCAT (Advanced Scatterometer) soil moisture products. Zheng et al. (2016) investigated the impacts of Noah model physics on catchment-scale run-off simulations. Zeng et al. (2016) combined the in situ soil moisture networks with the classification of climate zones to produce the in-situ-measured soil moisture climatology at the plateau scale. Zhao et al. (2018) studied the soil hydraulic and thermal properties of the 0.8 m top soil column. Zhuang et al. (2020) blended the surface soil moisture data from satellites, land data assimilation, and in situ measurements with the constraint of in situ data climatology and estimated the root zone soil moisture by scaling the blended surface soil moisture product. The present research focuses on the hydrogeological and hydrogeophysical aspects, complementing previous studies.

3 Materials and methods

Figure 2 shows the fieldwork workflow towards establishing a hydrogeological conceptual model, which is beyond the scope of this study and will be presented in a future study. The workflow includes the borehole core lithology analysis, soil thickness measurement, an altitude survey, hydrogeological surveys, and hydrogeophysical surveys (Table 1 and Fig. 2). Yellow boxes in Fig. 2 represent the fieldwork, and green boxes represent the results of fieldwork, which will finally contribute to the set-up of the hydrogeological conceptual model. The obtained information on lithology, soil thickness, and ground surface elevation provides basic knowledge in the study area. Hydrogeological measurements of water table depth and hydraulic conductivity provide important input that can be used to deduce the direction and rate of regional groundwater flow, which is beyond the scope of this study. For hydrogeophysical results, magnetic susceptibility measurements ensure the suitability of applying MRS, while underground resistivity measurements from ERT are used in the inversion of MRS for deriving water content and transmissivity. TEM provides deeper resistivity information than other methods and estimates the depth of the hydrogeological boundary. The locations of the surveys and measurements are shown in Figs. 3, 4, and 5.

Table 1 Methods, equipment, and timing for carrying out relevant measurements as in Fig. 2.

^*40 water table depths were measured instantaneously in 34 boreholes in 2018 and 2019.

Figure 2 Fieldwork workflow for setting up a hydrogeological conceptual model at Maqu catchment, where italics represent indirect technique (e.g. inversion type of retrieval) with unknown uncertainty, regular bold letters represent direct technique with low uncertainty, and regular letters do not convey uncertainty information.

Figure 3 Locations of the hydrogeological surveys, ground surface elevation measurements, core lithology analysis, and soil thickness measurements. The numbers from 1 to 46 (due to limited space, several numbers are not shown in the figure) indicate the measurement sequence of GPS-RTK.

Figure 4 Panels (a), (b), (c), (d), (e), and (f) are the exact locations of soil thickness measurements at the sites a, b, c, d, e, and f, respectively, shown in Fig. 3, in the .KML-formatted image from © Google Earth.

Figure 5 Location of hydrogeophysical surveys.

3.1 Borehole core lithology

The borehole core lithology is helpful in terms of understanding the formation of an area, assessing geometry of aquifers, and estimating hydrogeological parameters. Some boreholes were available for water table depth measurement in the study area, but information of borehole core lithology was only available in one borehole (ITC_Maqu_1; shown in Fig. 3 as core lithology analysis) drilled in 2017 down to the depth of 32 m from the ground surface. The lithology of the cores was determined based on particle size analysis using the sieving method with mesh sizes of 60, 40, 20, 10, 5, 2, 1, 0.5, 0.25, and 0.075 mm.

3.2 Soil thickness measurement

In the mountainous west, where transportation and set-up of hydrogeophysical surveys are difficult, we auger-sampled the thickness of the overlying soils (Fig. 3) to build the hydrogeological conceptual model and to validate simulations of spatially distributed soil thickness by the landscape evolution model LAPSUS (LandscApe ProcesS modeling at mUlti-dimensions and scaleS) (Schoorl et al., 2006, 2002) (relevant details will be presented in another study). The fieldwork was carried out at six sites (Fig. 4a–f). Measurements at sites 1 and 2 were conducted in 2018, while the rest were conducted in 2019. Soil thickness and slope of the ground surface were measured using an auger and a clinometer from Eijkelkamp Soil & Water Company (https://en.eijkelkamp.com, last access: 4 March 2021​​​​​​​). The exact measurement positions at each site (Fig. 4) were decided based on slope forms and surface pathways.

3.3 Altitude survey

The accuracy of ground surface elevation is crucial for the assessment of hydraulic heads, hydraulic gradients, and also groundwater flow and its direction; therefore it is also important for groundwater modelling. The dynamic type of GPS positioning technique, GPS-RTK, is able to achieve centimetre-level accuracy in real time, in terms of geolocation and elevation. The GPS-RTK instrument CHCNAV T4 from Shanghai Huace Navigation Technology Limited (https://www.chcnav.com, last access: 4 March 2021), with a vertical accuracy of 3 cm and a horizontal accuracy of 2 cm, was employed to measure elevations in 2019. Among the 46 ground surface elevation measurements made in total, 33 were located in the flat eastern area, and 13 were located in the mountainous area (Fig. 3). The ground data were intended to evaluate seven DEM datasets (Table 2) to select the most accurate DEM to be applied for the calculation of hydraulic heads where the ground-based altitude survey is not available and also as the top model boundary in groundwater modelling (that is not part of this paper). The seven DEMs are all open-access and were downloaded from websites of the United States Geological Survey (USGS), Japan Aerospace Exploration Agency (JAXA), and Alaska Satellite Facility (ASF).

Table 2 Seven DEM datasets.

3.4 Hydrogeological surveys

3.4.1 Water table depth measurement

Water table depth information is important for hydrology and hydrogeology. By subtracting the water table depth from ground surface elevation, a hydraulic head is obtained. A set of hydraulic heads distributed over the study area can be used to determine the regional groundwater piezometric map to enable a general understanding of the groundwater flow system in the study area. We measured 40 instantaneous water table depths in 34 boreholes during 5–8 August 2018 and 20 August–5 September 2019 using a dipper (Fig. 3). For each measurement, the depth was measured three to four times to ensure data quality. Eight level-loggers were installed to monitor the long-term groundwater level fluctuation (see Data.xlsx in repositories), but the data are not available yet.

3.4.2 Aquifer tests

Aquifer tests, including pumping tests and slug tests, were conducted to obtain aquifer hydraulic conductivity (Fig. 3). The first pumping test was done in 2017, in the borehole ITC_Maqu_1, where core lithology information was available. The pumping rate was constantly 55.6 m³ d⁻¹, measured with a flowmeter, and the pumping duration was about 30 min. The pumping rate was limited because the borehole ITC_Maqu_1 could easily collapse if the pumping rate were too high. The water level became stable soon after the start of pumping and was recorded every minute using a data logger (TD-Diver manufactured by Van Essen Instruments, with a measuring range of 10 m). Other tests were carried out in 2019, including two pumping tests and eight slug tests (Fig. 3). For the two pumping tests with the constant rate of 31.6 and 101.52 m³ d⁻¹, only water level recovery data were analysed because the water pump ran out of energy before reaching a stable water level. Also recovery tests were done in the eight slug tests, where in all the tests, groundwater levels were abruptly lowered by extracting 11.75 L of water. The recovering water levels were recorded every second or 2 s in the six slug tests and every 5 s or 20 s in pumping-recovery tests, all using a data logger (3001 Levelogger Edge manufactured by Solinst, with a range of 10 m).

The pumping test data acquired from the borehole ITC_Maqu_1 were analysed using the Boulton (1963) method as follows:

$$\begin{array}{}\text{(1)}& s={\displaystyle \frac{Q}{\mathrm{4}\mathit{\pi }T}}W,\end{array}$$

where s is drawdown (m), Q is pumping rate (m³ d⁻¹), T is transmissivity (m² d⁻¹), and W is a dimensionless variable obtained by Theis curve matching. The method can be used for pumping tests performed in an unconfined aquifer (isotropic or anisotropic) and for both fully or partially penetrating wells.

Slug test data were analysed using the Bouwer and Rice (1976) method for hydraulic conductivity in an unconfined aquifer as follows:

$$\begin{array}{}\text{(2)}& K={\displaystyle \frac{r^{\mathrm{2}}\ln \left(\frac{R_{\mathrm{e}}}{R}\right)}{\mathrm{2}L}}\cdot {\displaystyle \frac{\mathrm{1}}{t}}\cdot \ln \left({\displaystyle \frac{h_{\mathrm{0}}}{h_{\mathrm{t}}}}\right).\end{array}$$

In case of partially penetrating wells,

$$\begin{array}{}\text{(3)}& \ln \left({\displaystyle \frac{R_{\mathrm{e}}}{R}}\right)={\left[{\displaystyle \frac{\mathrm{1.1}}{\ln \left(H/R\right)}}+{\displaystyle \frac{A+B\ln \left[\left(D-H\right)/R\right]}{L/R}}\right]}^{-\mathrm{1}},\end{array}$$

where K is hydraulic conductivity (m d⁻¹), r is the radius of the borehole casing (m), R_e is the effective radial distance over which the head difference is dissipated (m), R is the radius measured from the borehole centre to the undisturbed aquifer (m), L is the length of the screen (m), t is time (d), h₀ is the water level at time 0 (m), h_t is the water level at time t (m), H is the length from the bottom of the well screen to the water table (m), D is the thickness of the aquifer (m), and A and B are dimensionless coefficients that can be obtained from empirical curves developed by Bouwer and Rice (1976).

Another two pumping test data were analysed using the Boulton and Agarwal (1980) method. The Agarwal method employs the following simple data transformations for recovery test data:

$$\begin{array}{}\text{(4)}& {\displaystyle }{S}_{\mathrm{r}}=h-h_{\mathrm{p}}\text{(5)}& {\displaystyle }{t}_{\mathrm{r}}=t-t_{\mathrm{p}},\end{array}$$

where S_r is the recovery drawdown (m), h is the head during the recovery period (m), h_p is the head at the end of the pumping period (m), t_r is the recovery time (d), t is the time since pumping started (d), and t_p is the duration of pumping (d). Applying the data transformations in the Boulton method accounts for the drawdown period and allows recovery plots to be analysed with curve matching (Thorne and Newcomer, 2002). Data were processed automatically in AquiferTest software (https://www.waterloohydrogeologic.com/products/aquifertest/, last access: 4 March 2021)

3.5 Hydrogeophysical surveys

3.5.1 MRS

MRS was conducted to define aquifer geometry and estimate transmissivity, hydraulic conductivity, and water content with depth that can be used to derive storage parameters. In total, 18 soundings (Fig. 5) were performed using the MRS instrument Numis Poly, the latest version of MRS equipment from the IRIS Instruments company (http://www.iris-instruments.com, last access: 4 March 2021). The Larmor frequency, measured with the proton magnetometer in the field, was set at 2241.8 Hz, and the inclination of the Earth's magnetic field was set at 52^∘ N. Square loops with a side length of 150 m or 100 m were used. Positions were measured with the UniStrong MG858s handheld GPS (global positioning system) instrument (http://www.unistrong.com, last access: 4 March 2021), with a horizontal and vertical accuracy of 30 cm.

To estimate hydraulic conductivity, the decay time constant T_d was used. There are three kinds of T_d: longitudinal decay time constant T₁, transverse decay time constant T₂, and free induction decay time constant $T_{\mathrm{2}}^{\ast }$. With the current instrument, only T₁ (actually an approximate value $T_{\mathrm{1}}^{\ast }$) and $T_{\mathrm{2}}^{\ast }$ were available. The Seevers equation (Seevers, 1966) (Eq. 6) and the Kenyon equation (Kenyon et al., 1989) (Eq. 7) were used for estimating hydraulic conductivity K (m d⁻¹):

$$\begin{array}{}\text{(6)}& {\displaystyle }K=C_p{\mathit{\theta }}_{\mathrm{MRS}}{T}_{\mathrm{d}}^{\mathrm{2}},\text{(7)}& {\displaystyle }K=C_p{\mathit{\theta }}_{\mathrm{MRS}}^{\mathrm{4}}{T}_{\mathrm{d}}^{\mathrm{2}},\end{array}$$

where C_p is the calibration coefficient, which is a lithology-dependent factor that needs to be calibrated from the pumping test (dimensionless). θ_MRS is the MRS-estimated water content (%). Compared to the Kenyon equation, the Seevers equation is more accurate (Plata and Rubio, 2008) and has been widely used (e.g. Legchenko et al., 2002; Vouillamoz et al., 2007; Nielsen et al., 2011) and is used in this study. Once K is estimated, the transmissivity T (m² d⁻¹) can be calculated using the equation

$$\begin{array}{}\text{(8)}& T=K\cdot \mathrm{\Delta }z,\end{array}$$

where Δz is the layer thickness (m) derived from MRS inversion.

MRS data were interpreted with an open-access software Samovar V6.6 from the IRIS Instruments company (http://www.iris-instruments.com​​​​​​​, last access: 4 March 2021), which is based on the Tikhonov regularization method (Legchenko and Shushakov, 1998). Samovar assumes the default calibration coefficient C_p of 7E–09 for sandy aquifers and aquifers composed of weathered and highly fractured rock based on MRS calibration experience in France (Legchenko et al., 2004). In this study, C_p was estimated using only pumping test data.

3.5.2 Magnetic susceptibility

The magnetic susceptibility of rocks changes the local geomagnetic field. The magnetic rocks, which lead to different gradients and intensities of the geomagnetic field, result in different Larmor frequencies and further can make the MRS signal undetectable (Lubczynski and Roy, 2007; Plata and Rubio, 2007). The MRS sounding is usually not possible when the magnetic susceptibility is larger than 10⁻² SI units but possible when it is lower than 10⁻³ SI units. The sounding may be or may not be possible within the interval between 10⁻² SI and 10⁻³ SI units, probably depending on the remanent magnetization of the material (Bernard, 2007). Therefore, it is always recommended to measure the magnetic susceptibility before embarking on an MRS survey (Roy et al., 2008). In this study, a portable magnetic susceptibility meter SM–20 was used to measure the magnetic susceptibility at 11 sites in the field (Fig. 5). At each site, an average magnetic susceptibility was obtained from three to five repeated measurements.

3.5.3 ERT

Subsurface resistivity depends on many different parameters, e.g. lithology, water content, and water conductivity. Its distribution in the subsurface can be visualized by 2D ERT. ERT was employed in this study because it not only provides a general understanding of the aquifer but also supports the analysis of MRS measurements.

We performed seven ERT surveys with the WGMD-9 ERT instrument manufactured by Chongqing Benteng Digital Control Technical Institute (http://www.cqbtsk.com.cn, last access: 4 March 2021​​​​​​​) in China using two standard configurations, Wenner and dipole–dipole. Wenner usually has a good signal-to-noise ratio ($S/N$) and is good at detecting vertical changes in resistivity, i.e. suitable for image horizontal structures. Dipole–dipole is sensitive to horizontal changes in resistivity, so it is ideal for vertical structure delineation. Multicore cables with a fixed electrode spacing of 10 m were used in the field. The length of cable was 890 m for ERT1–ERT4 and 810 m for ERT5–ERT7 (Fig. 5). Electrode positions were measured with a UniStrong MG858s handheld GPS (global positioning system) instrument (http://www.unistrong.com, last access: 4 March 2021), with a horizontal and vertical accuracy of 30 cm. The industry-standard RES2DINV (Loke, 1999) was employed for ERT inversion.

3.5.4 TEM

Compared to ERT, TEM also provides subsurface resistivity but with deeper penetrations, a relatively lower resolution, and a shorter time of data acquisition. The TEM instrument is usually operated in a 1D sounding mode as compared to the ERT 2D profiling mode. Since magnetic fields propagate faster in resistive media than in conductive ones, TEM is advantageous in low-resistivity media and in mapping deep conductive targets. Similar to MRS but with different constraints, there is a dead time between the excitation or transmitter function and the detection or receiver function, which are time-shared. Such TEM dead time is much shorter than in the case of MRS. TEM commonly involves placing a square loop on the targeted place and performing soundings. It generates a primary magnetic field that is abruptly interrupted to produce induced eddy currents in the subsurface. The eddy currents will lead to a secondary magnetic field, which can be detected by the loop on the ground surface. The received signals can be used to estimate subsurface resistivities by using appropriate inversion techniques (Nabighian and Macnae, 1991).

The TEM soundings were performed at 10 locations (Fig. 5) using the TEM-FAST 48 TEM instrument developed by Applied Electromagnetic Research Limited (http://www.aemr.net, last access: 4 March 2021). The TEM-FAST 48 is very small, compact, portable, and easy to deploy and apply in the field (Gonçalves, 2012). Only one TEM configuration was used, i.e. coincident square loop, consisting of one loop that combines functions of the transmitter and receiver. At each location, different loop sizes (3–95 m), time ranges (3–9 µs), stacks (5–10), and currents (0.7–1.1 A) were applied to select the optimal dataset to have the maximum investigation depth. If abrupt changes occurred in the obtained curve presenting the relation between apparent specific resistivity and time, the measurement was repeated to ensure data quality. After field collection, data were processed using TEM-Researcher proprietary software (http://www.aemr.net, last access: 4 March 2021) based on the solution of the inverse problem in time domain electromagnetic sounding.

4 Results and discussion

4.1 Borehole core lithology

The core lithology of the borehole ITC_Maqu_1 is shown in Fig. 6a, b, and c. The top layer is aeolian sand and loam. There are dunes that have been blown out of the river bed on top of the terraces. The deep layer is fluvial sediment. Based on the lithology information, the range of lithology-related parameters can be estimated. According to Chen et al. (1999), the Ruoergai basin was occupied by a large inland lake during the Quaternary before around 40 ka BP, while currently, it is a dry lake basin, with lake deposits exceeding 300 m in thickness. The extent of the ancient lake and Quaternary lake deposits is shown in Fig. 6c. Based on Fig. 6c and the log of the ITC_Maqu_1 borehole shown in Fig. 6a, it is assumed that deep sections of the Quaternary lacustrine deposits in the eastern part of the Maqu catchment are covered with the alluvial Yellow River deposits with a thickness of more than 32 m. This conclusion is generally consistent with the log of two other boreholes located to the east of the study area in Ruoergai basin, RM (33^∘57^′, 102^∘21^′) and RH (33^∘54^′, 102^∘33^′) (Fig. 6). RH is about 40 km east of the study area, with a depth of 120 m, not reaching bedrock. The top 12.4 m of coarse alluvial sediments, i.e. sands, was deposited by rivers, while the deeper deposits are lake sediments, mainly composed of silt clay, clay silt, and clay (Wang et al., 1995). RM is about 20 km east of the study area, with a depth of 310 m also not reaching bedrock. Like RH, RM core also reveals thick lake sediments, with 5.5 m river deposits on the top (Xue et al., 1998). Besides, the river deposits and alluvial sediments in RM and RH, located near a tributary of the Yellow River, are thinner than the fluvial sediment in borehole ITC_Maqu_1, which is located near the main stream of the Yellow River.

Figure 6 Borehole information: (a) the core lithology of borehole ITC_Maqu_1, (b) a picture of the core sediment when the borehole was drilled, (c) location of boreholes RM and RH (after Chen et al., 1999).

4.2 Soil thickness measurement

Results of soil thickness measurements indicate that in the mountainous west, feldspathic quartzose sandstone and sandy slate parent materials show variable soil depths related to landscape position, i.e. the soil thicknesses decrease as the slope increases, and are within 1.2 m in most cases (Fig. 7a). Based on the measurements, the relationship between the soil thickness and slope can be expressed using the equation

$$\begin{array}{}\text{(9)}& y=-\mathrm{1.1739}x+\mathrm{82}\left(\mathrm{0}\le x\le \mathrm{46}\right),\end{array}$$

where x is the slope (^∘), and y is soil thickness (cm). Equation (9) is a regression line from data obtained over residual soils in the west. The measured thickness is a result of in situ soil-forming processes, while in the east, a transported soil is observed, the thickness of which is controlled by different processes from those acting on residual soils. In general, assuming similar geology and except for the valley bottom, Eq. (9) would apply to the western study area (Fig. 7b).

Figure 7 Soil thickness: (a) soil thickness (cm) vs. slope (^∘), (b) estimated soil thickness using Eq. (9).

In the west, under the soil layer, a less weathered layer exists where water can also flow and needs to be taken into account in the conceptual model. In the field, the difference between the less weathered layer and the soil layer is that the less weathered layer contains partially weathered stones. According to the owners of three boreholes located in or near the valley (numbered 32–34 in Fig. 8), their depths are larger than 10 m and do not reach bedrock. By subtracting the estimated soil thickness (Fig. 7b) from available depth-to-bedrock estimates, for example from Yan et al. (2020) and Shangguan et al. (2017), the thickness of the less weathered layer can be estimated (Fig. 8). In the mountainous west, because the estimated depth to bedrock is often at least an order of magnitude larger than the soil thickness, the uncertainty in the less weathered layer thickness mainly depends on the uncertainty in the estimated depth to bedrock, which is high due to the lack of boreholes for appropriate training (Shangguan et al., 2017; Yan et al., 2020).

Figure 8 The estimated thickness of the less weathered layer in the west: (a) based on the ensemble model estimated depth to bedrock from Yan et al. (2020), (b) based on the depth to bedrock from Shangguan et al. (2017).

According to the study from Yan et al. (2020), the ensemble-model-estimated depth to bedrock in the mountainous west ranges from 6.2 to 70.0 m, and in the east it is generally within the range of 70 to 95 m, while according to Shangguan et al. (2017), the depth to bedrock in the mountainous west ranges from 6.4 to 35.0 m and in the east is generally within the range of 35 to 74 m. Both studies appear to underestimate depth to bedrock in the east, which is also due to the lack of boreholes for proper constraining. Though with high uncertainty, the calculated less weathered layer thickness in the west can be applied for characterizing the hydrogeological conceptual model since the west mainly plays a role in collecting water, and there are no other data available to determine the thickness of the less weathered layer.

4.3 Hydrogeological surveys

4.3.1 Water table depth measurement

For the altitude survey, ALOS RT1, with a spatial resolution of 12.5 m, performed better than other DEM products across the whole study area and had a higher resolution than the others. It was the most suitable DEM to be used in this study area for determining water table (WT) depths. For details, see Appendix A1.

There were 22 WT depths measured in 2018 and 18 in 2019 (Fig. 3). In the flat eastern area, the WT depths were interpolated with the software Surfer (https://www.goldensoftware.com/products/surfer, last access: 4 March 2021​) using the default ordinary kriging method with the linear variogram model (slope = 1, anisotropy ratio = 1, anisotropy angle = 0) (Cressie, 1990, 1991), which provides reasonable grids in most circumstances. Owing to the fact that most people living in the mountainous west use water from streams (via field survey), the need for groundwater is low, and only few boreholes exist. As such, only three boreholes numbered 32–34 were found in that western area (Fig. 9), and WT depths were measured. Normally, a good interpolation of WT depths or piezometric head over a large area needs and uses every measurement. But in this case, a reasonable WT depth map or piezometric head map in the mountainous west will need more than 100 borehole measurements (Hopkins and Anderson, 2016) because the ground surface elevation changes dramatically in the west, and so does the groundwater level. The three boreholes are far from enough to provide a reasonable WT depth map or piezometric head map and, therefore, were excluded from the interpolation. In contrast, the measured groundwater depths (and the interpolation) in the eastern study area can give a reasonable WT depth map or piezometric head map (Fig. 9a and b). Among the 34 boreholes with WT depth measurements, ground surface elevations measured by GPS-RTK were only available for 13 boreholes. For the remaining 21 boreholes, ground surface elevations were extracted from ALOS RT1. These elevations were used to derive hydraulic heads by subtracting the WT depths from the ground surface elevation. Similarly, using the kriging method, hydraulic heads were interpolated to obtain piezometric maps in the flat east (Fig. 9c and d).

Figure 9 Water table depths (m) and piezometric heads (m a.s.l.) of the eastern Maqu catchment. Panels (a) and (b) are water table depths (m) of the eastern Maqu catchment in 2018 and 2019, respectively; (c) and (d) are piezometric heads (m a.s.l.) of the eastern Maqu catchment in 2018 and 2019, respectively; (e) is the difference (m) in water table depth between 2018 and 2019; (f) is the profile 1 and profile 2 showing the hydraulic head and ground surface, with the profile locations shown in Fig. 9c and d. Numbers from 1 to 34 are identification numbers of boreholes.

In general, the range of WT depth was between 0.0 and 19.1 m in 2018 and between 0.7 and 18.0 m in 2019. In both 2018 and 2019, the interpolated WT depths (Fig. 9a and b) show a similar trend; i.e. the depth increases from the middle of the study area to the eastern boundary. The difference in WT depth in 2018 and 2019 (Fig. 9e) is probably caused by (1) different positions and number of control points, (2) the gates being open to reduce water storage in the reservoir (Fig. 9e) in 2019 to facilitate nearby constructions, (3) the interannual variation in precipitation and evapotranspiration. Nevertheless, in both 2018 and 2019, hydraulic heads (Fig. 9c and d) decrease from the middle of the study area to the eastern boundary, meaning that the groundwater flow is from the west to the east, with a hydraulic gradient of about 0.002 (dimensionless), recharging the Yellow River (Fig. 9f). This is consistent with the conclusion from Chang (2009). Ground surface elevations in Fig. 9f were extracted from ALOS RT1, and hydraulic heads were extracted from Fig. 9c and d. Some hydraulic heads are higher than the ground surface elevations as shown in Fig. 9f, which is due to (1) the accuracy of ALOS RT1 and (2) the lack of control points of hydraulic heads.

4.3.2 Aquifer tests

There were 11 aquifer tests conducted (Figs. 10 and 3) in unconfined aquifers, in partially penetrating boreholes. Eight slug tests were done in boreholes numbered 16, 21, 24, 26, 27, 32, 33, and 34 (Fig. 9); two pumping tests were carried out in 2019 at boreholes 6 and 23 (Fig. 9); and one pumping test was carried out in 2017 at the borehole ITC_Maqu_1. Data were processed automatically in AquiferTest software (https://www.waterloohydrogeologic.com/products/aquifertest/, last access: 4 March 2021) with assumptions made considering the average conditions in the study area: the aquifer is unconfined and 35 m thick; the borehole is partially penetrating; the screen radius is 0.27 m; the screen length is 15 m; the distance from aquifer top to screen bottom is 15 m; the casing radius is 0.25 m; the borehole radius is 0.3 m. Considering the spatial heterogeneity of hydraulic conductivities, they were not interpolated in the study area. As a result, the hydraulic conductivities ranged from 0.1 to 15.6 m d⁻¹ (Fig. 10). According to Healy et al. (2007), this range of hydraulic conductivity can be classified as fine silty sand to coarse clean sand. The slug test in borehole 24 provided lower estimates of hydraulic conductivity than the nearby pumping test in borehole ITC_Maqu_1, which likely underestimated the hydraulic conductivity as borehole 24 was not used for a period of time. Therefore, compared to the slug test, the hydraulic conductivity obtained from the pumping test is more accurate and is a volumetric average, which makes it more suitable to calibrate C_p because MRS results are also volumetric averages.

Figure 10 Hydraulic conductivity (m d⁻¹) obtained from aquifer tests in the Maqu catchment.

4.4 Hydrogeophysical surveys

4.4.1 MRS

ERT results were used to establish geoelectrical models for MRS inversion (see Appendix A3). The magnetic susceptibility measurements reveal very low susceptibility in the catchment, ensuring the suitability of applying MRS in the study area (see Appendix A2).

We did not use default inversion parameters because they sometimes result in abrupt changes or discontinuities of water content at two adjacent MRS sounding sites and lead to unrealistic results (Shah et al., 2008). Some excitations were excluded during inversion based on $S/N$ and the mismatch in terms of amplitude, Larmor frequency, and phase. The inversion parameters are listed in Table A7. The temperature of the water leads to different water densities and viscosities and therefore also influences hydraulic parameters. In Samovar V6.6, a default temperature of 20 ^∘C is used. But in the study area, the average groundwater temperature is 6.2 ^∘C. Therefore, it was necessary to take the true groundwater temperature into account when estimating hydraulic parameters. Thus, by correcting gravitational acceleration, water densities, and water viscosities using Eq. (10), a correction factor of 0.69 was used during the inversion process to improve accuracy.

$$\begin{array}{}\text{(10)}& K=k\mathit{\rho }g/\mathit{\eta },\end{array}$$

where K is hydraulic conductivity (m d⁻¹), k is the permeability of porous media (m²), ρ is water density (kg m⁻³), g is the gravitational acceleration (m s⁻²), and η is water viscosity (Pa s​​​​​​​).

MRS3-1 sounding (Fig. 5) was used to calculate the calibration coefficient C_p because it is the nearest MRS sounding to the borehole ITC_Maqu_1 for which pumping test data are available. Using a single point of calibration, the calibration coefficient C_p can be estimated with the uncertainty ≤150 % (Boucher et al., 2009). The calibrated C_p is $\mathrm{8.78}\times {\mathrm{10}}^{-\mathrm{9}}$ for T₁ and $\mathrm{8.13}\times {\mathrm{10}}^{-\mathrm{9}}$ ​​​ for $T_{\mathrm{2}}^{\ast }$. Figure 11 shows inversion results of water content and T₁ of MRS2-1, MRS3-1, and MRS3-2, and complete results are shown in Fig. A3 in the Appendix. Details are available in MRS inversion results.xlsx at the National Tibetan Plateau Data Center, including T₁, $T_{\mathrm{2}}^{\ast }$, water content, T₁- and $T_{\mathrm{2}}^{\ast }$-derived hydraulic conductivities $K_{T_{\mathrm{1}}}$ and $K_{T_{\mathrm{2}}^{\ast }}$, and T₁- and $T_{\mathrm{2}}^{\ast }$-derived transmissivities $T_{T_{\mathrm{1}}}$ and $T_{T_{\mathrm{2}}^{\ast }}$.

Figure 11 Water content and T₁ of MRS2-1, MRS3-1, and MRS3-2.

The water content distribution of MRS9-2, MRS7-2, MRS7-1, and MRS4-1 (Fig. A3) extends down to 150 m deep. Except for MRS4-1, soundings MRS9-2, MRS7-2, and MRS7-1 are adjacent, indicating that in the south-east, near the Yellow River, the groundwater extends to more than 150 m depth. So it is concluded that the flat east plays the main role in storing groundwater, and the groundwater can extend to more than 150 m depth.

Limiting values of 0.00 and 1000.00 ms for $T_{\mathrm{2}}^{\ast }$ and 0.00 and 3000.00 ms for T₁ are indicators that a valid numerical solution to the measured records (i.e. the inversion) was not reached, and no valid outcome is available. Except for invalid values, T₁-derived hydraulic conductivity ($K_{T_{\mathrm{1}}}$) ranges from 0.00 to 210.98 m d⁻¹, and $T_{\mathrm{2}}^{\ast }$-derived hydraulic conductivity ($K_{T_{\mathrm{2}}^{\ast }}$) ranges from 0.00 to 19.64 m d⁻¹. The value of 0.00 m d⁻¹ comes from the estimation of very low water content. Here, an order-of-magnitude difference is observed between the range of $K_{T_{\mathrm{1}}}$ and $K_{T_{\mathrm{2}}^{\ast }}$, which is due to the big difference between T₁ and $T_{\mathrm{2}}^{\ast }$. In theory, T₁ is less affected by magnetic heterogeneities and thus permits a better estimation of the hydraulic conductivity compared to $T_{\mathrm{2}}^{\ast }$. However, it is to note that no magnetic disturbance is expected in the Maqu catchment (Fig. A1). Furthermore in the case of T₁, because two timed delayed responses are compounded, any model mismatch, e.g. the MRS loop sampled volume being significantly different from a layered model parallel to the loop due to near-source river deposit media heterogeneity, can make the measured responses “doubly” distorted and may not fit a T₁ expected response. In both cases, T₁ and $T_{\mathrm{2}}^{\ast }$, a distortion is occurring. Nevertheless, according to specific circumstances, $T_{\mathrm{2}}^{\ast }$, which is evaluated from rest with a single pulse, may undergo less severe overall distortion. So $K_{T_{\mathrm{2}}^{\ast }}$ and $T_{T_{\mathrm{2}}^{\ast }}$ tend to be more reliable than $K_{T_{\mathrm{1}}}$ and $T_{T_{\mathrm{1}}}$ and should be used for future study. By checking the values of $K_{T_{\mathrm{2}}^{\ast }}$, it is concluded that there is an unconfined aquifer in the eastern study area. Based on $K_{T_{\mathrm{1}}}$ (and water content results), with a proper threshold to define aquifer and non-aquifer, the aquifer geometry can be defined.

MRS has its own limitations in that some of the in situ water information is missing and that the current “window of the technique” is only sensitive to the larger pore fraction of water content. The near-source river environment leads to the unknown mixture of varied lithology. Missing water is unknown, but accounting for a variety of lithology, including fine-pore lithology, from water table depth (Fig. 9) to the base of the aquifer (50 to 208 m range; see the following TEM section) may lead to well over 50 % missed water (Boucher et al., 2011). Besides, the invalid values for $T_{\mathrm{2}}^{\ast }$ and T₁ may be attributed to the hydrogeological conditions, such as highly heterogeneous lithology or too-low signal-to-noise ratio, and may be eased using an updated version of Samovar V6.6, such as Samovar V11.4, which not only improves the capability of signal analysis, for example, and allows the optimization of the number of inverted layers but also adds an uncertainty estimation function by incorporating singular value decomposition. Nevertheless, in highly heterogeneous environments, the indetermination of some parameters may remain with current technology. In terms of using default inversion parameters, part of the difficulty is in fitting the observed data to a too-large number of layers: i.e. partly fitting to the noise component of the records. The heterogeneity of the near-source river environment is also contributing to this difficulty. With more recent tools, like Samovar V11.4, the difficulty can be better handled (Legchenko et al., 2017).

In general, MRS provides preliminary and valuable information on water content, hydraulic conductivity, and transmissivity. Once the geometrical mapping and its fabric have been mapped, groundwater flow parameters and groundwater storage or volume can be better determined.

4.4.2 TEM

A variety of acquisition parameters were applied, and the signals were analysed to determine the optimal TEM measurements for each site. Detailed information of 10 optimal TEM measurements and inversion parameters is listed in Tables 3 and A8, respectively. The industrial noise filter was set at 50 Hz, and the amplifier was off. In the study area, using the square loop with a side length of 48 m or 95 m, the maximum time of 1 ms or 4 ms, stack between 5–10, and current of 0.8 A or 1.1 A, the TEM method can reach the maximum investigation depth ranging from 150 to more than 1000 m. For data processing, the invalid data points were first removed, and the field data were smoothed. Then the voltage data were inverted to obtain an apparent resistivity, which was further inverted to get true resistivity (Abiye and Haile, 2008; Cosentino et al., 2007). Applying a variety of acquisition parameters ensures that bad data can be revealed (Auken et al., 2006). Induced polarization (IP) and superparamagnetic (SPM) effects were not considered in the inversion process (Macnae, 2016). Because of the dead time and the fact that at most sites, a relatively dry layer of sediments exists near the ground surface with a corresponding high-resistivity depth interval, the upper 15 to 30 m of the sounding is lost, although subsequent layered Earth modelling attempts fill the gap.

Table 3 Acquisition parameters of optimal TEM data.

The root mean square (RMS) error in the inversion results shown in Fig. 12 is below 2 % in the flat area and below 10 % in the mountainous area. The results in the mountainous area, i.e. results of TEM6, TEM7, and TEM8, indicate that the resistivity becomes larger in the deep subsurface and is consistent with our understanding that the bedrock is located at relatively shallow depth from the ground surface. The maximum investigation depth of TEM6 is shallow; only 10 time windows were available and resulted in about 150 m investigation depth from the ground surface. This may be due to the local unknown geological condition. In addition to consolidated rock resistivity of the order of 2 to 4 kΩ m, TEM7 and TEM8 responses may show instances of fracturing, weathering, or faulting. To determine exactly what structure it is and the scope of the structure, further investigation is needed, for example, a systematic high-spatial-resolution geophysical survey with appropriate depth capability, such as the airborne electromagnetic survey, followed by systematic borehole drilling. The rest of the TEM measurements are scattered in the east, where it is likely that lake deposits are covered by river deposits on the top. Because the clay silt lithology has a lower resistivity than sand-rich lithology, and Chen et al. (1999) suggested that the ancient lake in Ruoergai basin was a freshwater or slightly saline lake for most of its life, the decrease in resistivity may indicate the change from river deposits to lake deposits. Table 4 listed the TEM-derived depth of the bottom of river deposits in the east. For TEM0, TEM1, TEM2, TEM3, TEM4, and TEM9, the bottoms of river deposits are deeper than 100 m, with lake deposits underneath. But for TEM5, located near the mountainous area, the bottom of river deposits is at 50 m deep, followed by 64 m thick lake deposits, with the bedrock at the bottom, and the nearest MRS sounding MRS6-1 indeed shows that there is no free water under 50 m depth.

Figure 12 Apparent resistivity with depth. The red triangles connected by the red line represent the measured apparent resistivity values, and the red line without triangles represents the inverted 1-D geoelectric model.

Table 4 TEM-derived depth of river deposits bottom in the east.

The large contact area between river and lake deposits in the eastern study area and the Ruoergai basin allows large water exchange. According to Fig. 12, the thickness of the lake deposits can reach more than 150 m. The thick lake deposits constitute a nearly infinite hydraulic buffer reservoir, allowing multi-year storage equivalent in the Maqu catchment. Thus it would be relevant to include the deep lake deposits in future hydrogeological models.

4.5 Uncertainties

As shown in Fig. 2, direct techniques, i.e. particle size analysis, altitude survey, soil thickness measurement, water table depth measurement, aquifer test, and magnetic susceptibility measurements, have low uncertainties. There are random errors for particle size analysis (Wang, 2011), but they are small and not expected to affect the final lithology result (ASTM D6913/D6913M-17, 2017) and thus can be neglected. For measured ground surface elevations, soil thicknesses, and water table depths, the uncertainty is supposed to be within a few centimetres considering the accuracy of equipment and errors during the measurement process (Burt, 2014; Cunningham and Schalk, 2011; Rydlund Jr. and Densmore, 2012). In terms of hydraulic heads derived from ALOS RT1 in boreholes, the uncertainty comes not only from water table depth measurement but also from ALOS RT1, which contains the mean absolute error of 4.4 m in the study area based on our results (Table A1). For hydraulic conductivities obtained from aquifer tests, the uncertainty mainly comes from data collection and processing. Though the duration of pumping in the borehole ITC_Maqu_1 did not reach 48 h, the water levels became steady very soon after the pumping started, so the uncertainty is estimated to be within 25 % according to studies from Brown et al. (1995) and Delnaz et al. (2019). For magnetic susceptibility, although the resolution of SM–20 is $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{6}}$ SI units, the actual reading accuracy is dependant on appropriate corrections, e.g. temperature, shape, volume, effective distance to sensor, etc. The corrections may reach a few orders of magnitude for volume and up to $\pm \sim \mathrm{50}$ % for shape (Hoffman, 2006). In the case of the Maqu catchment, these are so far from the levels at which the MRS problem may occur that, corrected or not, the final results will still be below the threshold for concerns.

In terms of indirect techniques, ERT, MRS, and TEM, performances of the raw data could be evaluated with parameters such as $S/N$ for MRS and number of bad data for ERT and TEM. Knowledge of the subsurface geometry and fabric would lead to the resolution of the main uncertainty issues for inverted data because there are implicit modelling assumptions for each method. For example, the assumption for MRS is that the subsurface is made of 1D planar layers parallel to the MRS loop with depth-increasing thickness. We cannot quantify the extent to which these assumptions are met and therefore also the extent to which the inversion data are accurate measurements of the site's hydrogeological parameters; thus appropriate uncertainty figures cannot be reliably generated for inverted data. The inverted data, as an illustration of what can be extracted from the raw data, are preliminary results with only inversion RMS errors quantified (ERT and TEM). Lake deposits, being far from the source, should not suffer from the near-source river deposits' heterogeneity, but their lithology makes it insensitive to MRS.

To further reduce uncertainties from indirect techniques, ERT, MRS, and TEM, it is important to determine subsurface geometry and its fabric. State-of-the-art airborne electromagnetic technology allows high-spatial-resolution mapping down to 500 m depth and is probably the most appropriate tool for now (Legault, 2015). After the site geometry is properly mapped, and the subsurface fabric is properly understood, optimum borehole drilling locations can be selected. When the detailed geometrical mapping of the subsurface and systematic borehole information are available, the inversion process can be better constrained and improved (Galazoulas et al., 2015; Vouillamoz et al., 2005; Wang et al., 2021).

5 Data availability

The raw dataset is archived and freely available in the DANS repository under the link https://doi.org/10.17026/dans-z6t-zpn7 (Li et al., 2020a), and the dataset containing the processed ERT, MRS, and TEM data is available at the National Tibetan Plateau Data Center with the link https://doi.org/10.11888/Hydro.tpdc.271221 (Li et al., 2020b).

6 Conclusion

We conducted borehole core lithology analysis, an altitude survey, soil thickness measurement, hydrogeological surveys, and hydrogeophysical surveys in the Maqu catchment of the Yellow River source region on the Tibetan Plateau, where little subsurface data are available. Seven DEMs were evaluated using GPS-RTK-measured ground surface elevations, and ALOS RT1 is suggested to be used in this study because of its better accuracy as well as higher resolution than other DEMs. The medium-deep lithology of the subsurface down to 32 m below the ground surface, mainly composed of sand, is available from one borehole (ITC_Maqu_1). It helps provide useful information for the cross-validation with hydrogeophysical survey results. According to the in situ measurements, soil thicknesses, except on the valley floor, are within 1.2 m depth in most cases in the west, and the soil thickness decreases as the slope increases. The hydrogeological surveys reveal that groundwater flows from the west to the east, recharging the Yellow River, and the hydraulic conductivity ranges from 0.2 to 12.4 m d⁻¹. The hydrogeophysical surveys demonstrate the presence of an unconfined aquifer in the east. The water content and hydraulic parameters of that aquifer were estimated at MRS sounding locations. The depth of the Yellow River deposits was derived at TEM sounding positions in the flat eastern area.

Although water table depths were only measured once or twice, and hydrogeophysical methods, like ERT, TEM, and MRS, have inherent non-uniqueness problems during the inversion process, they all provide valuable information, especially in such data-scarce areas as the TP. The data from direct techniques with low uncertainties, as shown in Fig. 2 (lithology, ground surface elevation, soil thickness, water table depth measurement, hydraulic conductivity from aquifer test, magnetic susceptibility), can contribute to related global or regional databases where the in situ data over the TP are scarce or be regarded as verification and validation data for groundwater modelling over the Maqu catchment. Indirect techniques such as ERT, MRS, and TEM are a rare, unique, and particularly rich training data source for geoscientists interested in the data processing and interpretation of the particular hydrogeological and hydrogeophysical techniques used here. It is a dynamic set where additional complementary data will gradually add constraints to the inversion processes. For example, a researcher developing new techniques for $S/N$ improvements of some of these techniques will get free and highly relevant data to work with.

In the Maqu catchment, a near-source river environment without adequate geometrical mapping, the representativity of the various sampled volumes is unknown as well as whether the sampled volume fits the models used for data inversion. This is much less the case farther away from the source where homogeneity and fitting of the model to the actual hydrogeological setting are achieved. In such an away-from-source case, pumping test data may be assumed to be representative of the tested formation, while in techniques such as MRS, depth and thickness information may be extracted from the datasets as well as hydraulic estimates. This is an ongoing project, and it may become available later if such above-mentioned mapping is completed. Any further similar surveys and borehole drillings would benefit from such geometrical mapping since their precise localization may then be optimized in view of proper data inversion and information gap filling.

Generally speaking, all the methods work well in the study area and have confirmed the presence of an unconfined fluvial aquifer within the 250 m below the surface and the presence of lake deposits with much finer pore lithology. By combining our dataset with available depth-to-bedrock datasets, a preliminary hydrogeological conceptual model can be established. If combining our dataset with detailed geometrical mapping of the subsurface and deep borehole information, a more complete and accurate conceptual model can be obtained. Furthermore, we will be monitoring the groundwater and surface water in the study area and aim for establishing a long-term monitoring network, which will eventually contribute to the verification and validation of future studies on groundwater modelling over the Maqu catchment. To our knowledge, this is the first time such detailed surveys were conducted in a TP catchment in order to set up a hydrogeological conceptual and numerical groundwater model. This paper is expected not only to contribute to the hydrogeological conceptual and numerical model of the Maqu catchment at the TP but also to provide data for hydrogeological and hydrogeophysical communities and promote interdisciplinary research.

Appendix A

A1 Altitude survey

A total of 46 ground surface elevations were measured (33 in the flat east, 13 in the mountainous west) and were used to evaluate the accuracies of seven DEM datasets (data available in DEM.xlsx in the National Tibetan Plateau Data Center), and the most accurate one was applied in this study. The statistical analysis results of the seven DEMs in the study area are shown in Table A1. The root mean square errors (RMSEs) of ALOS RT1 and ALOS RT2 are 5.695 and 5.477 m, respectively, much smaller than the RMSE of the other five DEMs. The correlation coefficient, the mean error, and the mean absolute errors in ALOS RT1 and ALOS RT2 also show better performance than those of the other five DEMs. Comparing ALOS RT1 with ALOS RT2, ALOS RT1 slightly outperforms ALOS RT2 with regards to RMSE, correlation coefficient, and the mean error. Tables A2 and A3 list the statistical analysis results of seven DEMs separately for the flat eastern area and the mountainous western area. All seven DEMs behave better in the west than the east in terms of the correlation coefficient. In the west, the correlation coefficients of seven DEMs are all larger than 0.94, while in the east, the correlation coefficients are all lower than 0.24. This is because the range of ground surface elevation in the flat east is much smaller than the range of elevation in the mountainous west. With regard to the RMSE, mean error, and mean absolute error, all seven DEMs have better behaviour in the east than in the west. In general, ALOS RT1 and ALOS RT2 also outperform the other five DEMs, according to Tables A2 and A3.

Since ALOS RT1 performs slightly better than ALOS RT2 in the whole study area and has a higher resolution than ALOS RT2, it is the most suitable DEM to use in this study area. For ALOS RT1 in the flat east, 52 % of errors (DEM value − GPS-RTK value) are within the range of −3 to 3 m, and 79 % of errors are within the scope of −5 to 5 m, while in the mountainous west, 54 % of errors are within the range of −8 to −12 m, and 46 % of errors are within the range of 0 to 7 m.

Previous TP works about DEM evaluation mainly focused on SRTM and ASTER. Our results are generally consistent with previous studies in terms of RMSE of SRTM. Nan et al. (2015) evaluated the height accuracy of SRTM and ASTER in the eastern TP with reference to the relatively high precision of the 1:50 000 scale DEM surveyed and mapped by the State Bureau of Surveying and Mapping in China. As a result, the RMSEs of SRTM and ASTER are 35.3 and 50.2 m, respectively. Ye et al. (2011) evaluated SRTM and ASTER in the Mt. Qomolangma (Mt. Everest) area on the TP by comparing 211 elevation checkpoints on the 1:50 000 topographic maps surveyed and mapped by the State Bureau of Surveying and Mapping in China, demonstrating an average height difference of 31.3 and 44.9 m for SRTM and ASTER, respectively. However, there are other studies that have different evaluation results. Fujita et al. (2008) found that the elevation differences between DEMs and ground survey data from differential GPS were 11.0 m for ASTER and 11.3 m for SRTM in the Lunana region, Bhutan Himalaya. The DEM evaluation results also indicated that in different places over the TP, the satellite DEM estimates are acquired with varying accuracy. This may be due to different topographic complexity in different areas.

The DEMs' quality can be influenced by several factors, such as sensor type, algorithm, terrain type, and grid spacing (Hebeler and Purves, 2009). In this study, grid spacings of DEMs are similar except for ALOS RT1, so the main factors that affect the accuracy of the DEMs should be sensor types and algorithms. For SRTM, the issue inherent to the production method is mast oscillations, while for ASTER and AW3D30, the issue is scene mismatch (Grohmann, 2018). As for radiometrically terrain-corrected (RTC) products ALOS RT1 and ALOS RT2, the quality is directly related to the quality of the source DEM SRTM which was used in the RTC process. This results in very similar correlation coefficients of SRTM, ALOS RT1, and ALOS RT2 and obvious improvements in RMSE, MAE, and ME (Table A1).

Table A1 Statistical analysis of seven DEMs in the study area.

^* Error = DEM value−GPS-RTK value.

Table A2 Statistical analysis of seven DEMs in the flat eastern area.

^* Error = DEM value−GPS-RTK value.

Table A3 Statistical analysis of seven DEMs in the mountainous western area.

^* Error = DEM value−GPS-RTK value.

A2 Magnetic susceptibility

The magnetic susceptibility measurements (Fig. A1) reveal very low susceptibility in the catchment, with susceptibility values all smaller than $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{5}}$ SI units, with an average of $\mathrm{3}\times {\mathrm{10}}^{-\mathrm{6}}$ SI units. A previous study from Chen et al. (1999) also reported low magnetic susceptibility of the RH core (Fig. 6c) with 120 m length located 40 km east of the study area in the Ruoergai basin. The low magnetic susceptibility ensured the suitability of applying MRS in the study area.

Figure A1 Magnetic susceptibility measurements (10⁻³ SI units) ensured the suitability of applying MRS in the study area.

A3 ERT

Detailed information of ERT profiles and inversion parameters is listed in Tables A4 and A5, respectively. For a specific pseudo-depth, the values between adjacent points generally vary smoothly. Bad data points can be easily identified as they appeared as outlier points in the pseudo-section plot in RES2DINV due to their too-high or too-low apparent resistivity values. The bad data points were filtered out based on the following criteria: (i) having negative apparent resistivity or small apparent resistivity close to 0 Ω m or (ii) having negative or positive pulse amplitude ratios <0.75 or >1.33 (a measure of waveform symmetry) (Slater et al., 2010; Wilkinson et al., 2016). After filtering, the least square method was used for the inversion. Results are shown in Fig. A2, with the root mean square (RMS) error less than 5 %. A pattern of roughly regular parallel-to-surface electrostratigraphy is observed in all ERT profiles, except 0–310 m of profile ERT5, where the pattern is dipping relative to the surface. This means that strata are likely to be stratified in most parts of the study area. For ERT2, ERT3, ERT5, and ERT6, three electrostratigraphic layers can be identified: the first layer with the highest resistivity, the second layer with the lowest resistivity, and the third layer with a medium resistivity. The second layer is likely to represent an aquifer. However, considering ERT4 and ERT7, there is a lack of marker electrostratum; i.e. a layer with high resistivity does not exist at the ground surface. This is probably due to high water content near the ground surface in the mountainous area where ERT4 and ERT7 were located. As for ERT1, rainfall occurred during the field measurement. Rainwater accumulations occurred next to some of the electrodes, causing abnormal current distribution during the ERT measurements, and about half of the data were filtered out. The ERT1 inversion results show a three-layer pattern similar to the one observed along the ERT2, ERT3, ERT5, and ERT6 profiles. One or more short-wavelength anomalies (<200 m) are observed along all profiles but particularly in the case of ERT1, ERT3, and ERT6. Short-wavelength anomalies along ERT1 may be due to data acquisition made during rainfall, while in the case of the other profiles, localized changes in water content or lithology variations are suspected.

Compared to the dipole–dipole configuration, the investigation depth of the Wenner configuration is deeper. So resistivity values obtained from the Wenner configuration were used to establish geoelectrical models for MRS inversion. For ERT2, ERT3, ERT5, and ERT6, three-layer geoelectrical models were extracted, while for ERT4 and ERT7, two-layer geoelectrical models were extracted. ERT1 was neglected due to the influence of rainfall. For ERT5, from 0 to 310 m, there is a topographic change; the ground surface elevation decreases from 3395 m a.s.l. and stabilizes at around 3390 m a.s.l. Ground surface with low resistivity exists along this 310 m transect. Because the MRS soundings were conducted in flat areas, only resistivity from 310 to 810 m was used for the first layer of the geoelectrical model. The geoelectrical models and corresponding MRS measurements are shown in Table A6. The depths of the last layer of geoelectrical models are extended to 1.5 times the MRS investigation depth since signal distortion due to subsurface resistivity is calculated down to that depth while making the MRS linear filter. In this particular version, MRS investigation depth was considered to be the MRS loop size, i.e. 150 and 100 m. Nevertheless, like other geophysical methods, ERT has equivalence problems, i.e. non-uniqueness of inversion results. Despite equivalence problems, the ERT method still provides important subsurface information in the Maqu catchment, where we have little fundamental information. This is the very first investigation in this area; when more lithology information becomes available later, ERT can be better constrained to reflect the subsurface lithology.

Table A4 Detailed information on ERT.

Table A5 Inversion parameters for ERT.

Table A6 Geoelectrical models used for MRS inversion.

^* The depth of the third layer is 150 m rather than 225 m.

Figure A2 ERT measurements and corresponding ground surface elevation and vertical exaggeration (VE).

A4 Inversion parameters for MRS

Table A7 Inversion parameters for MRS.

A5 MRS inversion results

Figure A3 Water content and T₁ derived from MRS measurements.

A6 TEM inversion parameters

Table A8 TEM inversion parameters.

A7 Abbreviations

ALOS PALSAR RT1	Advanced Land Observing Satellite – Phase Array type L-band Synthetic Aperture Radar – high terrain correction resolution
ALOS PALSAR RT2	Advanced Land Observing Satellite – Phase Array type L-band Synthetic Aperture Radar – low terrain correction resolution
AMSR-E	Advanced Microwave Scanning Radiometer for Earth Observing System
ASCAT	Advanced Scatterometer
ASF	Alaska Satellite Facility
ASTER	The Terra Advanced Spaceborne Thermal Emission and Reflection Radiometer
CAS	Chinese Academy of Science
CLM	Community Land Model
CPC	Climate Prediction Center
DEM	Digital elevation model
ERT	Electrical Resistivity Tomography
GDEM	Global digital elevation model
GLDAS	Global Land Data Assimilation System
GPS	Global positioning system
GPS-RTK	Real-time kinematic global positioning system
GRACE	Gravity Recovery and Climate Experiment
JAXA	Japan Aerospace Exploration Agency
LAPSUS	LandscApe ProcesS modeling at mUlti-dimensions and scaleS
ME	Mean error
MAE	Mean absolute error
MRI	Magnetic resonance imaging
MRS	Magnetic resonance sounding
NMR	Nuclear magnetic resonance
RMSE	Root mean square error
SRTM	Shuttle Radar Topography Mission
TDEM	Time domain electromagnetic
TEM	Transient Electromagnetic
TP	Tibetan Plateau
USGS	United States Geological Survey
VIC	Variable Infiltration Capacity model
WT	Water table
YRSR	Yellow River source region