はじめに：別府の地下はどうなっているのか
別府市の地面の下には、世界有数の地熱系が広がっている。地表では至るところから湯けむりが立ち上り、湧出量も泉源数も日本最大級だ。しかし、そこで暮らす人々の足元にある冷たい地下水はいったいどこから来て、どこへ行くのか——そしてなぜ年々水位が下がっているのか。
この連載では、京都大学地球熱学研究施設（別府市）で行った不圧地下水位の観測データを出発点に、信号処理・統計解析・Pythonコードを使いながら、地下水の動きを読み解いていく。
別府の地熱系：温泉水の流れ
別府地熱系の地下構造と温泉水・地下水の流動モデルを以下に示す（図 1）。
別府南部の地下は、大きく分けて3種類の温泉水が異なる深度で流れている。
| タイプ | 水質 | 深度（海面下） | 特徴 |
|---|---|---|---|
| A型 | Ca・Mg-HCO₃型 | −100〜−200 m | 浅部涵養型。雨水起源が強い。 |
| B型 | Na-Cl型 | 約 −250 m | 堀田断層より深部。熱水成分が増加。 |
| C型 | Na-HCO₃型 | 約 −300 m | 最深部。熱水貯留層と関連。 |
さらに地表近くでは、D型（蒸気加熱硫酸塩泉）が存在する。火山ガスが地下水と反応して酸性になったもので、別府の「地獄めぐり」の一部がこれにあたる。
地質の骨格を作っているのは、鶴見岳・伽藍岳の火山噴出物（凝灰岩・凝灰角礫岩）と、その上を覆う堆積層（砂・礫）だ。堀田断層・亀川断層が地下水の流れを区切る「壁」として機能しており、断層の南北で水質や温度が大きく異なる。
不圧地下水位はなぜ下がるのか
この地域で長年観測を続けてきた由佐（1979, 2001）らの研究によると、不圧地下水位は年々低下している。その原因として主に2つが挙げられている。
この2つの原因を区別して評価するためには、地下水位の周期的な変動パターンを詳しく解析する必要がある。そこで登場するのが、今回の研究で用いたFFT（高速フーリエ変換）と相互相関解析だ。
不圧地下水位を観測する意味
「なぜ温泉水頭を直接測らず、不圧地下水位を見るのか？」という疑問が浮かぶかもしれない。理由は明快だ。
- 間接的な温泉水頭の代理指標——不圧層と熱水貯留層は垂直方向に水理的に繋がっており、深部の圧力変化が浅部の水位変動として現れる。
- 地層の透水性を評価できる——水位の変動パターンから、帯水層の水理拡散係数（≈ 透水係数 × 厚さ ÷ 貯留係数）を推定できる。
- 観測が容易で連続記録が可能——深部の熱水は高温高圧で計器が傷みやすいが、浅部の冷地下水位は CTD-Diver（圧力式自動記録計）で安全に長期計測できる。
観測設定：どこで、何を、どのように測ったか
今回の解析に用いたデータは、京大地球熱学研究施設（別府市野口元町）の敷地内に設置した2本の観測井で収集した。
| 観測井 1 | 観測井 2 | |
|---|---|---|
| 地表標高 | 73 m | 80 m |
| 観測井間距離 | 110 m | |
| 平均水位（海抜） | 41.06 m | 41.25 m |
| 観測期間 | 2017年8月19日〜2018年7月23日 | |
| サンプリング間隔 | 1時間（計 8,136 点） | |
| 計器 | CTD-Diver（分解能 2 mm） | |
| 観測項目 | 地下水位・気圧・降水量 | |
2つの観測井は約110m離れているが、平均水位はほぼ同じ（海抜41m前後）で、不圧帯水層の水面が水平に近いことを示している。しかし細かく見ると、気圧変動に対する応答が2井で異なる——この謎が#4で明らかになる。
Pythonでデータを読み込んで最初の一歩
では早速、Pythonを使ったデータ可視化の最初の一歩を踏み出してみよう。
この記事の内容を誰でも手元で再現・学習しやすくするため、ここではPython（numpy）を使って人工的に合成した仮想データを用いて解説を行う。実際の観測データ（CSVなど）をお持ちの場合は、後半で紹介する pd.read_csv() を使ってご自身のデータに差し替えるだけで、全く同じように解析できる。
# 【仮想データ】別府の地下水位・気圧変動を模したダミーデータの生成と概観プロット
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
# フォント設定
plt.rcParams.update({
"font.family": "sans-serif",
"font.sans-serif": ["MS Gothic", "Noto Sans CJK JP", "DejaVu Sans"],
"axes.unicode_minus": False,
"figure.dpi": 150,
})
# ---- 仮想データ生成（実データがある場合はここをpd.read_csv()に差し替え） ----
np.random.seed(42)
n = 8136 # 1時間×8136 = 約1年
t = pd.date_range("2017-08-19", periods=n, freq="h")
# 季節変動（正弦波）＋気圧応答（高周波）＋ノイズ
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は応答が少し大きい
baro = 1010 + 8 * np.sin(2 * np.pi * np.arange(n) / (24 * 14)) \
+ 3 * np.random.randn(n)
# Figure
fig, axes = plt.subplots(2, 1, figsize=(13, 7), sharex=True)
fig.patch.set_facecolor("#F8F9FA")
# 上：地下水位
ax1 = axes[0]
ax1.set_facecolor("#F0F4F8")
ax1.plot(t, gwl_ow1, lw=0.8, color="#2563EB", label="観測井 1", alpha=0.9)
ax1.plot(t, gwl_ow2, lw=0.8, color="#DC2626", label="観測井 2", alpha=0.9)
ax1.set_ylabel("地下水位 (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("別府南部 不圧地下水位（2017–2018年）", fontsize=12)
# 下：気圧
ax2 = axes[1]
ax2.set_facecolor("#F0F4F8")
ax2.plot(t, baro, lw=0.8, color="#374151", label="気圧", alpha=0.85)
ax2.set_ylabel("気圧 (hPa)", fontsize=10)
ax2.set_xlabel("日付", 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.png", bbox_inches="tight")
plt.show()仮想データで見えてきたこと
上記のコードを実行して得られた生データの概観プロット（図 2）を眺めるだけでも、いくつかのことが読み取れる。長周期の季節変動に加えて、短周期の微小な変動が重畳していることが確認できる。
実測データのプロットと解釈
本節では、数値シミュレーションによる仮想データから離れ、実際の観測データを用いた解析を行う。
本解析で使用するデータは、京都大学地球熱学研究施設の敷地内に設置された観測井において、約1年間にわたり1時間間隔で連続計測された地下水位・気圧・降水量の実測値である。解析に用いたCSVデータは以下よりダウンロード可能である。
- 実観測データ: BP.csv
前述の仮想データ生成用Pythonスクリプトを実測データ用に変更するには、スクリプト内の仮想データ生成ブロック（L207〜L224）を、以下のように pandas.read_csv を用いた実データ読み込み処理へと置き換える。
# ---- 実データの読み込み ----
df = pd.read_csv("BP.csv", parse_dates=["datetime"], index_col="datetime")
gwl_ow1 = df["OW1"]
gwl_ow2 = df["OW2"]
baro = df["press"]※CSVファイル内の列名（OW1、OW2、press）に対応させてデータを代入する。
この変更を適用してスクリプトを実行し、実データを可視化した結果を以下に示す（図 3）。
この実測データのグラフからは、地下で生じている物理的なダイナミクスが明瞭に確認できる。
- ベースライン水位と局所的な勾配 両観測井の平均水位は海抜およそ 39.5〜40.0 m を推移しているが、年間を通じてOW2の方がわずかに水位が高い。これは局所的な地下水の流動方向（水力勾配）を示唆している。
- 台風通過などに伴う「気圧応答」 下段の気圧グラフを見ると、例えば9月中旬などに気圧が 985 hPa 付近まで急激に落ち込む谷がある（台風の通過）。この時、上段の地下水位グラフを見ると、気圧の低下と寸分違わず水位が急上昇していることがわかる。帯水層にかかる大気圧の荷重が減ったことで井戸内の水位が押し上げられる現象であり、不圧地下水であっても顕著な気圧応答を示すことが確認できる。
- 長周期の季節変動 人工データで模したのと同じく、数ヶ月スケールでの緩やかなうねり（降水や蒸発散のサイクルによるもの）が存在する。
- ノイズのように見える微小な振動 線の太さの中に隠れるように、1〜5 cmほどの微小な振動が絶えず続いている。これが地球潮汐や日々の微小な気圧変動による「短周期変動」である。
1〜5 cm という振幅は極めて微小であるが、この微小な変動には帯水層の水理特性や気圧応答、そして地球潮汐や海洋潮汐（海水位変動）による影響が包摂されている。
地下水位の短周期変動特性を理解するうえでの典型的な比較対象として、海洋潮汐（海水位）とその周波数成分の関係を以下に示す（図 4）。
図 4 の左図（青線）は、海洋潮汐の1ヶ月間にわたる実測時系列データである。約半日（12時間）および1日（24時間）周期の変動が卓越しているが、その振幅は月と太陽の相対位置（潮汐の干満サイクル）を反映し、約15日周期のうなりを伴いながら時間的に変化している。
この複雑な時間変動をFFT（高速フーリエ変換）によって周波数領域に変換した結果が右下の赤線である。これにより、特定の「主要分潮（Tidal Constituents）」に対応する鋭いスペクトルピークが明瞭に分離される。
- \(M_2\) (主太陰半日周潮, 周期約12.42時間): 月の起潮力に起因する、海洋潮汐において最もエネルギーが卓越する代表的な半日周潮。
- \(S_2\) (主太陽半日周潮, 周期12.00時間): 太陽の起潮力に起因する半日周潮。
- \(K_1\) (太陰太陽日周潮, 周期約23.93時間) / \(O_1\) (主太陰日周潮, 周期約25.82時間): 地球の日周回転と軌道傾斜角に関連する主要な日周潮。
- \(N_2\) (楕円太陰半日周潮, 周期約12.66時間): 月公転軌道の離心率に起因する半日周潮。
海に近い沿岸帯水層や、別府湾に面した別府南部地域においては、潮汐に伴う海水位変動による応力（荷重）が地盤を介して帯水層に伝播する。その結果、不圧地下水位であっても海洋潮汐と完全に同調した周期（特に半日周期の \(M_2\) 分潮など）で、数センチメートル規模の微小な水位昇降が誘起される。
次報では、別府の実測データ（BP.csv）に対してFFTを適用し、時系列データ内に混在する気圧由来および潮汐由来の周期成分の周波数応答を定量的に分離・抽出する手法について述べる。
まとめ
- 別府の地下には深度・水質の異なる複数の温泉水が流れており、断層が流れを支配している。
- 不圧地下水位は都市化と熱水貯留層の圧力低下により年々低下している。
- 不圧水位を観測することで、温泉水頭の変動を間接的に把握できる。
- 観測データには季節スケールの大きな変動と、1〜5 cm の短周期変動が重なっている。
- 短周期変動の原因（気圧？地球潮汐？）を分離するには、FFTと相互相関解析が必要。
8,136点の時系列データにFFTをかけると、肉眼では見えない周期成分が一目瞭然になる。O₁・K₁・M₂・S₂という潮汐分潮の名前と、それが不圧帯水層でどう現れるか（あるいは現れないか）を解説する。
参考文献
- 由佐悠紀 (1979) 別府温泉南部域の化学成分長期変化について, 大分県温泉調査研究会報告, 30, 10-18.
- 気象庁ホームページ, 過去の気象データ検索. (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.