Introduction: Hydrogeological Settings of the Beppu Geothermal System
One of the world’s most active geothermal systems resides beneath the city of Beppu. Geothermal manifestations, such as steaming ground and active hot springs, are widely distributed, and both the spring discharge rate and the number of active vents are among the highest in Japan. However, the origin and flow dynamics of the shallow meteoric groundwater, as well as the mechanisms governing its long-term decline, remain critical hydrogeological questions.
In this series, we analyze time-series data of unconfined groundwater levels measured at the Beppu Geothermal Research Center of Kyoto University. Using signal processing, statistical analysis, and Python-based scripts, we characterize the sub-surface flow dynamics.
The Beppu Geothermal System and Hydrothermal Flow Paths
The subsurface geological structure and the conceptual hydrogeological model of the hydrothermal flow system in southern Beppu are illustrated in Figure 1.
The subsurface flow system in southern Beppu is categorized into three major hydrothermal fluid types distributed across different depth intervals.
| Type | Water Quality / Facies | Depth Range (relative to sea level) | Key Characteristics |
|---|---|---|---|
| Type A | Ca・Mg-HCO₃ type | −100 to −200 m | Shallow meteoric water origin. Strong recharge signal. |
| Type B | Na-Cl type | approx. −250 m | Located south of the Horita Fault. Enriched in deep thermal fluid. |
| Type C | Na-HCO₃ type | approx. −300 m | Deepest zone. Directly connected to the hydrothermal reservoir. |
In addition to these types, Type D (steam-heated acid sulfate water) is present near the ground surface. This acidic water is generated by the reaction between volcanic gases and shallow groundwater, contributing to the famous “Hell Hot Springs” (Jigoku) in Beppu.
The geological framework of this aquifer system consists of volcanic products (tuff and tuff breccia) from Mt. Tsurumi and Mt. Garandake, overlain by unconsolidated sediments (sand and gravel). The Horita and Kamegawa faults act as low-permeability boundaries (barriers) that compartmentalize the groundwater flow, leading to significant variations in water quality and temperature across the fault blocks.
Mechanisms of Long-Term Water Level Decline in Unconfined Aquifers
Long-term hydrological monitoring (Yusa, 1979; 2001) has revealed a steady decline in the unconfined groundwater levels in this region. This downward trend is primarily attributed to two distinct anthropogenic and environmental factors:
Evaluating and distinguishing between these two mechanisms requires a rigorous analysis of the periodic fluctuations in the groundwater level. In this study, we utilize Fast Fourier Transform (FFT) and cross-correlation analysis to analyze these signals.
Hydrogeological Significance of Unconfined Groundwater Level Monitoring
Monitoring the shallow unconfined groundwater table, rather than directly measuring the hydraulic head of deep geothermal reservoirs, offers several distinct advantages:
- Proxy for Deep Reservoir Pressure: The unconfined aquifer and the deep geothermal reservoirs are hydraulically connected in the vertical direction. Consequently, pressure variations in the deep reservoirs manifest as water table fluctuations in the shallow zone.
- Evaluation of Aquifer Hydraulic Properties: Analyzing the attenuation and lag of water table responses allows for the estimation of aquifer hydraulic diffusivity (defined as transmissivity divided by storativity).
- Feasibility of Long-Term Monitoring: Unlike deep geothermal fluids, which present harsh conditions (high temperature, high pressure, and mineral scaling) that damage sensors, the shallow unconfined aquifer can be safely monitored over multi-year periods using standard CTD-Diver pressure loggers.
Monitoring Setup and Dataset Characteristics
The time-series dataset was collected from two observation wells located within the campus of the Beppu Geothermal Research Center, Kyoto University (Noguchimoto-machi, Beppu).
| Parameter | Observation Well 1 (OW1) | Observation Well 2 (OW2) |
|---|---|---|
| Ground Elevation | 73 m a.s.l. | 80 m a.s.l. |
| Well Distance | 110 m | |
| Mean Water Level (Elevation) | 41.06 m a.s.l. | 41.25 m a.s.l. |
| Monitoring Period | August 19, 2017 – July 23, 2018 | |
| Sampling Interval | 1 hour ($N = 8,136$ data points) | |
| Instrumentation | CTD-Diver (Resolution: 2 mm) | |
| Monitored Variables | Groundwater Level, Barometric Pressure, Precipitation | |
Although the two observation wells are separated by approximately 110 m, their mean water levels are nearly identical (approximately 41 m a.s.l.), indicating a relatively flat water table in this unconfined aquifer. However, high-resolution analysis reveals distinct differences in their response to barometric pressure variations, which is discussed in Part 4 of this series.
First Step in Data Analysis: Visualization with Python
We first implement a Python script to import and visualize the time-series data.
To ensure reproducibility and ease of learning, we initially utilize a synthetic (dummy) dataset generated using numpy. When applying this analysis to field measurements, the synthetic generator section can be replaced with pd.read_csv() as described in the subsequent section.
# [Synthetic Data] Python script for generating and plotting synthetic groundwater level and barometric pressure data
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
# Plot settings
plt.rcParams.update({
"font.family": "sans-serif",
"font.sans-serif": ["Arial", "DejaVu Sans"],
"axes.unicode_minus": False,
"figure.dpi": 150,
})
# ---- Generate Synthetic Dataset (Replace with pd.read_csv() for field data) ----
np.random.seed(42)
n = 8136 # Hourly sampling for ~1 year
t = pd.date_range("2017-08-19", periods=n, freq="h")
# Define seasonal trend, barometric response, and random noise
seasonal = 2.0 * np.sin(2 * np.pi * np.arange(n) / (24 * 365))
barometric = 0.03 * np.sin(2 * np.pi * np.arange(n) / 24)
noise_ow1 = 0.005 * np.random.randn(n)
noise_ow2 = 0.020 * np.random.randn(n)
gwl_ow1 = 41.06 + seasonal + barometric + noise_ow1
gwl_ow2 = 41.25 + seasonal + 1.05 * barometric + noise_ow2 # OW2 exhibits larger response amplitude
baro = 1010 + 8 * np.sin(2 * np.pi * np.arange(n) / (24 * 14)) \
+ 3 * np.random.randn(n)
# Figure configuration
fig, axes = plt.subplots(2, 1, figsize=(13, 7), sharex=True)
fig.patch.set_facecolor("#F8F9FA")
# Upper panel: Groundwater level
ax1 = axes[0]
ax1.set_facecolor("#F0F4F8")
ax1.plot(t, gwl_ow1, lw=0.8, color="#2563EB", label="Well 1 (OW1)", alpha=0.9)
ax1.plot(t, gwl_ow2, lw=0.8, color="#DC2626", label="Well 2 (OW2)", alpha=0.9)
ax1.set_ylabel("Groundwater Level (m a.s.l.)", fontsize=10)
ax1.legend(fontsize=9, loc="upper right")
ax1.grid(True, ls="--", lw=0.4, color="#CBD5E1")
ax1.set_title("Unconfined Groundwater Levels in Southern Beppu (Synthetic: 2017–2018)", fontsize=12)
# Lower panel: Barometric pressure
ax2 = axes[1]
ax2.set_facecolor("#F0F4F8")
ax2.plot(t, baro, lw=0.8, color="#374151", label="Barometric Pressure", alpha=0.85)
ax2.set_ylabel("Barometric Pressure (hPa)", fontsize=10)
ax2.set_xlabel("Date", fontsize=10)
ax2.legend(fontsize=9, loc="upper right")
ax2.grid(True, ls="--", lw=0.4, color="#CBD5E1")
ax2.xaxis.set_major_formatter(mdates.DateFormatter("%m/%d\n%Y"))
plt.tight_layout()
plt.savefig("beppu_gwl_overview-en.png", bbox_inches="tight")
plt.show()Insights from the Synthetic Dataset
Visual inspection of the raw synthetic time-series plot (Figure 2) reveals several distinct features. Specifically, transient high-frequency fluctuations are superimposed on the long-term seasonal trend.
Analysis of Field Measurements
In this section, we transition from the numerical simulation using synthetic data to the analysis of actual field measurements.
The dataset analyzed herein consists of unconfined groundwater levels, barometric pressure, and precipitation continuously monitored at hourly intervals for approximately one year at the observation wells of the Beppu Geothermal Research Center, Kyoto University. The field measurement dataset is available for download below:
- Field Measurement Data: BP.csv
To adapt the previously defined Python script for the field measurements, the synthetic data generation block (lines 207–224) should be replaced with the following data import routing utilizing pandas.read_csv:
# ---- Load Field Measurements ----
df = pd.read_csv("BP.csv", parse_dates=["datetime"], index_col="datetime")
gwl_ow1 = df["OW1"]
gwl_ow2 = df["OW2"]
baro = df["press"]Note: The variables are assigned corresponding to the respective column headers (OW1, OW2, and press) in the CSV file.
Applying this modification and executing the script yields the field data plot presented in Figure 3.
The field measurement time-series (Figure 3) reveals several complex physical dynamics:
- Baseline Water Levels and Local Gradients: The mean water levels of OW1 and OW2 fluctuate within the range of 39.5 to 40.0 m a.s.l. The water level in OW2 remains consistently higher than in OW1 throughout the year, reflecting a stable local hydraulic gradient and indicating the direction of groundwater flow.
- Barometric Response during Typhoon Events: The barometric pressure drops rapidly (below 985 hPa) during storm events, such as mid-September typhoons. Simultaneously, the groundwater levels in both wells show a sharp increase. This immediate, in-phase rise in water level in response to a decline in barometric pressure is a classic barometric response. It demonstrates that the unconfined aquifer behaves elastically under changing atmospheric load.
- Long-Term Seasonal Trends: Similar to the synthetic dataset, a long-term seasonal fluctuation with a period of several months is observed. This represents the catchment-scale precipitation-recharge and evapotranspiration cycle.
- High-Frequency Micro-Fluctuations: Continuous micro-fluctuations with amplitudes of 1 to 5 cm are superimposed on the baseline. These represent the combined effects of earth tides and diurnal/semidiurnal barometric variations.
Although the 1 to 5 cm amplitude of the short-period fluctuations is extremely small, these micro-scale signals encapsulate crucial hydrogeological information, including aquifer hydraulic properties, atmospheric loading response, and tidal forcing (oceanic and earth tides).
To establish a baseline for analyzing these short-period fluctuations, the relationship between ocean tide (sea level) and its constituent frequency components is presented in Figure 4.
The blue curve in the left panel of Figure 4 represents a 30-day time-series of ocean tide measurements. The sea level fluctuates periodically at diurnal (approx. 24 h) and semidiurnal (approx. 12 h) scales. The amplitude of these fluctuations is modulated over a fortnightly (approx. 15-day) cycle corresponding to the spring-neap tidal variations.
Applying the FFT to this complex time-series transforms the data into the frequency domain, yielding the power spectrum shown in the right panel. This transformation isolates sharp spectral peaks corresponding to specific tidal constituents:
- \(M_2\) (Principal lunar semidiurnal constituent, period of 12.42 h): A semidiurnal tide induced by lunar gravitational forcing, typically exhibiting the highest spectral energy in most marine tidal records.
- \(S_2\) (Principal solar semidiurnal constituent, period of 12.00 h): A semidiurnal tide induced by solar gravitational forcing.
- \(K_1\) (Luni-solar declinational diurnal constituent, period of 23.93 h) / \(O_1\) (Principal lunar diurnal constituent, period of 25.82 h): Major diurnal tidal constituents associated with the Earth’s daily rotation and orbital inclination.
- \(N_2\) (Larger lunar elliptic semidiurnal constituent, period of 12.66 h): A semidiurnal constituent that accounts for the eccentricity of the lunar orbit.
In coastal aquifers and regions adjacent to the sea, such as southern Beppu facing Beppu Bay, stress variations induced by tidal sea-level fluctuations propagate through the geologic formations. Consequently, even in unconfined aquifers, micro-scale water table fluctuations (on the order of centimeters) are induced in phase with these marine tides, exhibiting identical periodicities (particularly the semidiurnal \(M_2\) component).
In the subsequent part of this series, we apply the Fast Fourier Transform (FFT) to the Beppu field dataset (BP.csv) to quantitatively separate and isolate the barometric and tidal periodic components in the frequency domain.
Summary
- The subsurface structure of southern Beppu hosts multiple hydrothermal flow paths categorized by water chemistry, with faults acting as low-permeability boundaries.
- Long-term decline in the unconfined groundwater level is driven by reduced recharge due to urbanization and vertical leakage induced by reservoir pressure drawdown.
- Shallow water level monitoring serves as a reliable proxy for deep reservoir pressures while offering sensor stability.
- Groundwater time-series contain both long-term seasonal trends and micro-scale (1–5 cm) short-period fluctuations.
- Decomposing these short-period fluctuations into barometric and tidal components requires frequency-domain methods such as FFT and cross-correlation analysis.
By applying the Fast Fourier Transform (FFT) to the 8,136 hourly data points, we identify periodic components that are otherwise invisible in the time-series domain. We will introduce the tidal constituents (\(O_1, K_1, M_2, S_2\)) and examine their manifestations in this unconfined aquifer system.
References
- Yusa Y. (1979) Long-term variations of chemical components in the southern Beppu hot spring area, Reports of the Beppu Hot Spring Research Society, Oita Prefecture, 30, 10-18. (in Japanese)
- Japan Meteorological Agency (JMA) Historical Weather Data Search. (https://www.data.jma.go.jp/obd/stats/etrn/index.php)
- Merritt M. L. (2004) Estimating hydraulic properties of the Floridan Aquifer System by analysis of earth-tide, ocean-tide, and barometric effects, Collier and Hendry counties, Florida. U.S. Geological Survey Water-Resources Investigations Report 2003-4267.
- Yang H.*, Shibata T. (2020) Aquifer classification and pneumatic diffusivity estimation using periodic groundwater level changes induced by barometric pressure. Hydrological Research Letters 14(3): 111–116.