Regional Assessment of Groundwater Recharge in the Lower Mekong Basin

Groundwater recharge remains almost totally unknown across the Mekong River Basin, hindering the evaluation of groundwater potential for irrigation. A regional regression model was developed to map groundwater recharge across the Lower Mekong Basin where agricultural water demand is increasing, especially during the dry season. The model was calibrated with baseflow computed with the local-minimum flow separation method applied to streamflow recorded in 65 unregulated sub-catchments since 1951. Our results, in agreement with previous local studies, indicate that spatial variations in groundwater recharge are predominantly controlled by the climate (rainfall and evapotranspiration) while aquifer characteristics seem to play a secondary role at this regional scale. While this analysis suggests large scope for expanding agricultural groundwater use, the map derived from this study provides a simple way to assess the limits of groundwater-fed irrigation development. Further data measurements to capture local variations in hydrogeology will be required to refine the evaluation of recharge rates to support practical implementations.


Introduction
Groundwater is a primary resource for agricultural, domestic and industrial uses globally [1].However, it remains largely underutilized in the Mekong Basin [2] and particularly in parts of Laos and Cambodia where the primary source for irrigation is surface water [3,4].Groundwater resources are better protected from pollution than surface water, and allow for more reliable abstraction over time due to smoother inter-annual and inter-seasonal variations than the surface water strongly influenced by the monsoonal climate [5].However, over abstraction of groundwater resources can lead to rapid depletion of aquifers, with deleterious effects on ecosystems and human communities [6].Total utilizable groundwater volume in an aquifer is controlled largely by the seasonal variations in groundwater recharge, the drainage from and into rivers, and pumping.A sustainable level of resource utilization implies that outflow (drainage to the rivers and pumping) does not exceed inflow (recharge) over a multi-year period [7].One of the first steps required to assess the sustainability of an aquifer for potential water uses involves the assessment of the groundwater recharge rate.At broad scale, groundwater recharge occurs either in diffuse form (from rainfall infiltration via the unsaturated soil zone), or more locally via water seepage through riverbeds and other water bodies such as lakes and reservoirs.Very few studies have attempted to quantify groundwater recharge in the Mekong Basin [3,[8][9][10].
Several methods can be used to quantify groundwater recharge using different types of data (geological, hydrological and geochemical indicators) [11] and models [12].For broad-scale assessments, river flow records are convenient because they integrate information over areas as large as the catchments of the gauged rivers, and they are a good alternative to other approaches whose data requirement is not compatible with the overall data scarcity characterizing the Mekong Basin.Flow separation methods are applied in catchments where baseflow is assumed to be equivalent to groundwater recharge over multi-year periods.This assumption, based on water budget, is valid only if there is no major pumping and reservoir (e.g., wetlands, lakes, hydropower dams) in the catchment, likely to enhance or reduce baseflow naturally originating from groundwater return flow.Flow separation methods are particularly appropriate for application in the Mekong Basin due to its humid climate and tropical hydrology [13].Humid regions are usually characterized by shallow water tables replenished at the end of the wet season mainly through diffuse recharge from rainwater infiltration, while focused recharge from water bodies plays a minor role in the overall groundwater balance [11].For catchments greater than 200 km 2 in humid tropical settings, groundwater usually discharges to rivers all year round, resulting in gaining streams, compatible with groundwater recharge estimation by stream hydrograph separation [11,14].However, baseflow separation methods inevitably involve subjective choices on the mathematical algorithms that cannot fully capture the nonlinearity of processes controlling surface-groundwater exchanges: riverbank storage, spatial variability in evaporation, recharge and storage capacity of the aquifer [15].Consequently, values of baseflow estimated with different algorithms can vary by a factor of two [5].There is no consensus on which geomorphic and land-cover characteristics are most closely linked to subsurface storage and baseflow [16].Factors that promote infiltration and recharge to subsurface storages will increase baseflow, while factors associated with higher evapotranspiration will reduce baseflow.
The objective of this study is to map diffuse groundwater recharge rates across the Lower Mekong Basin (LMB).The climatic, hydrological, and geomorphological characteristics of the LMB are described in Section 2. As detailed in Section 3, the methodological approach involves three steps: (i) quality-control and selection of suitable flow records available from 72 gauging stations along the main tributaries of the Mekong River; (ii) application of a mathematical digital filter to separate baseflow from total flow in the selected catchments where baseflow is equivalent to groundwater recharge over a multi-year period; and (iii) interpolation of groundwater recharge rates in un-gauged areas of the Lower Mekong River Basin, using multiple regression analyses of catchment geomorphologic and climatic characteristics.The results of this analysis are presented in Section 4 and discussed in Section 5.

Study Site
The Mekong is one of the world's most important trans-boundary rivers in terms of total length (4800 km), catchment area (about 795,000 km 2 ) and mean annual discharge (475 km 3 ) [13].Originating from the Tibetan Plateau in China, this river crosses Myanmar, Laos, Thailand, Cambodia and finally Vietnam where it flows into the South China Sea.The LMB includes the portion of the basin located downstream of China.It accounts for 79% of the total Mekong basin area and 84% of the flow volume [13].

Climate and Hydrology
Mean annual rainfall in the LMB ranges from less than 1000 mm on the Khorat Plateau in Thailand to more than 2500 mm on the highlands of Laos.Mean potential evapotranspiration is maximal in Southern Cambodia where it exceeds 1400 mm/year and the lowest over the mountain ranges of Laos where it falls below 1050 mm/year [17].Rainfall patterns are seasonal with about 80% of annual totals occurring during the wet season between May and October.The Mekong River discharges 75% of its annual flow between July and October which corresponds to the rainiest period of the year, caused by the tropical monsoonal climate [13].

Topography
Four major topographic regions are distinguishable in the LMB: the Northern Highlands, the Khorat Plateau, the Tonle Sap Basin and the Mekong Delta (Figure 1).The Northern Highlands with elevations more than 600 m above the mean sea level, cover Northern Thailand, Laos and the northern section of the Annamite Range in Vietnam.They are characterized by a dense forest coverage.The Mekong mainstream and its tributaries flow through steep, rock-cut valleys.The Khorat Plateau extends over Northeast Thailand and includes much of the lowland areas of Laos.It is a sandstone plateau undulating between 100 and 500 m above the mean sea level, characterized by erratic rainfall and poor coarse-textured sandy and unevenly distributed soils, often saline.The drainage network in this relatively flat landscape is mostly dendritic.Further downstream, the Tonle Sap Basin covering most of Cambodia and Southern Laos is made of hills surrounding an alluvial plain including the large freshwater lake Tonle Sap, which reverses its flow direction into the Mekong during the wet season.Elevation is usually lower than 100 m above the mean sea level.The Mekong Delta begins near Phnom Penh in Cambodia where the river becomes a delta less than 10 m above the mean sea level, expanding across Southern Vietnam and covering 62,520 km 2 of mangrove, swamps, sand dunes, spits, tidal flats, and irrigated rice paddy fields.

Topography
Four major topographic regions are distinguishable in the LMB: the Northern Highlands, the Khorat Plateau, the Tonle Sap Basin and the Mekong Delta (Figure 1).The Northern Highlands with elevations more than 600 m above the mean sea level, cover Northern Thailand, Laos and the northern section of the Annamite Range in Vietnam.They are characterized by a dense forest coverage.The Mekong mainstream and its tributaries flow through steep, rock-cut valleys.The Khorat Plateau extends over Northeast Thailand and includes much of the lowland areas of Laos.It is a sandstone plateau undulating between 100 and 500 m above the mean sea level, characterized by erratic rainfall and poor coarse-textured sandy and unevenly distributed soils, often saline.The drainage network in this relatively flat landscape is mostly dendritic.Further downstream, the Tonle Sap Basin covering most of Cambodia and Southern Laos is made of hills surrounding an alluvial plain including the large freshwater lake Tonle Sap, which reverses its flow direction into the Mekong during the wet season.Elevation is usually lower than 100 m above the mean sea level.The Mekong Delta begins near Phnom Penh in Cambodia where the river becomes a delta less than 10 m above the mean sea level, expanding across Southern Vietnam and covering 62,520 km 2 of mangrove, swamps, sand dunes, spits, tidal flats, and irrigated rice paddy fields.

Hydrogeology
Four main hydrogeological units can be delineated in the LMB (Figure 2) [3,18,19].(1) Along the eastern and southeastern border of the LMB, volcanic and granitic rocks with water-bearing features (joints, faults, and weathering zones) are overlapped by cemented early Paleozoic metasedimentary rocks with reduced porosity and permeability.However, the basalt flows of Southern Lao PDR, the Central Highlands of Vietnam, and Southeastern Cambodia have greater permeability; (2) In the Northern LMB, the porous and permeable late Paleozoic sedimentary rocks, dissected into relatively small blocks by subsequent orogeny, and topped by Mesozoic deposits, supports local groundwater

Hydrogeology
Four main hydrogeological units can be delineated in the LMB (Figure 2) [3,18,19].(1) Along the eastern and southeastern border of the LMB, volcanic and granitic rocks with water-bearing features (joints, faults, and weathering zones) are overlapped by cemented early Paleozoic metasedimentary rocks with reduced porosity and permeability.However, the basalt flows of Southern Lao PDR, the Central Highlands of Vietnam, and Southeastern Cambodia have greater permeability; (2) In the Northern LMB, the porous and permeable late Paleozoic sedimentary rocks, dissected into relatively small blocks by subsequent orogeny, and topped by Mesozoic deposits, supports local groundwater flow systems locally discharging into tributaries of the Mekong River.In some locations, the Permian-aged karstic limestone are considered to be more productive aquifers; (3) In the lowlands, and particularly in Northeast Thailand, deep confined and shallow unconfined aquifers from the Mesozoic are comprised of sandstones.Outcrops of early Paleozoic rock bounded by Upper Jurassic sandstones have lower water yields; (4) In the Mekong delta, Cenozoic alluvial and deltaic sediments of up to 800 m thick form both unconfined and confined aquifers.The alluvial deposit in the upstream narrower valley rarely exceeds 100 m in thickness.Recharge in Pleistocene through Miocene deposits occurs in terraces along the margins of the Holocene alluvial plain.These terraces are generally coarser-textured and more productive aquifers than fine-textured Holocene surficial sediments including large amount of clay.While these observations provide an overall view of the hydrogeology across the LMB, the level of understanding remains largely qualitative.The transmissivities and water storage capacities of aquifers in the LMB have yet to be measured in a systematic way [3].Only localized pumping tests allowed estimating aquifer hydraulic properties over a few areas [20,21].

Materials and Methods
After identifying suitable gauged catchments and computing daily baseflow, multiple regression analyses were applied to various climatic and geomorphologic characteristics computed for each catchment of the study sample.These variables were then used to predict baseflow values (i.e., groundwater recharge rates) in un-gauged areas of the LMB.

Catchment Selection and Baseflow Computation
Daily streamflow records acquired from the Mekong River Commission were collected from 72 gauged sub-catchments of the LMB, corresponding to 1st to 5th order tributaries of the Mekong River.Visual inspection and comparison of streamflow hydrographs allowed selecting 65 gauging stations located along 50 rivers (Figure 3), on the basis that they provide records which were not subject to dam regulation, gaps, and questionable values [22].The location and year of commission of

Materials and Methods
After identifying suitable gauged catchments and computing daily baseflow, multiple regression analyses were applied to various climatic and geomorphologic characteristics computed for each catchment of the study sample.These variables were then used to predict baseflow values (i.e., groundwater recharge rates) in un-gauged areas of the LMB.

Catchment Selection and Baseflow Computation
Daily streamflow records acquired from the Mekong River Commission were collected from 72 gauged sub-catchments of the LMB, corresponding to 1st to 5th order tributaries of the Mekong River.Visual inspection and comparison of streamflow hydrographs allowed selecting 65 gauging stations located along 50 rivers (Figure 3), on the basis that they provide records which were not subject to dam regulation, gaps, and questionable values [22].The location and year of commission of hydropower dams were inspected for that purpose [23,24].Depending on the catchments, the selected records include between 1 and 41 years of daily flow values with a median length of 17 years, circumscribed between January 1951 and December 2007.hydropower dams were inspected for that purpose [23,24].Depending on the catchments, the selected records include between 1 and 41 years of daily flow values with a median length of 17 years, circumscribed between January 1951 and December 2007.Daily baseflow at each of the 65 gauging stations was calculated using the local minimum, a filtering separation method available in the BFI+ 3.0 software package [25,26].This method has already been applied locally in the Mekong Basin to assess groundwater recharge over the basaltic plateau of Dak Lak province in Vietnam [27].This method is suitable for processing long records of discharge data at many locations because it requires only a single parameter (N) derived from the catchment area, to be determined, and can be easily translated in computer code for fast processing.Baseflow is calculated by connecting the local minimums of total daily streamflow selected in a sliding time interval moved along the flow time series by 1-day increments.The duration of this time interval is calculated as follows.The theoretical average number N of days between a flood peak and the time when surface runoff ceases, is first computed with Equation N = A 0.2 where A is the catchment area in square miles [28].The time interval includes (2N*−1)/2 days where 2N* is the odd integer between 3 and 11 nearest to 2N.In each of the 65 catchments, specific annual baseflow (mm/year) was computed for each hydrological year (1st April-31st March) available in the records and the median annual value QB estim was selected as the independent variable to perform the multiple linear regression analysis (cf.Section 3.2).
Although in some situations, the local minimum method can overestimate baseflow [29], it mostly tends to provide a lower bound of digital filtering evaluations in general [30], thus providing conservative estimations of recharge.This is convenient for two reasons: (i) given the wide range of baseflow estimations provided by different flow separation algorithms, it is informative to know that our results will be a lower bound; and (ii) risks of overestimating the available water resources with the resultant water stresses are reduced.Daily baseflow at each of the 65 gauging stations was calculated using the local minimum, a filtering separation method available in the BFI+ 3.0 software package [25,26].This method has already been applied locally in the Mekong Basin to assess groundwater recharge over the basaltic plateau of Dak Lak province in Vietnam [27].This method is suitable for processing long records of discharge data at many locations because it requires only a single parameter (N) derived from the catchment area, to be determined, and can be easily translated in computer code for fast processing.Baseflow is calculated by connecting the local minimums of total daily streamflow selected in a sliding time interval moved along the flow time series by 1-day increments.The duration of this time interval is calculated as follows.The theoretical average number N of days between a flood peak and the time when surface runoff ceases, is first computed with Equation N = A 0.2 where A is the catchment area in square miles [28].The time interval includes (2N*−1)/2 days where 2N* is the odd integer between 3 and 11 nearest to 2N.In each of the 65 catchments, specific annual baseflow (mm/year) was computed for each hydrological year (1st April-31st March) available in the records and the median annual value Q B estim was selected as the independent variable to perform the multiple linear regression analysis (cf.Section 3.2).

Prediction of Baseflow Across the Lower Mekong Basin
Although in some situations, the local minimum method can overestimate baseflow [29], it mostly tends to provide a lower bound of digital filtering evaluations in general [30], thus providing conservative estimations of recharge.This is convenient for two reasons: (i) given the wide range of baseflow estimations provided by different flow separation algorithms, it is informative to know that our results will be a lower bound; and (ii) risks of overestimating the available water resources with the resultant water stresses are reduced.

Prediction of Baseflow Across the Lower Mekong Basin
To map groundwater recharge across the LMB, multiple regressions are used to identify multi-variate statistical relationships between measured baseflows and their catchments characteristics.These relationships are then used to predict baseflow in any un-gauged sub-catchment of the LMB, based on its characteristics.Candidate catchment characteristics for inclusion in the regression equation (presented in Section 3.3) are selected based on data availability across the whole LMB, and according to their potential influence on groundwater recharge.Instead of using linear regressions, power-law equations (Equation ( 1)), frequently used in environmental predictive statistics because of their greater performance compared to linear relationships [31], are used to compute Q B predict that predicts Q B estim from m catchment characteristics X i .Their logarithmic transformation produces a linear model (Equation ( 2)) whose m + 1 coefficients β i can be determined with multiple linear regressions.
β 0 is the intercept term of the model.ν and ε are the log-normally and normally distributed errors of the models, respectively.Normality in ε distribution is usually easier to obtain than in not log-transformed linear model, hence the advantage of power-law equations.The logarithm function being defined for strictly positive values only, adding one to catchment characteristics X i including zero values allows a correct mapping between the value of ln(X i+1 ) and X i [32].The selection of the catchment characteristics X i that best predict Q B estim , and the calculation of their respective coefficients β i are performed by weighted least squares regressions applied to the 65 values of Q B estim j calculated in the 65 catchments (j = 1, . . ., 65), and their respective catchment characteristics X ij .Unlike ordinary least square regressions treating all observations Q B estim j equally, weighted least square regression enables the varying number k j of hydrological years used to calculate each value of Q B estim j to be accounted [33].Values of Q B estim j derived from a greater number of hydrological years are more precise (have lower variance) and thus should have a greater weight in the regression.However, this reliability decreases as the variance of Q B estim j increases.To account for these two counteracting factors, weights (w j ) were calculated as follows: where Stdev(Q B estim j ) is the standard deviation of Q B estim j .The explanatory variables X i (i.e., catchment characteristics) that best predict baseflow were identified among a set of 15 candidate variables (described in Section 3.3 and listed in Table 1) using the two selection algorithms "best subsets regression" and "step-wise regression" available in MiniTab 16.This selection was intended to maximize the prediction R-squared (R 2 pred ) calculated by leave-one-out cross-validations.Unlike the classical R-squared the maximization of which can lead to model over-fitting and loss of robustness, R 2  pred reflects the ability of the model to predict observations which were not used in the model calibration.Maximizing R 2  pred generally leads to greater parsimony in the number of explanatory variables.An explanatory variable was considered to be statistically significantly different from zero if its p-value, derived from the Student's t-test, was lower than 0.05.The required homoscedasticity (homogeneity of variance) of the model residuals ε was verified by visual inspection of the residual plots.Possible multi-collinearity among the explanatory variables was controlled with the variance inflation factor (VIF).The influence statistic Cooks D was used to identify and remove outlier catchments exhibiting high influence on the estimation of the model coefficients [34].The predictive power of the regression model was assessed by three criteria: the Nash-Sutcliffe efficiency coefficients (NSEC), the adjusted R-squared (R 2 adj ), and R 2 pred .While R 2 pred assesses how well Equation (2) predicts responses to new observations, R 2 adj allows comparing the performance of linear regressions including different numbers of explanatory variables.While the value of the classical R 2 systematically increases when a new explanatory variable is added in Equation ( 2), R 2 adj will increase only if the new term improves the performance of the linear regression more than what would be expected by chance alone [35].While R 2 adj and R 2 pred estimate the strength of the linear association between the estimations and predictions, NSEC measures the goodness of fit of linear and non-linear models (including power law models, i.e., Equation (1)), thus allowing performance comparison with any hydrological model.NSEC is computed as follows: where Q B estim j is the median value of specific annual baseflow computed in catchment j with the local minimum method and Q B pred j is the corresponding value predicted with the power-law model (i.e., Equation ( 1)) in the same catchment.Q B estim is the spatial mean of the estimated baseflow Q B estim j across all catchments.

Catchment Characteristics
Statistics for each catchment characteristic are listed in Table 1.

Climate
Catchment areal rainfall was computed using "Aphrodite", a 0.25 • × 0.25 • grid of daily rainfall covering monsoonal Asia over the period 1951-2007 [36].Aphrodite is one of the most reliable precipitation datasets to model discharge in the large river basins in Asia [37].Several rainfall variables were tested for correlation with Q B estim : annual and monthly rainfall, rainfall cumulated over the n-day (n = 5, 10, and 15) rainiest periods of the hydrological year.Annual median rainfall, exhibiting the greatest correlation coefficient with Q B estim , was included as the only candidate rainfall variable in the regressions (Table 1).Median annual values of catchment areal temperature and standard evapotranspiration (ET 0 ) were computed using 0.5 • × 0.5 • gridded monthly values from the Climate Research Unit [38].ET 0 was calculated with the FAO grass reference evapotranspiration Equation applied to climate variables from the same data source covering the period 1901-2009 [39].These three climate median values and Q B estim j were derived from the same k j hydrological years in each selected gauged catchment.In addition to ET 0 , median annual values of actual catchment areal evapotranspiration were computed using the land surface evapotranspiration product MODIS 16 available at daily time step with 1 km 2 resolution for the period 2000-2012 [40].

Geomorphology and Geographic Coordinates
Five geomorphological characteristics, likely to control hydrology, were derived from HydroSHEDS, a quality-controlled 90-m digital elevation model [41].These catchment characteristics, including drainage area, drainage density, mean elevation, mean slope and perimeter, were computed with ArcMap 10.0.The drainage density is the cumulative length of the stream network within the catchment, normalized by its drainage area.The stream network is made of outlet points draining an area greater than 40 km 2 .This threshold value was selected to best capture the variability of drainage densities among the studied catchments.The geographic coordinates of the catchment centroid (latitude and longitude) were selected as two additional candidate explanatory variables to capture any longitudinal and latitudinal gradients in incomputable environmental variables possibly influencing groundwater recharge across the LMB (e.g., aquifer properties).

Soil
The top-soil texture and soil depth, two soil characteristics likely to control hydrological processes and groundwater recharge were quantified using a four-level scale suggested by the Mekong River Commission (Table 2) [42].Mean values of each soil characteristic in each catchment were obtained by weighting each scale level by the respective area covered in the catchment.

Land Cover
Two land covers, likely to exert contrasted influences on hydrological processes across the LMB [43] and at the studied catchment scales [44], were selected as candidate explanatory variables: forest and rain-fed lowland rice paddy fields.Forest cover was produced by merging four forest types available as separate land-cover classes in the published map of the Mekong River Commission [42]: "coniferous forest", "deciduous forest", "evergreen forest" and "forest plantation".The land cover "rain-fed lowland rice paddy field" was directly available as an individual class in the map.While forest is usually characterized by high infiltration and evapotranspiration rates, rain-fed lowland rice paddy fields limit deep percolation due to an impermeable soil layer at the basement of the rice root zone [43].These two variables correspond to the percentage area of each land cover in each catchment.
Variables related to the hydrogeology (e.g., transmissivity, storage capacity of aquifers) could not be included in the list of candidate explanatory variables because of the absence of broad scale quantitative data in the LMB.This limitation is discussed in Section 5.

Baseflow Estimations
The median value of specific annual baseflow computed in each catchment with the local minimum method (Q B estim ) varies between 53 mm/year (catchment of the Nam Mun River at Rasi Salai station, Thailand) and about 1000 mm/year (catchment of the Nam Sane River at Muong Borikhan station, Laos) with a median value of 439 mm/year.Circles depicted in Figure 4 show the spatial distribution of these estimated recharge rates.Strikingly, extreme values are grouped regionally.Greatest values are observed in Central Laos, ranging from 625 to 1000 mm/year, and Southern Laos, ranging from 600 to 1000 mm/year.Lowest values, between 83 and 190 mm/year, are grouped in Northeast Thailand.

Baseflow Estimations
The median value of specific annual baseflow computed in each catchment with the local minimum method (QB estim) varies between 53 mm/year (catchment of the Nam Mun River at Rasi Salai station, Thailand) and about 1000 mm/year (catchment of the Nam Sane River at Muong Borikhan station, Laos) with a median value of 439 mm/year.Circles depicted in Figure 4 show the spatial distribution of these estimated recharge rates.Strikingly, extreme values are grouped regionally.Greatest values are observed in Central Laos, ranging from 625 to 1000 mm/year, and Southern Laos, ranging from 600 to 1000 mm/year.Lowest values, between 83 and 190 mm/year, are grouped in Northeast Thailand.

Prediction of Groundwater Recharge
The influence statistic Cooks D allowed identifying three outlier catchments in the original dataset: Nam Loei River at Ban Wang Sai station, Thailand; Stung Sangker River at Treng Station, Cambodia, and Krong Kno River at Duc Xuyen Station, Vietnam.Their removal from the regression analyses allowed increasing R 2 pred percentage value by 27 points.A combination of four explanatory variables selected among the 15 variables listed in Table 1 is sufficient to predict Q B estim with the following performances: R 2 pred = 66.36%,R 2 adj = 70.30%,and NSEC = 63.70% (Equation ( 5)).
where Q B pred is the independent variable predicting Q B estim .The dependent variables Rain (medium annual rainfall), ET 0 , Lat (latitude) and Long (longitude) are listed in Equation ( 5) according to their decreasing explanatory power (T-ratios = 7.51; −4.35; −2.96, and −2.48, respectively).The coefficient of Rain is much greater than unity, indicating that an increase of x% in annual rainfall would induce an >x% increase in baseflow (i.e., groundwater recharge).Consistently, the coefficient of ET 0 is negative, reflecting the moderating effect of evapotranspiration on groundwater recharge.Unlike ET 0 , actual land surface evapotranspiration had no explanatory power.5) is greater than unity, reflecting the non-linear relationship between rainfall depth and groundwater recharge, and highlighting the dominant role of annual rainfall in the control of groundwater recharge rates.This ratio varies between less than 15% in Northeast Thailand, up to more than 50% in Central and Southern Laos.  5) is greater than unity, reflecting the non-linear relationship between rainfall depth and groundwater recharge, and highlighting the dominant role of annual rainfall in the control of groundwater recharge rates.This ratio varies between less than 15% in Northeast Thailand, up to more than 50% in Central and Southern Laos.Groundwater recharge rates are often limited by soil properties.In addition, land cover and land use influence groundwater recharge by controlling actual evapotranspiration rates and altering soil properties, which in turn influence infiltration rates.However, no variables related to land cover and soil are included in Equation ( 5) because their explanatory power is lower than that of Rain, ET0, Lat and Long.

Model Performance
QB pred values depicted by square grid cells in Figure 4 consistently exhibit local maximums and minimums in areas where QB estim show similar extremes (e.g., Central and Southern Laos, and Northeast Thailand, respectively).Figure 7 compares QB estim and QB pred to assess the performance of the power-law model.The scatter plots align well along the first bisector with more than half of the catchments having an absolute normalized error (ANE = |QB estim−QB pred|/QB estim) lower than 30%.These errors result from the assumptions of the local minimum and multiple linear regression methods, and from possible inaccuracies in the original flow values.Even though cross-validation has been performed, extrapolation to un-gauged catchments still adds non-measurable uncertainty.Figure 8 illustrates how ANE varies according to the aridity index (i.e., ET0/Rain) and the drainage area of the studied catchments.ANE ranges from 0.26 to 0.64 across the studied catchments, with a median of 0.45, typical of humid tropical areas where regression models predicting flow are known to perform best [45].Although flow prediction is usually hampered by greater hydrological variability and higher presence of ephemeral rivers in drier areas, the power-law model predicting baseflow is not influenced by the aridity index in the LMB (Figure 8a).In contrast, ANE varies according to the drainage area of the catchments and exhibits a maximum in medium-size catchments (5000-10,000 km 2 ) (Figure 8b).Groundwater recharge rates are often limited by soil properties.In addition, land cover and land use influence groundwater recharge by controlling actual evapotranspiration rates and altering soil properties, which in turn influence infiltration rates.However, no variables related to land cover and soil are included in Equation ( 5) because their explanatory power is lower than that of Rain, ET 0 , Lat and Long.

Model Performance
Q B pred values depicted by square grid cells in Figure 4 consistently exhibit local maximums and minimums in areas where Q B estim show similar extremes (e.g., Central and Southern Laos, and Northeast Thailand, respectively).Figure 7 compares Q B estim and Q B pred to assess the performance of the power-law model.The scatter plots align well along the first bisector with more than half of the catchments having an absolute normalized error (ANE = |Q B estim −Q B pred |/Q B estim ) lower than 30%.These errors result from the assumptions of the local minimum and multiple linear regression methods, and from possible inaccuracies in the original flow values.Even though cross-validation has been performed, extrapolation to un-gauged catchments still adds non-measurable uncertainty.Figure 8 illustrates how ANE varies according to the aridity index (i.e., ET 0 /Rain) and the drainage area of the studied catchments.ANE ranges from 0.26 to 0.64 across the studied catchments, with a median of 0.45, typical of humid tropical areas where regression models predicting flow are known to perform best [45].Although flow prediction is usually hampered by greater hydrological variability and higher presence of ephemeral rivers in drier areas, the power-law model predicting baseflow is not influenced by the aridity index in the LMB (Figure 8a).In contrast, ANE varies according to the drainage area of the catchments and exhibits a maximum in medium-size catchments (5000-10,000 km 2 ) (Figure 8b).While ANE allows the predictive performance of the models to be assessed for individual catchments and to determine how it relates to the catchments characteristics, the three criteria (R 2 pred = 66.36%;R 2 adj = 70.30%;NSEC = 63.70%)assess how well the power-law model described in Equation ( 5) performs, allowing comparisons with regional regression models developed in other parts of the world.Our NSEC value, equivalent to the classical R 2 reported in [45], falls within the range of values typically observed for regression models predicting low flows in other humid regions of the world [45].Finally, it should be noted that we re-applied the regression analysis using the original values of the variables (not log-transformed) and verified that the power-law structure outperformed the linear one, likely because of the typical non-linear relationship between rainfall and flow.

Factors Determining Groundwater Recharge
Like many other studies under temperate or tropical climates [46][47][48], our results indicate that spatial variations in baseflow, used as a proxy of groundwater recharge, are predominantly controlled by rainfall and ET0.While ET0 is the second most powerful predictive variable of Equation ( 5), actual land surface evapotranspiration had no predictive power.Likely explanations include: (i). the mismatch between periods used to compute QB estim and actual evapotranspiration derived from MODIS 16, and available since 2000 only, considering that recent land-cover changes have occurred across the region; (ii).the incomplete validation of MODIS 16 product for tropical Southeast Asia [49].While ANE allows the predictive performance of the models to be assessed for individual catchments and to determine how it relates to the catchments characteristics, the three criteria (R 2 pred = 66.36%;R 2 adj = 70.30%;NSEC = 63.70%)assess how well the power-law model described in Equation ( 5) performs, allowing comparisons with regional regression models developed in other parts of the world.Our NSEC value, equivalent to the classical R 2 reported in [45], falls within the range of values typically observed for regression models predicting low flows in other humid regions of the world [45].Finally, it should be noted that we re-applied the regression analysis using the original values of the variables (not log-transformed) and verified that the power-law structure outperformed the linear one, likely because of the typical non-linear relationship between rainfall and flow.

Factors Determining Groundwater Recharge
Like many other studies under temperate or tropical climates [46][47][48], our results indicate that spatial variations in baseflow, used as a proxy of groundwater recharge, are predominantly controlled by rainfall and ET0.While ET0 is the second most powerful predictive variable of Equation ( 5), actual land surface evapotranspiration had no predictive power.Likely explanations include: (i). the mismatch between periods used to compute QB estim and actual evapotranspiration derived from MODIS 16, and available since 2000 only, considering that recent land-cover changes have occurred across the region; (ii).the incomplete validation of MODIS 16 product for tropical Southeast Asia [49].While ANE allows the predictive performance of the models to be assessed for individual catchments and to determine how it relates to the catchments characteristics, the three criteria (R 2 pred = 66.36%;R 2 adj = 70.30%;NSEC = 63.70%)assess how well the power-law model described in Equation ( 5) performs, allowing comparisons with regional regression models developed in other parts of the world.Our NSEC value, equivalent to the classical R 2 reported in [45], falls within the range of values typically observed for regression models predicting low flows in other humid regions of the world [45].Finally, it should be noted that we re-applied the regression analysis using the original values of the variables (not log-transformed) and verified that the power-law structure outperformed the linear one, likely because of the typical non-linear relationship between rainfall and flow.

Factors Determining Groundwater Recharge
Like many other studies under temperate or tropical climates [46][47][48], our results indicate that spatial variations in baseflow, used as a proxy of groundwater recharge, are predominantly controlled by rainfall and ET 0 .While ET 0 is the second most powerful predictive variable of Equation ( 5), actual land surface evapotranspiration had no predictive power.Likely explanations include: (i) the mismatch between periods used to compute Q B estim and actual evapotranspiration derived from MODIS 16, and available since 2000 only, considering that recent land-cover changes have occurred across the region; (ii) the incomplete validation of MODIS 16 product for tropical Southeast Asia [49].Latitude and longitude are the third and fourth most powerful predictive variables of Table 1.They most likely act as surrogates for environmental processes controlling baseflow, exhibiting latitudinal and longitudinal gradients, and not listed in Table 1 since independency between the explanatory variables is a prerequisite for inclusion in the regression equation.Regional characteristics of the aquifers (Figure 2) exhibit such gradients: the southwest part of the LMB corresponds to sedimentary sandy deposits of the Khorat Plateau in Northeast Thailand and the alluvial plain of the Mekong River in Cambodia and Southern Vietnam.These sedimentary deposits are usually characterized by greater permeability than the more compact metamorphic rocks prevalent along the Annamite Range from the North to the Southeast of the LMB [3,18,19].Although this contrast is reversed in a few and confined locations (e.g., the basalt flows in Southern Laos with greater permeability; the early Paleozoic rocks in Northeast Thailand with reduced permeability), these local contrasts have likely minor influence on groundwater recharge at the scale of the studied catchments.Based on these observations, and accounting for the limited quantitative hydrogeological data explained in Section 2, we hypothesize that Lat and Long in Equation ( 5) are surrogate variables for the permeability and transmissivity of the aquifers that likely increase from North to South and from East to West across the LMB, in accordance with the signs of the coefficients of Lat and Long in Equation (5).It should be noted that a power-law Equation with only Rain and ET 0 as explanatory variables yields an R 2 pred value of 63.69%, about 3 points lower than the R 2 pred value of Equation ( 5).This comparison indicates that: (i) climate explains more than half of Q B estim variability across the LMB; and (ii) assuming that Lat and Long are surrogates for aquifer properties, the regional geology explains at least 3% of Q B pred variability.
The non-inclusion of land-cover and soil characteristics in Equation ( 5) is due to their explanatory power, lower than that of the four selected variables Rain, ET 0 , Long and Lat.The lower hydrological influence associated with land covers is consistent with the usually moderate hydrological effect of land-cover changes in catchments with mixed land covers and an area greater than 1000 km 2 .Over such large areas, the combinations of various land covers, with counteracting changes, generally render their individual hydrological effects difficult to detect [44].The exclusion of soil characteristics from Equation ( 5) may be related to the poor accuracy of the soil maps used in this assessment, and/or to the surrogating effect of the latitude and longitude.

Comparison with Previous Studies
The range of recharge rates that we estimated and predicted in the LMB (Figure 4) is similar to that provided in the global assessment calculated by the WaterGAP Global Hydrology Model WGHM [5], and varying between 100 and 1000 mm/year across mainland Southeast Asia.Our estimations of recharge rates can also be compared with more localized studies [8][9][10][50][51][52][53][54].Baseline recharge rates estimated with the HELP3 hydrologic model and projected under the ECHAM GCM A2 and B2 scenarios were used to predict groundwater recharge rates in the Huai Khamrian catchment in Northeast Thailand, yielding values of 253, 351, and 329 mm/year, respectively [8].While these estimations exceed our conservative estimates of less than 170 mm/year in Northeast Thailand (Figure 4), they confirm that recharge rates in this region are positively correlated to rainfall, in agreement with Equation (5).The ratios (recharge/rainfall) obtained with the HELP3 model in the Huai Khamrian catchment (around 16%) are well aligned with our estimations mapped in Figure 6, suggesting that rainfall values used by [8] are greater than those in Aphrodite product, used in Equation (5).In Northern Thailand, groundwater recharge rates corresponding to different land uses were modeled in a 6500 km 2 catchment with the WetSpass model [50].Resulting annual recharge rates averaged 360 mm/year, mostly influenced by rainfall and evapotranspiration.In Southeast Vietnam, measurements from 10 monitored wells were used to infer groundwater recharge using finite difference methods [10].Estimated recharge rates varied between 307 and 325 mm/year, slightly below our estimations in this part of Vietnam (500 mm/year).Similar consistency was observed in Northern Vietnam where our estimates (315 mm/year) were moderately exceeded by a mean rate of 477 mm/year derived from the rainfall infiltration breakthrough model calibrated with measurements of rainfall and groundwater levels [9].In the Day river sub-basin of the Red River in Northern Vietnam, Van der Wolf [51] calibrated the SWAT model using hydro-meteorological observations and detailed maps of land uses, topography and soils.The SCS curve number method was used to model surface runoff and to infer infiltration rates.The resulting recharge rates averaged 248 mm/year, ranging between 37 and 601 mm/year, in agreement with our results (213-318 mm/year) for this area.Broader-scale groundwater recharge assessments was performed over 15,000 km 2 in central Cambodia [52], yielding 448 mm/year from the SCS runoff curve method, aligned with Q B estim (464 mm/year) and slightly greater than Q B pred (278-357 mm/year) in this area.In Northwest Cambodia, groundwater recharge rates were estimated over 3375 km 2 of sandstones from Upper Triassic to Lower Cretaceous, using the water-table fluctuation method and the stable isotopes analysis from 12 piezometers [53].Recharge rates of 10 to 70 mm/year, lower than our estimations in this area (200 mm/year), were explained by clayey soils overlying sandstone, whose presence can highly vary at the scale of few tens of meters.Similar low estimates (20 mm/year) were derived from the groundwater flow model MODFLOW in Southeast Cambodia, much lower than our estimations in this location (230 mm/year), and attributed to the presence of a surficial clay aquitard [54].

Limitations of the Study
The comparison of our results with previous studies undertaken in the LMB (Section 5.2) has highlighted an overall general agreement: absolute relative differences between these previous local assessments and our recharge rates in these locations are yielding an average of 32% excluding the two case studies in North and Southeast Cambodia where this percentage exceeds 100 [53,54].These higher discrepancies likely result from two main factors: (i) the presence of local clay deposits in the surveyed sites; (ii) the sparse network of river gauging stations in Cambodia not monitoring these specific sites, thus not allowing our model to capture this sub-regional recharge constraint.Nonetheless, the differences pointed out in this comparison could also originate from inaccuracies in the reported studies that cannot be readily addressed.
Two counteracting processes can influence the relationship between the ANE values of the predicted recharge rates and catchment size (Figure 8b).ANE of low flow models is usually lower in larger catchments due to greater space-time aggregation of runoff processes [38].However, in larger catchments, specific baseflow (mm/year) tends to increase because of increased seepage [55].Other scaling issues include groundwater lateral fluxes at the margin of the studied catchments, which cannot be accounted in a vertical analysis.These trends possibly explain the nonlinear relationship between ANE and the catchment drainage areas (Figure 8b).
Although the underlying physical processes associated to the geographic coordinates could not be ascertained, their inclusions as explanatory variables significantly improves the predictive power of Equation ( 5) and the accuracy of the groundwater recharge map in Figure 4.This map gives a reasonably reliable picture of the regional variations in groundwater recharge and their main controlling factors, namely annual rainfall and standard evapotranspiration.

Groundwater Potential for Irrigation
The recharge rates mapped in Figure 4 can help assess potential groundwater resources available for irrigation and other purposes.Groundwater is a largely untapped resource for agricultural development in the LMB and particularly in Laos where groundwater irrigation in 2010 represented just 0.1% of the total irrigated area in the country [2].However, in areas remote from reliable surface water supplies, farmers are increasingly resorting to this resource to irrigate high-value crops [4,56].Although there is considerable scope for extending groundwater-fed irrigations schemes, setting an upper limit can help preserve groundwater resources and maintain their eco-systemic functionalities.Indeed, over-abstraction depletes water tables leading to endogenic contamination (e.g., salinization on the Khorat Plateau) and groundwater shortages.Due to the continuum between surface and groundwater resources, particular attention must be paid to abnormally low baseflow threatening downstream ecosystems and other downstream water uses, influenced by aquifer levels in upstream areas.While even moderate groundwater extraction will systematically impact water table levels, it is useful to assess a tolerance threshold under which the socio-economic benefits derived from groundwater use outweigh the costs.
We assessed a first-order approximation of an upper limit of irrigable area using a simple groundwater balance [57] applied to two regions of the LMB exhibiting high potential for irrigation development: the Vientiane Plain in Laos and Northeast Thailand.These two regions present similarly flat topography, limited surface water supply during the dry season, and good road connections to local markets.Median groundwater recharge rates yield 150 mm/year in Northeast Thailand and 350 mm/year in the Vientiane Plain (Figure 4).Referring to the typical crops grown in the region (e.g., rice, sugar cane, vegetables), a conservative (i.e., higher estimate) crop water demand for one cycle in the dry season would amounts about 1000 mm/year according to FAO [58].Assuming that half of the annual ground recharge rate should be reserved to service other needs (environment, domestic consumptions, industries, livestock), the percentage of irrigable area I % can be estimated as follows: I % = 0.5 × Q B /1000 [57].In Northeast Thailand, I % = 7.5%.This percentage area is equivalent to the fraction of agricultural land actually irrigated with surface water in this region [59], suggesting that irrigation could potentially be doubled by improving groundwater access.However, any change in the local groundwater balance can have detrimental effects locally (e.g., increased groundwater salinity), and downstream (e.g., altered water uses and ecosystems).In the Vientiane Plain, I % = 17.5% while the current percentage area of irrigated land in the Vientiane Prefecture, which largely covers the southern boundary of the Vientiane Plain, is around 10% (unpublished sources from Department of Irrigation, Ministry of Agriculture and Forestry, Laos).This figure, although conservative, demonstrates the considerable groundwater potential for developing irrigation, albeit with numerous technical and non-technical issues which severely constrain development [60].

Conclusions
Based on estimations of baseflow derived from the local-minimum filtering method, a regional regression model was developed to map the spatial distribution of groundwater recharge rates in the LMB.Results indicate that spatial variations in groundwater recharge are predominantly controlled by the climate (rainfall and evapotranspiration) at the scale of the LMB.While this study confirms that large areas exist for agricultural groundwater development, the proposed map provides a simple way to assess the likely limit of sustainable groundwater-fed irrigation, useful for broad-scale water resources planning.Compared to regional regression models developed in other parts of the world, the power-law model that we developed to predict groundwater recharge rates performs reasonably well.Recharge estimates compare favorably with estimates from the few other local studies that used different recharge assessment approaches, though with greater inaccuracies in some instances, due to local-scale heterogeneities in soils and geology.While extensive field measurements of biophysical variables across the LMB will help to improve the performance and physical basis of the power-law model presented in this paper, further local-scale analyses are required for more detailed assessments and policy development.

Figure 1 .
Figure 1.Topography of the Lower Mekong Basin.Elevations are derived from HydroSHEDS.

Figure 1 .
Figure 1.Topography of the Lower Mekong Basin.Elevations are derived from HydroSHEDS.
discharging into tributaries of the Mekong River.In some locations, the Permianaged karstic limestone are considered to be more productive aquifers; (3) In the lowlands, and particularly in Northeast Thailand, deep confined and shallow unconfined aquifers from the Mesozoic are comprised of sandstones.Outcrops of early Paleozoic rock bounded by Upper Jurassic sandstones have lower water yields; (4) In the Mekong delta, Cenozoic alluvial and deltaic sediments of up to 800 m thick form both unconfined and confined aquifers.The alluvial deposit in the upstream narrower valley rarely exceeds 100 m in thickness.Recharge in Pleistocene through Miocene deposits occurs in terraces along the margins of the Holocene alluvial plain.These terraces are generally coarser-textured and more productive aquifers than fine-textured Holocene surficial sediments including large amount of clay.While these observations provide an overall view of the hydrogeology across the LMB, the level of understanding remains largely qualitative.The transmissivities and water storage capacities of aquifers in the LMB have yet to be measured in a systematic way[3].Only localized pumping tests allowed estimating aquifer hydraulic properties over a few areas[20,21].

Figure 3 .
Figure 3. Gauged sub-catchments of the Lower Mekong Basin.

Figure 3 .
Figure 3. Gauged sub-catchments of the Lower Mekong Basin.

Figure 4 .
Figure 4. Median annual groundwater recharge in the Lower Mekong Basin.Circles are located in the centroid of the gauged catchments and their corresponding values are estimated with the local minimum method (QB estim).Un-gauged values in square grid cells are predicted by multiple log-linear regressions (QB pred).The graduated color-scale indicates the values for both QB estim and QB pred.

4. 2 .
Multiple Regressions Analysis 4.2.1.Prediction of Groundwater Recharge The influence statistic Cooks D allowed identifying three outlier catchments in the original dataset: Nam Loei River at Ban Wang Sai station, Thailand; Stung Sangker River at Treng Station, Cambodia, and Krong Kno River at Duc Xuyen Station, Vietnam.Their removal from the regression

Figure 4 .
Figure 4. Median annual groundwater recharge in the Lower Mekong Basin.Circles are located in the centroid of the gauged catchments and their corresponding values are estimated with the local minimum method (Q B estim ).Un-gauged values in square grid cells are predicted by multiple log-linear regressions (Q B pred ).The graduated color-scale indicates the values for both Q B estim and Q B pred .
Lat and Long are negatively correlated to baseflow.The exclusion of actual ET and inclusion of Lat and Long in Equation (5) are discussed in Section 5. Predicted values of groundwater recharge Q B pred derived from Equation (5) are represented with colored grid cells in Figure 4. To prepare this map at 0.25 • × 0.25 • spatial resolution, each value of ET 0 available at 0.5 • × 0.5 • spatial resolution was first replicated in the four corresponding 0.25 • × 0.25 • grid cells.Figure 5 displays regional variations in Rain and ET 0 .Visual comparison with Figure 4 confirms that Rain is the main driver of Q B pred , with maximum and minimum Q B pred values observed in the rainiest and driest locations, respectively.Anti-correlation between Q B pred and ET 0 is less obvious, though highly significant (F test p-value = 0.00), as confirmed by the correlation coefficients R Ln computed on the logarithms of the variables: R Ln (Q B pred , Rain) = 79.13%;R Ln (Q B pred , ET 0 ) = −32.10%.

Figure 6
maps spatial variations in the ratio between Q B pred and Rain.These variations follow the regional variations of Rain Figure5a, R Ln (Rain, Q B pred /Rain) = 52.06%,because the coefficient of the dependent variable Rain in Equation (

Hydrology 2017, 4 , 60 10 of 18 decreasing
explanatory power (T-ratios = 7.51; −4.35; −2.96, and −2.48, respectively).The coefficient of Rain is much greater than unity, indicating that an increase of x% in annual rainfall would induce an >x% increase in baseflow (i.e., groundwater recharge).Consistently, the coefficient of ET0 is negative, reflecting the moderating effect of evapotranspiration on groundwater recharge.Unlike ET0, actual land surface evapotranspiration had no explanatory power.Lat and Long are negatively correlated to baseflow.The exclusion of actual ET and inclusion of Lat and Long in Equation (5) are discussed in Section 5. Predicted values of groundwater recharge QB pred derived from Equation (5) are represented with colored grid cells in Figure4.To prepare this map at 0.25° × 0.25° spatial resolution, each value of ET0 available at 0.5° × 0.5° spatial resolution was first replicated in the four corresponding 0.25° × 0.25° grid cells.Figure5displays regional variations in Rain and ET0.Visual comparison with Figure4confirms that Rain is the main driver of QB pred, with maximum and minimum QB pred values observed in the rainiest and driest locations, respectively.Anti-correlation between QB pred and ET0 is less obvious, though highly significant (F test p-value = 0.00), as confirmed by the correlation coefficients RLn computed on the logarithms of the variables: RLn (QB pred, Rain) = 79.13%;RLn (QB pred, ET0) = −32.10%.Figure6maps spatial variations in the ratio between QB pred and Rain.These variations follow the regional variations of Rain Figure (5)a, RLn (Rain, QB pred/Rain) = 52.06%,because the coefficient of the dependent variable Rain in Equation (

Figure 5 .
Figure 5. Spatial variability of mean climate variables across the Lower Mekong Basin.(a) Median annual rainfall derived from Aphrodite [36].(b) Standard evapotranspiration derived from CRU [38].

Figure 5 .
Figure 5. Spatial variability of mean climate variables across the Lower Mekong Basin.(a) Median annual rainfall derived from Aphrodite [36].(b) Standard evapotranspiration derived from CRU [38].

Figure 6 .
Figure 6.Ratio between estimated median annual recharge (QB estim) and median annual rainfall (Rain) in the Lower Mekong Basin.

Figure 6 .
Figure 6.Ratio between estimated median annual recharge (Q B estim ) and median annual rainfall (Rain) in the Lower Mekong Basin.

Figure 7 .
Figure 7.Comparison of observed (QB estim j) and predicted (QB pred j) specific median annual baseflow in each gauged catchment j.

Figure 7 . 18 Figure 7 .
Figure 7.Comparison of observed (Q B estim j ) and predicted (Q B pred j ) specific median annual baseflow in each gauged catchment j.

Table 1 .
Candidate explanatory variables considered in the multiple regression analyses: Variation ranges across the 65 catchments.
* See Table2for the description of the 4-unit scale.

Table 2 .
The 4-unit scale of the two soil characteristics.