Aquifer classification and pneumatic diffusivity estimation using periodic groundwater level changes induced by barometric pressure

: To classify the aquifer type and to estimate the pneu‐ matic diffusivity using harmonic analysis and analytical expressions, we measured groundwater levels and baro‐ metric pressure in an unconfined aquifer with a thick un‐ saturated zone at a southern part of the Beppu volcanic fan, Japan. The groundwater level was inversely related to the barometric pressure, with little or no time lag behind baro‐ metric pressure. The groundwater level and barometric pressure exhibit periodic changes with prominent spectra at K 1 (lunar-solar diurnal) and S 2 (main solar semidiurnal) tidal constituents but not at O 1 (main lunar diurnal) and M 2 (main lunar semidiurnal) tidal constituents, which are the motions of solar and lunar relative to the earth. It indicates that the aquifer is an unconfined aquifer type with no earth tide effect. Pneumatic diffusivity in the unsaturated zone was estimated as 7.6 × 10 –2 to 5.7 × 10 –3 m 2 /s using changes between the groundwater level and the barometric pressure. The pneumatic diffusivity and the unsaturated zone thick‐ ness strongly affect amplitude ratios and phase lags between the groundwater level and the barometric pressure in the unconfined aquifer.


INTRODUCTION
Periodic changes of groundwater level in an aquifer include barometric pressures and earth tides (Allègre et al., 2016;Rahi and Hilihan, 2013). Changes in the barometric pressure and earth tides force dilatation or compression of aquifer materials, which is related to an increase or decrease in pore fluid pressure. Those changes induce deformation of the fluid-saturated rock (e.g. changes in volumetric strain), which engenders changes in the groundwater pressure in a well (Bredehoeft, 1967;McMillan et al., 2019). The changes in volumetric strain specifically appear in confined aquifers, but to a lesser degree in unconfined aquifers because the pore pressure balances barometric pressure (i.e. drained condition; Koizumi, 1994). In unconfined aquifers, the barometric pressure has greater influence on groundwater level changes than that of the earth tides (Bredehoeft, 1967;Rahi and Hilihan, 2013).
Correspondence to: Heejun Yang, Institute for Geothermal Sciences, Graduate School of Science, Kyoto University, 3088-176 Noguchibaru, Beppu, Oita 874-0903, Japan. E-mail: heejun@bep.vgs.kyoto-u.ac.jp An unconfined aquifer is formed below an unsaturated zone in the soil. The unsaturated zone includes air and some gases. The latter plays a crucially important role in considering unconfined groundwater. As the unsaturated zone thickens, a difference arises between groundwater pressure in wells and pore pressure in aquifers (Hosoya and Tokunaga, 2003;Weeks, 1979;Yusa, 1979). The difference appears as amplitude ratio changes and as phase lag delays between the barometric pressure and groundwater level fluctuations. Those differences are useful to infer the aquifer type classification as unconfined, semi-confined, or confined aquifers (McMillan et al., 2019;Rahi and Hilihan, 2013). Additionally, pneumatic diffusivity in the unsaturated zone can be estimated using harmonic analysis and analytical expressions (Rahi and Hilihan, 2013;Weeks, 1979).
Groundwater levels in an unconfined aquifer with a thick unsaturated zone at two observation wells were measured at the southern Beppu volcanic fan, Japan. Yusa (1979) studied the groundwater levels and found influences of the barometric pressure on the groundwater levels, but their periodicity and the exact aquifer type remain unclear. The present study was conducted to ascertain changes in groundwater levels associated with barometric pressure, and to estimate the possibility of being semi-confined aquifer and the pneumatic diffusivity using harmonic analysis and analytical expressions between the groundwater level and barometric pressure oscillations.

Study site
Beppu, situated on the eastern edge of the Beppu-Shimabara graben at the east coast of central Kyushu Island, Japan, is one of Japan's most famous hydrothermal areas ( Figure 1a). The hydrothermal activity spreads on a volcanic fan from the Garan-dake volcano to Beppu Bay, which is bounded by Kamegawa fault and Horita-Asamigawa faults with respective strikes of almost E-W orientation on the northern and southern edges of the fan (Yusa and Ohsawa, 2000). The fan deposits, which are about 50-150 m thick in the southern and northern areas, mainly comprise debris avalanches from Tsurumi-dake volcano and overlie middle Pleistocene volcanic lavas and tuff breccia (Kitaoka, 1990).
Observation wells of OW1 and OW2 are located at the southern area in the fan (Figure 1b). Their respective elevations are 73 m and 80 m above sea level ( Figure 1c). The respective depths of OW1 and OW2 are 60 m and 50 m, with a 11-m-long strainer in the unconfined aquifer that consists of a sand and gravel layer (Yusa et al., 1994), but the precise lithology at OW1 and OW2 remains unclear. Contours of the groundwater level decline from the mountains to the Beppu Bay: the groundwater flows from west to east. Transmissivity inferred as based on results of pumping tests in well OW1 was 6.0 × 10 -4 m 2 /s (Yusa et al., 1994).

Observations and analysis of groundwater level
Groundwater levels of the two wells were measured using loggers (CTD-diver, ±0.2% accuracy and 0.2 cm resolution; Daiki Rika Kogyo Co., Ltd.) installed at 3 m below groundwater level in the wells. They were monitored every 10 min. The barometric pressure logger (Baro-diver, ±0.5 cm accuracy and 0.2 cm resolution; Daiki Rika Kogyo Co., Ltd.) was set below the cap of well OW1. It was measured every 10 min.
Groundwater levels of the two wells were analyzed to verify the periodic oscillations during August 19, 2017 through July 23, 2018. The observation period was selected because the continuous data only exist within this period, which covers a hydrological year (i.e. the groundwater level rises and falls for a year). Tidal constituents of O 1 , K 1 , M 2 , and S 2 are the main diurnal and semidiurnal compo- . Cross-section map on the well OW1 and OW2 (c). p a stands for barometric pressure, p l denotes the soil gas pressure at water table in an aquifer, and Δh represents groundwater level changes in the wells. We designated the water level as the water table in the aquifer and groundwater level in the well nents of the ocean and earth tides, which are periodic as are the motions of solar and lunar relative to the earth. The K 1 (lunar-solar diurnal) and S 2 (main solar semidiurnal) tidal constituents have the same frequency as the barometric pressure variations. Consequently, they can be regarded as a signature for barometric pressure effects in observed groundwater levels, whereas the O 1 (main lunar diurnal) and M 2 (main lunar semidiurnal) tidal constituents explain a signature for earth tide effects (Rahi and Hilihan, 2013). We calculated the amplitude spectra of the groundwater levels and the barometric pressure to identify important tidal constituents using fast Fourier transform (FFT) with 1 h interval data of 8088 length. We specifically examined the periodic components in the groundwater level, so that the trend and long-term components in the groundwater level were removed using an FFT filter (Van Camp and Vauterin, 2005). The cut-off frequency was 0.03 cycle/h, which was chosen to extract semidiurnal and diurnal components of the groundwater levels and barometric pressure (i.e. high-pass filter; Yang et al., 2020). The FFT high-pass filter was applied using software (Origin Pro 2015; OriginLab Corp.). Then we analyzed the filtered data in the following manner to reduce noise. First, time-series data of groundwater levels and barometric pressure were subtracted by each averaged value. Then the subtracted data were divided by three segments (2696 data for one segment) before applying FFT. Secondly, zero padding was conducted before and after one segment to produce a power-of-two data structure, which became 4096 data for one segment. Finally, the amplitude spectrum of each timeseries was obtained by averaging the amplitudes of four segments on the same frequencies found from FFT.
Cross-correlation analysis was conducted to identify time lags between the groundwater levels and the barometric pressure using 10 min interval data. The cross-covariance function C xy (k) is given as (e.g. Larocque et al., 1998) where n stands for the time-series length, x − and y − respectively denote the means of the time-series of x t and y t , k represents the time lag, and r xy (k) is the cross-correlation coefficient. In addition, σ x and σ y respectively denote the standard deviations of the x t and y t time-series. A significant time lag was defined as one between the zero time lag (k = 0) and the maximum time lag indicating the maximum cross correlation coefficient (Yang et al., 2020). We designated the time-series of the barometric pressure and the groundwater level respectively as x − and y − . Pneumatic diffusivity can be estimated from changes between barometric pressure (p a ) and groundwater level in wells (Δh) (Weeks, 1979;Yusa, 1969). Soil gas in the unsaturated zone flows in isothermal and vertical directions (Carslaw and Jaeger, 1959). The amplitude and phase lag of soil gas pressure (p l ) change at any depth in the unsaturated zone. We used the groundwater level in the well and the water table in the aquifer. Periodic vertical flow of soil gas from the land surface to the water table is governed by a simple diffusion equation (Weeks, 1979) as subject to the following boundary conditions.
Therein, α represents the pneumatic diffusivity (m 2 /s), p stands for barometric pressure (m), τ denotes the period of barometric pressure (s), z denotes the depth (m), t is time (s), l expresses the distance from water table to land surface (m), and A 0 signifies the amplitude of the barometric pressure (m). From Equations (3), (4), and (5), a solution of the soil gas pressure at z = l is expressed (Yusa, 1969) as shown below.
In those equations, φ expresses the phase lag of the soil gas pressure at the water table behind that at the land surface. The groundwater level change is expressed as Δh = p ap l (Yusa, 1969). Therefore, the following expressions are obtained from Equations (4 to 9) (Yusa, 1969).
Therein, θ represents the phase lag between the groundwater level in wells and the barometric fluctuation in the land surface. For this study, pneumatic diffusivity was estimated using the amplitude ratio (Equation 12) and phase lag (Equation 11), respectively, by fitting the type curves of Δh/p a vs. kl and θ vs. kl. The value of kl includes information of the pneumatic diffusivity, unsaturated zone thickness, and periodicity. The amplitude ratios obtained using the amplitude spectrum at K 1 tidal constitute the data. The phase lag (θ) was used with the time lag obtained using the cross correlation analysis. The parameters were l = 32 m at well OW1 and 40 m at well OW2, τ = 1 day, and α = 200-500,000 m 2 / day.

RESULTS
Groundwater levels have seasonal variation with swelling in November-April (Figure 2a), indicating that they are affected by rainfall infiltration (Yusa, 1979). The average groundwater levels of wells OW1 and OW2 were, respectively, 41.1 m and 41.3 m above sea level. The groundwater level in the unconfined aquifer is typically found a few meters below the land surface, but those of the wells were measured at about 35 m below the land surface. They have a thick unsaturated zone. The filtered levels of OW1 and OW2 (< 1 day period) were found to be related to the filtered barometric pressure. The mean ranges of OW1, OW2, and barometric pressure which were converted into cmH 2 O from hPa were, respectively, 7 cm, 10 cm, and 10 cm (Figure 2b).
The FFT high-pass filtered data (cut-off frequency is 0.03 cycle/h), periodic components of the groundwater level, and the barometric pressure were used for FFT and cross-correlation analysis (Figure 3). The amplitude spectra of the groundwater levels and the barometric pressure show two strong signals at frequencies of 0.0417 and 0.0833 cycle/h, which are equal to the K 1 and S 2 tidal constituents (Figure 3a). Those constituents are typically associated with heating and cooling of the atmosphere at the ground surface caused by solar radiation (Spane, 2002). The amplitude spectra of O 1 and M 2 tidal constituents are lower than those of K 1 and S 2 tidal constituents. Positive time lags indicate that the groundwater levels are delayed from the barometric pressure. The amplitude ratio was calculated using the K 1 tidal constituent from the FFT result because the long-period wave propagates farther than the shorter one does (e.g. Yang et al., 2015). The amplitude spectrum was 0.62 cm at OW1, 1.99 cm at OW2, and 1.73 cm at the barometric pressure. Consequently, the amplitude ratio is 0.35 at OW1 and 1.13 at OW2. The time lag at the minimum coefficients showed a 1.3-h delay (39°p hase lag) at OW1 and 0-h delay at OW2 (Figure 3b). The minimum cross-correlation coefficient of the groundwater levels showed negative values, reflecting an inverse relation with the barometric pressure. Type curves of the amplitude ratio of Δh/p a vs. kl ( Figure  4a) and the phase lag θ vs. kl (Figure 4b) are drawn with three parameters: pneumatic diffusivity (α) = 200-500,000 m 2 /day); aquifer thickness (l) = 10, 32, 40, and 80 m; and periodicity = 1 day. The type curve finds the most appropriate pneumatic diffusivity and aquifer thickness. The kl values calculated from the amplitude ratio and phase lag were, respectively, 0.7 and 2.2 at OW1, and 1.9 and 3.2 at OW2. The kl values were higher at OW2 than at OW1. The pneumatic diffusivity values, calculated using α = l 2 π/τ (kl) 2 based on kl values, were, respectively 7.7 × 10 -3 and 7.6 × 10 -2 m 2 /s at OW1, and 5.7 × 10 -3 and 1.6 × 10 -2 m 2 /s at OW2. Although the phase lag shows positive radians (Equation 10), the pseudo-advance wave resulted from the calculation process (Hosoya and Tokunaga, 2003;Weeks, 1979). Furthermore, we found that the groundwater level was delayed from barometric pressure using crosscorrelation analysis. Therefore, we consider the groundwater level delays. Figure 3. Amplitude spectra of the groundwater levels and the barometric pressure (a) resulting from the FFT. Four harmonic components (O 1 , K 1 , M 2 , and S 2 ) are marked with gray dotted lines, with data given in cycles/h (cph). Crosscorrelation analysis of the groundwater levels are shown as outputs against the barometric pressure as an input for the whole period with 10 min interval data (b)

DISCUSSION
The groundwater levels of OW1 and OW2 have periodic fluctuations of K 1 and S 2 tidal constituents. They respectively change after 1.3 h and 0 h behind barometric pressure. In the case of an unconfined aquifer, the delay change of groundwater level derives from the difference between barometric change in the well and pressure propagation through the unsaturated sediment layer. This lag is expected to generate a barometric signature of K 1 and S 2 constituents (Weeks, 1979). The spectra of O 1 and M 2 constituents in the wells were not found to be significant, indicating no confinement (Rahi and Hilihan, 2013), and meaning that both wells are the unconfined aquifer type.
The amplitude ratio and phase lag were, respectively, high at OW2 (1.13 and 0°) and low at OW1 (0.35 and 39°). Therefore, the pneumatic diffusivity at OW2 is lower than that at OW1. As pointed out by Weeks (1979), the wide range of the pneumatic diffusivity estimated from the amplitude ratio and phase lag makes it difficult to ascertain the pneumatic diffusivity accurately, but it is readily apparent that the pneumatic diffusivity values of the two wells differ. Weeks (1979) reported that, in an unconfined aquifer system, low pneumatic diffusivity and thick unsaturated Figure 4. Graph and type curves showing relations between the amplitude ratio (a) and phase lag (b) of the groundwater level fluctuations in an unconfined aquifer to barometric fluctuations as a function of kl that includes information of the pneumatic diffusivity, unsaturated zone thickness, and periodicity. The grey numbers on the left axis of the amplitude ratio and phase lag are values calculated using FFT and cross-correlation analysis zone cause the zero phase lag between the groundwater level and barometric pressure. In addition, Equation (9) shows that the kl value increases as the pneumatic diffusivity decreases. It also shows that the unsaturated zone thickness increases, which indicates that the propagation velocity of the barometric pressure from land surface toward the water table in the unsaturated zone is low at OW2 because of the low pneumatic diffusivity. Moreover, it shows the thicker unsaturated zone that derives from the fact that the unsaturated zone thickness of OW2 is greater by about 8 m than that of OW1. Figure 5 shows how the pneumatic diffusivity and unsaturated zone thickness affect the amplitude ratio and phase lag. The figure was drawn using two constant pneumatic diffusivities: 4.2 × 10 -2 m 2 /s ( Figure 5a) and 1.1 × 10 -2 m 2 /s (Figure 5b), which respectively represent the average values of OW1 and OW2, and unsaturated zone thicknesses of 1-120 m. The amplitude ratio and the phase lag increase and decrease, respectively, as the unsaturated zone grows, but the high amplitude ratio and the zero phase lag appear at the thinner unsaturated zone thickness when the pneumatic diffusivity is lower (Figure 5b). These results demonstrate that the low pneumatic diffusivity also causes less propagation velocity of the barometric pressure. Consequently, if the soil gas pressure changes could be ignored near the water table, then the phase lag between the barometric pressure and groundwater level becomes zero and the amplitude ratio would become high (Figure 5b), which well explains the results of the amplitude ratio and the phase lag at OW2. Therefore, we infer that the low pneumatic diffusivity and the thicker unsaturated zone at OW2 would engender to the high amplitude ratio and the zero phase lag.
It is noteworthy that the groundwater flow in or out of the well must be considered. The flow depends on the permeability and storage coefficient in aquifers (McMillan et al., 2019;Rojstaczer, 1988). Although differences of the permeability and storage coefficient of between OW1 and OW2 in the unconfined aquifer were assumed as negligible because two wells are situated on the same aquifer and are mutually close (105 m), the results of this study imply that the pneumatic diffusivity between OW1 and OW2 differs, indicating that the permeability and the storage coefficient might be different. Additionally, Yusa (1969) reported that the hornblende andesite rocks distributed at 5 m depth below the land surface with 40 cm thickness at OW1 disturb the passage of air and water. Consequently, one might assume from the result of the low pneumatic diffusivity that the rock sedimentation at OW2 is thicker. To elucidate differences of the amplitude ratio and the phase lag between the wells, further analysis such as a numerical modeling is necessary, including that of groundwater flow in the unconfined aquifer and the existence of a low pneumatic diffusivity layer in the unsaturated zone. In addition, the lithology survey might be necessary to provide solid evidence of the existence of sediments showing low pneumatic diffusivity at OW2.

CONCLUSION
Periodic changes in groundwater levels associated with the barometric pressure in the unconfined aquifer with thick unsaturated zone were analyzed at the southern Beppu volcanic fan, Japan. The aquifer type and pneumatic diffusivity were estimated using harmonic analysis and analytical expression.
Groundwater levels showed seasonal variation: high in November-April and low in May-October. Significant spectra were found at K 1 and S 2 constituents but not at the O 1 and M 2 constituents, reflecting the unconfined aquifer type. The cross-correlation coefficient, which was used to calculate the barometric pressure as an input and the groundwater level as an output, fluctuated twice each day. The time lags at the maximum coefficients showed a 1.3 h delay at OW1 and a 0 h delay at OW2 with the inverse relation. The pneumatic diffusivity values which were found using the analytical expression were, respectively, 7.7 × 10 −3 and 7.6 × 10 −2 m 2 /s at OW1 and 5.7 × 10 −3 and 1.6 × 10 −2 m 2 /s at OW2. They were lower at OW2 than at OW1.
The groundwater level at OW1 and OW2 oscillates by inducement from barometric pressure changes. The amplitude ratio and phase lag in OW1 and OW2 are controlled by the pneumatic diffusivity and thickness of the unsaturated zone. Results also show that pneumatic diffusivity is applicable as a parameter for groundwater simulations as a fully coupled surface and subsurface flow calculator.