Estimating Baseﬂow and Baseﬂow Index in Ungauged Basins Using Spatial Interpolation Techniques: A Case Study of the Southern River Basin of Thailand

: This research aims to estimate baseﬂow (BF) and baseﬂow index (BFI) in ungauged basins in the southern part of Thailand. Three spatial interpolation methods (namely, inverse distance weighting (IDW), kriging, and spline) were utilized and compared in regard to their performance. Two baseﬂow separation methods, i.e., the local minimum method (LM) and the Eckhardt ﬁlter method (EF), were investigated. Runoff data were collected from 65 runoff stations. These runoff stations were randomly selected and divided into two parts: 75% and 25% for the calibration and validation stages, respectively, with a total of 36 study cases. Four statistical indices including mean absolute error (MAE), root mean squared error (RMSE), correlation coefﬁcient (r), and combined accuracy (CA), were applied for the performance evaluation. The ﬁndings revealed that monthly and annual BF and BFI calculated by EF were mostly lower than those calculated by LM. Furthermore, IDW gave the best performance among the three spatial interpolation techniques by providing the highest r-value and the lowest MAE, RMSE, and CA values for both the calibration and validation stages, followed by kriging and spline, respectively. We also provided monthly and annual BF and BFI maps to beneﬁt water resource management.


Introduction
Baseflow (BF) is a complete streamflow portion that slowly flows into a stream from the saturated soil or groundwater storage [1,2] and predominantly contributes to streamflow during the dry season. It is crucial to understand the hydrological characteristics, especially the spatiotemporal variation of BF availability in the watershed, to plan and monitor water resources and ecological systems. BF also helps us understand the hydrology of the watershed in terms of surface and subsurface water interactions, urbanization effects on runoff generation, and healthy aquatic habitats within a stream. Owing to the lack of directly measured BF data in general, the separation method has been recognized for determining BF. Recently, many researchers have applied and compared the performance of several different separation methods to obtain BF [3][4][5][6]. Eckhardt [3] recommended a two-parameter filter to be more reasonable than a one-parameter filter and indicated that the maximum baseflow index (BFI max ) values depended on the watershed's hydrological and hydrogeological characteristics. Seven different methods of BF separation (namely, HYSEP1, HYSEP2, HYSEP3, PART, BFLOW, UKIH, and Eckhardt) were investigated by water resource issues. Shyamala et al. [25] studied spatial interpolation techniques for groundwater pollution. They noticed that kriging was the best approach compared to IDW and spline. Yan et al. [26] applied the ordinary kriging interpolation method to characterize soil thickness on various land use types in the Yimeng Mountain area of China. Some research have shown that regression kriging can significantly boost spatial prediction accuracy when providing a weakly correlated auxiliary variable [27]. Ly et al. [28] reviewed various approaches for the spatial interpolation of rainfall data for operational hydrology and watershed modeling. They indicated that for point-by-point assessment, interpolation using the IDW and kriging methods is more efficient than simple Thiessen and spline techniques, particularly for a monthly time period. Li and Heap [29] presented guidance and suggestions for the application of three types of spatial interpolation techniques to environmental tasks including, the non-geostatistical interpolation approach, geo-statistical interpolation approach and combined approach. Apart from the quantitative approach, Wu and Hung [30] suggested a visualization approach for the evaluation of spatial interpolation techniques.
From reviewing the literature for the applicability of spatial interpolation techniques in estimating BF and BFI in ungauged watersheds, we found only one study that used the IDW and kriging approaches to interpolate BFI with a grid resolution of 1000 m for the conterminous United States [31]. The actual and interpolated BFI comparison was undertaken, and the analysis results showed that the error could vary from 12% to 22%, depending on the US region. Our research work was to present and expand the study case of applying and finding the most suitable spatial interpolation techniques to estimate BF and BFI in ungauged basins. To achieve our goal, we had three main objectives as follows.

1.
To compare the two separation methods, namely, the local minimum method (LM) and the Eckhardt filter method (EF), for estimating BF and BFI.

2.
To evaluate the efficacy of three spatial interpolation approaches (IDW, kriging, and spline) in estimating BF and BFI in ungauged basins.

3.
To create maps showing the BF and BFI's spatial and temporal variation as useful information for supporting water management-related agencies.
Our investigation results gave more accurate information on the spatial and temporal variation of BF and BFI, which is valuable information for water resource planning and management in our region.

Study Area and Data Used
Our study focused on the southern river basin of Thailand, which covers five main river basins: the Peninsular-East Coast, the Peninsular-West Coast, Mae Nam Tapi, Thale Sap Songkhla, and Mae Nam Pattani, as depicted in Figure 1. These river basins have areas ranging from 13 to 6713 km 2 in size. This area is geographically a peninsula between the Andaman Sea's western side and the South China Sea's eastern side. It is also extended in the northern and central parts by the long ridge on the west of the mountains. This section is divided into two regions by the Phuket ridge along the west coast and the Nakhon Si Thammarat ridge in the central lower portion of the southern backbone: the east and west coasts. The climate variability on both sides of the river basins was primarily dominated by the northeast monsoon winds and southwest monsoon. The southwestern monsoon wind usually begins in mid-May and ends in mid-October, while the northeastern monsoon typically starts in mid-October and ends in mid-February. The average annual rainfall in this area is approximately 2291 mm, and the average monthly rainfall is between 9 mm and 670 mm. The runoff data used in this analysis were obtained from the daily reports of 65 runoff stations with periods ranging from 4 to 12 years. The annual runoff average is approximately 853 MCM, and the monthly runoff average is between 0.29 MCM and 716 MCM. Table 1 shows the meteorological and hydrological characteristics of the study area. mm and 670 mm. The runoff data used in this analysis were obtained from the daily reports of 65 runoff stations with periods ranging from 4 to 12 years. The annual runoff average is approximately 853 MCM, and the monthly runoff average is between 0.29 MCM and 716 MCM. Table 1 shows the meteorological and hydrological characteristics of the study area.

Baseflow Separation
In our study, the Web-based Hydrograph Analysis Tool (WHAT) was applied as a method for performing BF separation [32]. It includes one graphical method (local minimum approach) and two digital filter approaches (a digital filter with a one-parameter approach [33] and a digital filter with a two-parameter approach, widely recognized as the Eckhardt filter method). This research applied and compared the results of using the local minimum method and the Eckhardt filter method. Previously, there was research

Baseflow Separation
In our study, the Web-based Hydrograph Analysis Tool (WHAT) was applied as a method for performing BF separation [32]. It includes one graphical method (local minimum approach) and two digital filter approaches (a digital filter with a one-parameter approach [33] and a digital filter with a two-parameter approach, widely recognized as the Eckhardt filter method). This research applied and compared the results of using the local minimum method and the Eckhardt filter method. Previously, there was research trying to observe the comparison between the Eckhardt filter method and seven other digital filter methods [4]. The description of the local minimum method and the Eckhardt filter method is as follows: The LM is one of the graphics approaches used to separate the BF from the total streamflow. Graphics approaches include fixed interval, sliding interval, and local minimum, as introduced by Sloto and Crouse [34]. For each time interval, the LM uses the idea of searching for the two minimum flows and then connecting them with a straight line to separate the BF section from the total runoff. In this regard, the time interval is calculated to determine the lowest discharge within one half the interval minus one day [(2N − 1)/2d] [15] before and after the considered day. The N value is equivalent to A 0.2 , where A is the basin area in square miles [35].

Eckhardt Filter Method (EF)
The focus of the EF is on recursive digital filters (RDFs) that apply signal analysis and processing theory to separate the low-frequency signal (BF) from the high-frequency signal (quick flow). There are two parameters, consisting of the filter parameter and BFI max , for the Eckhardt approach. The filter parameter means the rate at which the streamflow decreases and can be determined directly by a recession analysis after a recharging event. BFI max is the maximum BFI that a recursive digital filter algorithm can use to create a model. Mathematically, the EF can be expressed as: where Q b,t and Q b,t−1 are the baseflow at step t and t − 1 time, Q s,t is the complete streamflow at step t time; and an is the filter parameter. For the first step in the calculation, the value of Q b,t−1 was assumed to be 50% of the total flow as suggested by Zhang, Ahiablame, Engel and Liu [21]. However, empirical judgement would be tried to represent physical reality. Eckhardt [3] recommended representative BFI max values for three hydrological and hydrogeological conditions: BFI max = 0.80, for perennial streams with porous aquifers; BFI max = 0.50, for ephemeral streams with porous aquifers; and BFI max = 0.25, for perennial streams with hard rock aquifers. The default BFI max of 0.80 and filter parameters of 0.98 were used in our study, as recommended by WHAT, referring to the hydrological and geological characteristics of the southern river basin of Thailand, a perennial stream of porous aquifers.

Inverse Distance Weighting (IDW)
The inverse distance weighting (IDW) method is based on the principle that locations are closely related to each other, spatially, calculating the value at the desired location where the closest station position will have greater importance. It is an approximation of an unknown point from the linear sum of the known values and the weighted point to be limited by distance. With the distance from the unknown point to the next known point, this weighting is altered and is formulated as follows: (2) where Z represents the position of a known point (i) and the predicted value at the unknown point (j), n denotes the exponent that determines the rate of weight reduction as the distance increases (e.g., 1, 2, 3), and d ij indicates the distance from a known point (i) to an unknown point (j).

Kriging
Kriging proceeds the estimation by giving the weights of the averaged input values. This is similar to the moving-average method for calculating the weight. A semi-variogram Water 2021, 13, 3113 6 of 15 model is used to determine the data's spatial relationship by analyzing the relation between the squared differences of pairs of point values and distance. Therefore, it is necessary to test whether the data are suitable for the model. Kriging is not distance-weighted between the known and unknown positions. Rather, it is the grouping of known positions according to the spatial relationship characteristics that are related to each other and finding the variation to be used as a weight. There are several equations for modulating variograms. Each equation has an initial value of variation (nugget) and the level of the variogram ends. Alternatively, the default constant (sill) and the distance from each point to the sill (range) vary.
The approximation from the variogram was used to estimate the distance-weighted mean when estimating the spatial data. The resulting value is the sum of the known points' weighted values, depending on the distance between the estimated point and the known point. The weights selected for estimation were not biased. In this study, we herein applied a universal kriging approach. Universal kriging has the form of deterministic interpolation. It hypothesizes that the spatial variation in z values is combined and has spatial relationships with known points. In addition, universal kriging is a method of modulating the area curvature to integrate a plane surface with a quadratic surface, which takes the form of a polynomial equation. It can be described as: where M is the weight that has a relationship between the point to be approximated and the known point (x i ), y i is the distance between the points and b 1 and b 2 are the number of pairs of points each by distance h.

Spline
Using a mathematical function to minimize the overall surface curvature, a spline is a deterministic technique used to define two-dimensional curves on three-dimensional surfaces, resulting in a smooth surface that precisely passes through the input points or known points. It can create fairly distinct and visually appealing features using just a few sample points. However, the disadvantages are the resulting surfaces with different minimum and maximum values from the collection of input data. It is prone to outliers, and there is no error sign, which is similar to IDW. For surface interpolation, the algorithm used for the spline utility uses the following formula: where j = 1, 2,..., N; N is the number of points; λ j are coefficients found by the solution of a linear equation system, and r j is the distance from point (x, y) to the jth point. T(x, y) and R(r) are defined differently depending on the selected option. The spline method is a mathematical equation suitable for surfaces with gradual changes, such as the height and depth of the water surface, etc. We chose the regularized spline in this study. A regularized spline is a technique for smoothing the results, where the value of information has a more gradual increase or decrease. The optimal weight value should be between 0-0.5. The regularized spline function is taken as the basis function: T(x, y) = a 1 +a 2 x + a 3 y where a i is the coefficient found from the solution of a system of linear equations, and where r is the distance between the point and the sample, τ 2 is the weight parameter, K 0 is the modified Bessel function and "C" is a constant equal to 0.577215.

Results and Discussion
The present study deals with estimating BF and BFI in ungauged basins in the southern river basin of Thailand. Three major results, as indicated in our research objective, are presented and discussed as follows: (1) comparative results of estimating BF and BFI using two separation methods (i.e., LM and EF methods); (2) evaluation of the efficacy of three spatial interpolation approaches (i.e., IDW, kriging, and spline) in estimating BF and BFI in ungauged basins; and (3) maps showing the spatial and temporal variations of BF and BFI. The details are presented as follows:

Local Minimum Method vs. Eckhardt Filter Method for Estimating BF and BFI
The spatial variation of the annual runoff and BF performed by the LM and EF methods for the 65 runoff stations is shown in Figure 2. The BF computed by the LM method had a bit higher value than that calculated by the EF method. Additionally, Figure 3 shows the BFI for 65 runoff stations calculated using the LM and EF methods. We generally found the same trend as the BF; that is, the BFI calculated by the LM method had a higher value than that calculated by the EF method. It is because LM's determination process of the lowest discharge to separate baseflow from total runoff gives the local minimums and depends on the calculated time interval as a function of the basin area. Additionally, the baseflow values for each day between local minimums are estimated by linear interpolations. Therefore, it may result in overestimation in BF and BFI. On the other hand, the EF considers the shorter interval (day by day in this experiment), obtaining the smoother and more refined baseflow values.
where r is the distance between the point and the sample, τ 2 is the weight parameter, K0 is the modified Bessel function and "C" is a constant equal to 0.577215.

Results and Discussion
The present study deals with estimating BF and BFI in ungauged basins in the southern river basin of Thailand. Three major results, as indicated in our research objective, are presented and discussed as follows: (1) comparative results of estimating BF and BFI using two separation methods (i.e., LM and EF methods); (2) evaluation of the efficacy of three spatial interpolation approaches (i.e., IDW, kriging, and spline) in estimating BF and BFI in ungauged basins; and (3) maps showing the spatial and temporal variations of BF and BFI. The details are presented as follows:

Local Minimum Method vs. Eckhardt Filter Method for Estimating BF and BFI
The spatial variation of the annual runoff and BF performed by the LM and EF methods for the 65 runoff stations is shown in Figure 2. The BF computed by the LM method had a bit higher value than that calculated by the EF method. Additionally, Figure 3 shows the BFI for 65 runoff stations calculated using the LM and EF methods. We generally found the same trend as the BF; that is, the BFI calculated by the LM method had a higher value than that calculated by the EF method. It is because LM's determination process of the lowest discharge to separate baseflow from total runoff gives the local minimums and depends on the calculated time interval as a function of the basin area. Additionally, the baseflow values for each day between local minimums are estimated by linear interpolations. Therefore, it may result in overestimation in BF and BFI. On the other hand, the EF considers the shorter interval (day by day in this experiment), obtaining the smoother and more refined baseflow values.   5 show an example of the temporal variation of the total runoff, BF, and BFI at the X.248 runoff station. We observed that almost every month of the year, BF, and BFI were given by LM had higher values than those provided by EF, except in July, August, and November. During these three months, X.248 runoff station is influenced by the southwest monsoon wind (July and August) and the northeast monsoon wind (November), resulting in a high runoff. As explained previously, the EF method gives smoother and more refined baseflow values than those provided by the LM method. As a result, the LM method may produce underestimation in BF and BFI since the baseflow values for each day between local minimums are estimated by linear interpolations.  5 show an example of the temporal variation of the total runoff, BF, and BFI at the X.248 runoff station. We observed that almost every month of the year, BF, and BFI were given by LM had higher values than those provided by EF, except in July, August, and November. During these three months, X.248 runoff station is influenced by the southwest monsoon wind (July and August) and the northeast monsoon wind (November), resulting in a high runoff. As explained previously, the EF method gives smoother and more refined baseflow values than those provided by the LM method. As a result, the LM method may produce underestimation in BF and BFI since the baseflow values for each day between local minimums are estimated by linear interpolations.   Using two statistical indices (namely, mean difference (MD) and correlation coefficient (r)), a distinction was made between these two annual and monthly BF and BFI methods, as shown in Table 2. This indicates that the BF and BFI, as computed by the EF method, had lower values than those values as computed by the LM method for both annual and monthly periods. The results are similar compared to those found by Schilling and Jones [7]. It can be noted from the acquisition of MD's negative values that the discrepancy between the values given by EF and those given by LM, except for the BF result in April, gave a positive value (0.06 MCM) that was an inverse meaning, as discussed previously. In October, the maximum difference in value of the BF and BFI was −4.31 MCM and −0.09, respectively. In addition, there was a very high correlation between the two approaches, i.e., r ≥ 0.95 for BF and r ≥ 0.70 for BFI.  Using two statistical indices (namely, mean difference (MD) and correlation coefficient (r)), a distinction was made between these two annual and monthly BF and BFI methods, as shown in Table 2. This indicates that the BF and BFI, as computed by the EF method, had lower values than those values as computed by the LM method for both annual and monthly periods. The results are similar compared to those found by Schilling

Performance Comparison of Spatial Interpolation Techniques
To investigate the performance of three spatial interpolation techniques (namely, IDW, kriging, and spline), we randomly selected our data sets three times with 75% and 25% of the total 65 runoff stations for the calibration and validation stages, respectively. Thus, we had three different repeated samples for each considered period, i.e., April representing a dry season, November representing a monsoon season, and annual, for BF and BFI. Using two separation methods, i.e., LM and EF, we finally had a total of 36 study cases (3 cases by randomly selecting three times from the whole data sets, 3 cases from three considered periods, 2 cases by having BF and BFI, and 2 cases for two baseflow separation methods) in our experiment. The information in Tables 2 and 3 shows a summary of the comparative statistical metrics, that is, mean absolute error (MAE), root mean squared error (RMSE), correlation coefficient (r), and combined accuracy (CA) for the three spatial interpolation techniques in estimating BF and BFI. These were the average values from three randomly selected times, as mentioned above. The higher performance of interpolation technique gives the lower MAE and RMSE, and r approaching 1.0 [36]. The combined accuracy (CA) [37] refers to the combination of RMSE, MAE, and r 2 , which equals 0.33 *(RMSE + MAE + (1 − r 2 )) and should be as small as possible (optimally 0).  Tables 3 and 4 show comparisons of the statistical indices for spatial interpolation techniques in estimating BF and BFI, respectively. The findings revealed that among the three spatial interpolation techniques, IDW gave the best performance. It provided the highest r-value and the lowest values of MAE, RMSE, and CA for both the calibration and validation stages. We can notice the bold number showing the best values of r, MAE, RMSE, and CA in Tables 3 and 4. Most of these indices pointed out IDW as giving the best performance, followed by kriging and spline, respectively. It was similar to the results found by Ly, Charles and Degré [28] for estimating monthly rainfall data. We further noticed that the use of the IDW and spline techniques overfitted in our data sets because the values of r, MAE, RMSE, and CA for the calibration stage were significantly different from those values in the validation stages. In contrast, the kriging technique gave less difference in the values between the two stages. We found that the kriging technique gave a lower r-value and higher MAE, RMSE, and CA values than the values provided by the spline technique for the calibration stage. However, it also mostly gave higher r-values and lower MAE, RMSE, and CA values than those given by the spline technique for the validation stage. It indicated that overfitting occurred for Spline technique for estimating BF and BFI in this experiment. This made the kriging technique more suitable or generalized than the spline technique for estimating BF and BFI in our study region. The reason why the computed BF by the LM and EF methods in the dry season for the validation stage gave better results (higher r-value, lower MAE, RMSE, and CA values) because less runoff fluctuation and less amount of runoff were noticed in the dry season for all runoff stations resulting in the data range was narrower as compared to those in the monsoon season. The data range in April used for the validation stage for the LM method was narrower than that for the EF method. As a result, the LM method's r, MAE, RMSE, and CA values were better than the EF method's. The data range in November used for the validation stage for the EF method was narrower than that for the LM method. As a result, the EK method's r, MAE, RMSE, and CA values were better than the LM method's.

Spatio-Temporal Variation of BF and BFI Using IDW Method
As the IDW method was the best spatial interpolation technique for our study, it was used to create the spatio-temporal variation of monthly and annual BF and BFI maps with a resolution of 30 × 30 m 2 . Figures 6 and 7 show the monthly and annual BF and BFI, calculated by the EF and LM methods, respectively. We can explain and discuss the following: the maximum average annual BF was found at the Mae Nam Tapi river basin, where the BF value varied in the range of 80 MCM to more than 100 MCM, followed by the Mae Nam Pattani river basin and the upper part of the Peninsular-West Coast river basin, where the BF value varied in the range of 80 MCM to 100 MCM. The minimum average annual BF was found in the Thale Sap Songkhla river basin, where the BF value was lower than 20 MCM. For the average annual BFI, we found that the maximum value ranging from 0.75 to 0.80 was located at the upper part of the Peninsular-West Coast River basin, followed by the Mae Nam Tapi river basin and the lower part of the Peninsular-West Coast river basin, where the BFI value varied in the range of 0.70 to 0.75. The minimum average annual BFI was found in the lower part of the Peninsular-East Coast river basin, where the BFI value was lower than 0.65. Considering the monthly BF, we found that the maximum BF value occurred between October and January, which is during the monsoon season. In addition, we found the Mae Nam Tapi river basin gave the maximum BF value during October and November, whereas the Mae Nam Pattani river basin gave the maximum BF value ranging from 80 MCM to more than 100 MCM. As naturally found in the dry season, the minimum BF value occurred between February and June, and most of the southern river basins had this values lower than 20 MCM. For the monthly BFI, the maximum BFI value happened in February, from the end of the monsoon season to the beginning of the dry season, for most of the southern river basin, where the BFI value varied in the range of 0.90 to 1.0. For the other months, we clearly found a difference in the high value of the BFI between the Peninsular-East Coast and Peninsular-West Coast areas because of the effect of the northeast monsoon and southwest monsoon winds, respectively. During April and September, a high value of BFI ranging between 0.70 and 0.90 was found in the Peninsular-East Coast area, that is, the Peninsular-East Coast, Mae Nam Tapi, Thale Sap Songkhla, and Mae Nam Pattani river basins. However, during October and January September, a high value of BFI ranging between 0.70 and 0.90 was found in the Peninsular-West Coast river basin. We also found that the minimum BFI value of less than 0.70 happened in March, which is in the dry season, for most of the southern river basins. Water 2021, 13, x FOR PEER REVIEW 13 of 16 Figure 6. Spatio-temporal variation of BF and BFI by using the EF method. Figure 6. Spatio-temporal variation of BF and BFI by using the EF method.

Conclusions
In this study, we mainly aimed to estimate the monthly and annual BF and BFI. Three spatial interpolation approaches (i.e., IDW, kriging, and spline) were applied to assess their efficacy. In addition, we compared two baseflow separation methods (namely, the LM and EF methods) for estimating BF and BFI. Our study area included 65 runoff gauged stations located in the southern river basins of Thailand. The findings revealed that the BF

Conclusions
In this study, we mainly aimed to estimate the monthly and annual BF and BFI. Three spatial interpolation approaches (i.e., IDW, kriging, and spline) were applied to assess their efficacy. In addition, we compared two baseflow separation methods (namely, the LM and EF methods) for estimating BF and BFI. Our study area included 65 runoff gauged stations located in the southern river basins of Thailand. The findings revealed that the BF and BFI, calculated using the EF method, had lower values than those calculated using the LM method for both annual and monthly periods. More than 75% of the basins in the southern part of Thailand had a BFI of more than 0.60, indicating quite considerable water source potential in this region. For our experiment, IDW provided the best performance among the three interpolation techniques, followed by kriging and spline, respectively. Finally, monthly and annual BF and BFI maps with a resolution of 30 × 30 m 2 were created in this research work. They provide valuable information for water management agencies in regard to planning and managing water resources in the southern river basin of Thailand.  Data Availability Statement: The datasets used and/or analyzed during the current study are available from the corresponding author upon reasonable request.