Taiwan is bisected by The Tropic of Cancer, where its northern and central regions are subtropical while the southern regions are tropical. This distinct condition leads the weather to vary considerably across the island (Chen et al., 2021; Dibaj et al., 2020; Hung et al., 2017). The accumulated precipitation in the study area during 2018 is shown in Figure 1. It fluctuated during the dry and wet seasons with the highest and the lowest accumulated precipitation in August and December with 449.5 mm (the wet season) and 2.0 mm (the dry season), respectively.

The study area is located in the Taichung-Nantou Basin along the Wu River in Central Taiwan (Figure 2). This basin covers parts of Taichung city, Nantou County, and Changhua County. The study area encompasses the Maoluo River in the south, the Touke Hill lands of Nantou in the east, the southern district of Taichung City in the north, and the Bagua Plateau in the west. The average ground elevation ranges from 30 to 65 m above sea level, with the highest elevation found in the southeastern region and gradually decreasing towards the northwest of the study area.

Furthermore, the geological setting of the study area is mostly composite by the alluvial deposits, in which the western side of the survey area faces Bagua Hill, which includes the Toukoshan formation and the terrace deposits, as well as the Pakushan Anticline. The Tokushan formation is classified into three divisions: the lower division (up to 900 m thick) consists of sandstone and shale with a pebbly horizon, the middle division (50–100 m) consists of an alternating bed of sand, clay, and gravel that containing fresh-water, the upper-division consists of several hundred meters of massive conglomerate with a few thin beds of crudely cross-bedded sandstone. The terrace deposits consist mainly of unconsolidated gravel with flat-lying sandy or silty lenses. The west side of the survey area is Touke Hill, where the Chelungpu Fault is found along the edge of the hill. The area is also traversed by two rivers, the Maoluo river and the Wu River, whose branches confluence on the northwest right. The alluvial deposits mainly consist of thick gravel and sand with good permeability, making the area highly potential for groundwater recharge (Ho and Chen, 2000).

To investigate the groundwater condition, we conducted ten Electrical Resistivity Imaging (ERI) survey lines (shown in red in Figure 3) around three observation wells (shown as magenta circles). In order to analyze the influence of seasonal conditions on the trends of ERI results, we performed time-lapse surveys at consistent survey sites during both the dry and wet seasons. A comprehensive dataset spanning five distinct periods was collected, encompassing February (dry season), May (transition from dry to wet seasons), July (wet season), September (wet season), and October (transition from wet to dry seasons). This methodological approach enabled us to meticulously examine and analyze the discernible patterns exhibited by the ERI results across diverse seasonal contexts.

The fundamental idea of the resistivity method involves injecting current into the subsurface using a pair of electrodes. This current generates a potential difference in the subsurface, which is measured and monitored by another pair of electrodes, resulting in the measurement of the apparent resistivity. It’s important to note that the apparent resistivity represents the collective measurements of the electrical properties of all the subsurface layers within the corresponding electrode configuration. Therefore, an additional inversion process is required to determine the true resistivity value, which can provide valuable information about the subsurface lithology or materials. Figure 4 summarizes the workflow of the present study.

This study utilized the 4-point light 10 W resistivity meter and the active electrode (ActEle) system (Lippmann Geophysical Instruments) (Lippman, 2014) in Wenner-Schlumberger configuration with a 1 m electrode spacing for the field data measurements. Wenner-Schlumberger configuration is less susceptible to noise than the other arrays, which means that it yields a high signal-to-noise ratio with better sensitivity to horizontal structure (Dahlin and Zhou, 2004). The raw data was later inverted with the EarthImager2D™ version 2.4.2.627 (Advanced Geosciences, 2006). For a detailed review of the resistivity method inversion techniques, we forward readers to Sharma and Verma (2015).

After the inverted resistivity section was obtained, we carried out the groundwater level estimation of the study area. First, we calculated the water contents in the five vertical profiles in the central part of each resistivity survey line for the estimation of water content by adopting a similar process to Dietrich et al. (2018). We may deduce from Archie’s Law (Archie, 1942) how formation resistivity relates to porosity, saturation, and pore water resistivity in a clay-free matrix: ρ = α ρ w ϕ − n S w − m (1) where ρ is the formation resistivity, ρ w is the pore water resistivity, ϕ is the porosity, S w is the saturation, and α , m, and n are empirical parameters. Parameters m and n are usually termed saturation exponent and cementation factor, respectively. For general homogeneous rocks and soils, m ranges from 1.8 to 2.2, and n is about 2; so, m = n ≅ 2. On the other hand, it should be noted that the empirical factor α is not included in its original form of Archie’s law [Detailed in Glover. (2015)].

As it is known, the volumetric water content ( θ ) at the saturation point is equivalent to the total soil pore space, expressed as θ = ϕ − n S w − m (2) considering the value of m and n, and substituting Eq. 2 into Eq. 1, ρ f will yield ρ = ρ w θ − 2 (3)

If we assume that the resistivity value varies only with the water saturation, we can calculate the normalized relative saturation (S _r ) at different depths in the vadose zone as S r = θ u θ s = ρ s ρ u (4) where θ u , θ s , ρ u , and ρ s are unsaturated layer volume water content, saturated layer volume water content, unsaturated layer formation resistivity, and saturated layer formation resistivity, respectively. In this case, the lowest resistivity value obtained from the ERI of each site was set as ρ s , and the data of the topmost 2 m are eliminated because it represents the properties of the surface soil layer instead of the properties of the gravel layer. Next, we can obtain the normalized volumetric water content by multiplying S r and ϕ A (Average Porosity, 0.26), expressed as θ = S r ϕ A (5)

Figure 5 depicts how the normalized volumetric water content varies from the ground surface to the saturated layer in an unconfined aquifer, as determined from selected resistivity measurements.

We observed in Figure 5 that the vertical change in the water content shows a similar tendency to the Soil-Water Characteristic Curve (SWCC) obtained in lab tests, e.g., (Chin et al., 2010; Fredlund et al., 2011; Habasimbi and Nishimura, 2018; Azmi et al., 2019; Eyo et al., 2022). Hence, if we assume that the suction head is linearly proportional to the elevation of the groundwater level in the unconfined aquifer as discussed in Krahn and Fredlund (1972) and Stephens (2018), we then can estimate the groundwater level quantitatively and the relative hydraulic parameters in the vadose zone with the SWCC model. We applied the Van Genuchten (VG) model (Genuchten, 1980) that describes the physical relationship between the water contents and suction in the vadose zone through the following equation: θ h = θ r + θ s − θ r 1 + α h n m ;   w h e r e ;   m = 1 − 1 / n (6) Eq. 6 has four independent parameters, such as θ s , θ r , α , and n that need to be estimated from observed SWCC. Where θ s is the saturated water content (L³/L³); θ r is the residual water content (L³/L³); α is a parameter related to the air entry pressure (L⁻¹), and n is a dimensionless parameter related to the pore-size distribution of soil and h is the suction head (L). Furthermore, with the VG model of SWCC, we were able to estimate the air entry suction head relative to the same base (h _a ), if we designate a saturated base with a depth of (H _s ) and assume that the section head is proportional to the height of the presumed saturated base. Thus, it is possible to calculate the depth of the groundwater table (D): D = H s − h a (7)

We could invert the VG parameters and the air entry suction or heights from the base depth by minimizing the root mean square differences between the estimated and measured water content. To optimize the minimum object equation in the inversion, we utilized the Generalized Reduced Gradient (GRG) nonlinear program. Once we determined the groundwater depth, we advanced to estimate the groundwater level (GWL) by subtracting D from the ground level (G _L ), expressed as: G W L = G l − D (8)

In addition, we can obtain the theoretical specific yield, S _y (Eq. 9) by subtracting the residual water content from the saturated water content and the specific yield capacity, S _c (Eq. 10) with the VG model in Eq. 6, as: S y = θ s − θ r (9) S c = 1 D ∫ 0 D θ s − θ h d D ≅ 1 D ∑ 0 D θ s − θ h ∆ D (10) where θ h is the volumetric water content at a different depth, h, and ∆ D is the incremental depth (Chang et al., 2022). We will discuss the S _y and S _c later in the discussion.

Time series clustering (TSC) is an unsupervised machine-learning technique for identifying and grouping similar objects or data points into structures that are easier to understand. This technique is useful and capable of enhancing the interpretation and provides an intuitive explanation of relevant aspects of given time series datasets. Accordingly, it may lead us to discover interesting patterns that can be either frequent or rare patterns (Aghabozorgi et al., 2015; Rodriguez et al., 2019). A primary benefit associated with a TSC solution is its ability to automatically recover from failures, thereby eliminating the need for user intervention, compared to classification, the TSC not required prior information as label data. However, TSC entails certain drawbacks, namely, its inherent complexity and the inability to recover from database corruption. To address this limitation, the PCA (Principal Component Analysis) technique is utilized to reduce complexity, and preliminary data preparation steps are applied to prevent database corruption, such as removing blank data. The outcome of this method is the generation of a clustered resistivity value, which is determined based on its similarity to other values. This clustering allows for the identification and analysis of noteworthy patterns.

In this study, we attempted to use Hierarchical Agglomerative Clustering (HAC) by using the Orange data mining toolbox in python (Demsa et al., 2013). The toolbox effectively aids our aim to group synchronous and linearly correlated series or similarities in time. We were able to accomplish two tasks by utilizing such a technique: first, we managed to explain the trend of each site based on the VG model parameters over time in the area, and second, to describe the variation of groundwater level changes over time in the area. Nevertheless, TSC can be challenging due to the enormous complexity or dimensionality of the time series datasets. For instance, the second task dataset has a higher dimension (7200 × 4) compared to the first task dataset (50 × 3), therefore dimensionality reduction is necessary to conduct for the second task before applying the HAC, i.e., Principal Component Analysis (PCA ).

Five period of ERI surveys were carried out at ten sites in the Taichung-Nantou Basin along the Wu River in 2018. Figure 6 shows the inverted resistivity images for the ten sites, and Figure 7 shows the time-lapse resistivity images collected at the WS05 site.

WS01 to WS04 sites are located on the northern part of the Wu River, while WS05 to WS10 sites are located on the southern part of the Wu River and the northern part of the Maoluo River, as shown in Figure 3. In general, the resistivity values vary in the range of 40 to 1,400 Ωm. The topmost 2 m layer with low resistivity anomaly from 40 to 100 Ωm (bluish color) represents the soil layer of the paddy field. Thus, we excluded this layer in order to estimate the normalized water content (see Section 2.3). Below this layer, the resistivity increases to over 400 Ωm until a certain depth (yellowish to reddish colors) that mostly lies around 8 m depth. At some point, the resistivity decreases from the peak of over 1,000 Ωm to a level between 200 and 300 Ωm. Such trend may reflect the effect of the increasing water content from the lower vadose zone to the saturation zone (transition zone). This resistivity trend is detected in every survey line. This implies that the study area lies on the same lithological formation, i.e., Sand, Clay and Gravel (see Section 2.1). Interestingly, the inverted resistivity images located nearby the confluence of two rivers (WS01 and WS05) are dominated by blue and green contour colors (low resistivity), while the WS04 and WS10 located on the upper course are dominated by yellow and red contour colors (high resistivity).

The inverted time-lapse resistivity images of WS05 collected in 2018 reveal a significant resistivity variation in the vadose zone between the dry and the wet seasons. The region with a high resistivity value (>500 Ωm) shrank from February to July indicating the water content in the vadose zone increased due to an increase in rainfall recharge. Meanwhile, the region with a low resistivity value (<150 Ωm) expanded from July to October owing to limited rainfall. The inverted resistivity image on July (wet season) shows the lowest resistivity distribution compared to others since the accumulated precipitation was the highest at about 334 m·mm/month (see Section 2.1). To observe the significant difference in the resistivity distribution, the resistivity values for each month were subtracted from the resistivity value of the dry season (February), and the results are shown in Figure 8. These results highlighted that the resistivity value plummeted from February to July (>50%) especially in the vadose zone owing to significant accumulated precipitation. On the other hand, May-February only shows slight resistivity changes up to 25% in the top-right corner of the section. The same trend is also observed in September-February and October-February with resistivity change of up to 30% and 15%, respectively. The resistivity variations in the dry and wet seasons can be linked to changes of the water content in the vadose zone, as was what stated in Archie’s law in Eq. 1. Thus, we were able to estimate further the variation of water contents and the groundwater table with the inverted resistivity images collected in different seasons.

Figure 9 presents the fitted VG model for different months at the survey sites in the study area. Generally, the fitted curves during the wet season (blue, green and purple dashed curves) have higher air entry heights than those in the dry season (red and orange curves). Table 1 presents an example of the fitted parameters of the VG curve at the WS05. a varies from 7.33 to 10.46 m, n varies in the range of 1.52–3.12, and θ r is about 0.05–0.1. In addition, we assumed that the saturated water content would be equal to the average porosity, 0.26, for each measurement.

Table 2 lists the estimated groundwater levels (GWL) measured at different survey sites during the observation period in 2018. WS01-WS10 are the resistivity survey lines, while 98GH-06, 98GH-09, and 98GH-11 are the available observation wells around the survey area. The GWL covers a depth range of 22.9–51.4 m. The majority of the GWL are gradually increasing both in the observation wells and resistivity measurements from February to July due to an increase in recharge from the rainfall, in contrast to the GWL from July to October, which decreased owing to a decrease in recharge from the rainfall. To observe the spatial distribution trend of the GWL, the GWL were plotted onto the map. Figures 10A–E shows the spatial distribution of GWL for February, May, July, September, and October, respectively. The blue contour color stands for the high GWL, while the red contour color represents the low GWL.

Aside from seasonal variations in groundwater levels, it is also vital to understand the water storage capacity of the region for potential groundwater reservoirs. Using Eqs 9, 10, we were able to estimate the theoretical specific yield, S _y , and the specific yield capacity, S _c . The results are listed in Table 3. In this case, the maximum and average for each estimation from all measurements time were calculated. The maximum and average of the theoretical specific yield as well as those of the specific yield capacity are in the range of 0.16–0.25, 0.14–0.24, 0.06–0.24, and 0.04–0.17, respectively. Figure 11 shows the spatial distribution of these estimations.