行政院國家科學委員會補助專題研究計畫成果報告:地下水觀測站網水位急遽變化與地震發生關連性之先期研究：地震波對地下水觀測井孔隙水壓力學變化

行政院國家科學委員會補助專題研究計畫成果報告

※※※※※※※※※※※※※※※※※※※※※※※※※※

※

地震波對地下水觀測井孔隙水壓力學變化影響之研究 ※

※ The Investigation Resear ch of the Seismic Por e Pr essur e

※

※

Mechanics Var iation Influence on the Monitor ing Wells

※

※※※※※※※※※※※※※※※※※※※※※※※※※※

計畫類別：□個別型計畫

□整合型計畫

計畫編號：NSC 89 – 2116 – M – 053 – EAF

執行期間：89 年 08 月 01 日至 90 年 07 月 31 日

計畫主持人：譚義績 教授

共同主持人：郭嘉真 副研究員

計畫參與人員：

林允斌 博士生

余化龍 碩士生

馬國宸 碩士生

本成果報告包括以下應繳交之附件：

□赴國外出差或研習心得報告一份

□赴大陸地區出差或研習心得報告一份

□出席國際學術會議心得報告及發表之論文各一份

□國際合作研究計畫國外研究報告書一份

執行單位：國立台灣大學農業工程學研究所

中

華

民

國

90

年

7

月

31

日

地下水觀測站網水位急遽變化與

地震發生關連性之先期研究

地震波對地下水觀測井孔隙水壓力學變化影響之研究

The Investigation Resear ch of the Seismic Por e Pr essur e

Mechanics Var iation Influence on the Monitor ing Wells

計畫編號：NSC 89-2116-M-053-EAF

執行期限：89 年 8 月 1 日至 90 年 7 月 31 日

主持人：譚義績教授 國立台灣大學農業工程學研究所

共同主持人：郭嘉真副研究員 國家高速電腦中心

計畫參與人員：林允斌博士生、余化龍碩士生、馬國宸碩士生

國立台灣大學農業工程學研究所

一、中文摘要 本研究利用簡化的黏彈性模式，描述 含水層中粗顆粒礫石沈積物與細顆粒砂石 沈積物交界面處至拘限含水層尖滅端，於 地震力作用下液相之動力行為。假設拘限 含水層中砂石沈積物所在研究區域於震動 時為液相可壓密且孔隙介質不排水狀況 下，將液相及固相耦合的波動方程式分 離，僅考慮液相控制方程式，推導出解析 解，同時套配 921 集集大地震中於濁水溪 沖積扇中所觀測到之水井水位變化資料， 獲致良好的結果。本研究亦參照地質井取 樣調查所描繪的水文地質剖面圖，描述 921 集集大地震中發生土壤液化地區，與含水 層中粗顆粒礫石沈積物與細顆粒砂石沈積 物交界面處位置關係，據觀察紀錄顯示， 地震所引起之水井水位變化最大處皆位於 含水層中礫石與砂石交界面處，而導致該 區域靠近扇頂部分於地震後數日陸續發生 土壤液化之現象，顯示於大地震中評估土 壤液化潛能，除現有工程評估方法外，還 需考慮工程所在水文地質環境之影響。 關鍵詞：集集大地震、水井水位變化、黏 彈性力學、孔隙介質、土壤液化 Abstr act

The simplified visco-elastic model is adopted to describe the earthquake induced dynamic behavior of the liquid phase from

the interface between the gravels and the sands to the pin-off in the confined aquifer. We assume the porous media is undrained and the liquid phase is compressible in the sanded confined aquifer during the vibration, and then decouple the momentum equation of the liquid phase to find the analytical solution. The presented model is derived to match the captured data of the high dense network of hydrologic monitoring wells after the Chi-chi earthquake and does fairly well. Consulting the hydrogeological profiles identified from the drilled cores, the highest well water level changes are located at the interface between the gravels and the sands in the confined aquifer and induce the liquefaction at the upstream area in the next few days. The field observation emphasizes the importance of the hydrogeological influences in evaluating the liquefaction potential.

Keywor ds: The Chi-Chi earthquake, Well

water level changes,

Visco-elasticity, Porous media, Liquefaction 二、緣由與目的 緣由 車籠埔斷層發生嚴重的錯動所引發之 九二一集集大地震，斷層的錯動造成巨大 的直接傷害，而伴隨地震發生後之土讓液 化則擴大受災之面積。藉由現代化觀測儀 器，此次大地震所產生大量之地下水及地

質資料被記錄下來，若能善加探討分析， 則對未來台灣地區的地震防治工作將有所 助益。現今地震防治相關研究的最終目 的，皆冀望在強震來臨前，能爭取到半小 時至一小時的預警時間，以降低人員傷亡 及財產損失。而分析此次集集地震所產生 的眾多資料，其中最具預報潛力的資料即 為「台灣地區地下水觀測站網」於車籠埔 斷層附近地區的地下水觀測井水位資料， 該資料之重要性在於其為現今世界上存在 有關沖積扇地下水受地震影響之唯一完整 及長時距高密度之記錄。 台灣位於歐亞大陸板塊與太平洋的菲 律賓海板塊交界之處，屬世界上有感地震 最頻發的地區之一，地震發生時，除因地 震波及地表永久性位移所造成對建築物結 構的直接傷害外，伴隨地震而引發的大規 模地區之土壤液化所造成的災害，其範圍 與規模並不亞於地震所造成之直接傷害。 目的 強烈地震所造成之土壤孔隙水壓增加 是造成土壤液化的主要原因，若水力傳導 係數低，則土壤水分為吸收地震波極佳的 材質，其阻尼效應十分良好，因震動所產 生的土壤高孔隙水壓多發生於砏土及砂質 土壤交界區域，地下水位因此急遽變化， 而藉由地下水觀測站網可以很輕易的觀測 孔隙水壓隨地下水位之變化而改變，同時 根據觀測站網所紀錄的資料顯示，而地震 後土壤孔隙水壓之消散，造成地下水觀測 井觀測水位隨時間之變化。 本研究擬先針對地震波與井群及含水 層間關連性進行基礎研究，其後藉由對地 震後地下水觀測井觀測水位資料隨時間變 化之分析，可以充分掌握地下水位變化與 地震發生之關連性。水文地質研判及其參 數之取得，以往皆利用地下水觀測井所觀 測之資料進行內差，往往無法掌握整體趨 勢之變化，透過對地震前後地下水觀測井 觀測水位資料空間的分析，可以彌補以往 判識上之誤差。研究目的乃期望能將研究 成果應用於震災潛勢分析與預警系統。 三、研究內容

前言

經整理[Roeloffs（1996）]為數眾多紀 錄及研究有關地震引起的水文變化現象之 學 術 論 文 ， 皆 利 用 孔 隙 彈 性 力 學 （poroelastic）之原理來描述地震所引起局 部區域之孔隙水壓變化，並試圖以地下水 位 變 化 做 為 地 殼 變 形 之 指 標 [Ohno and Wakita（1997）]，唯因拘限於觀測井之數 目 無 法 建 立 整 體 的 關 連 性 [Roeloffs （1998）]，或因模擬尺度的關係導致相關 性欠佳 [Grecksch et al.（1999）]，另外 Cooper[Cooper et al.（1965）]以水井為控制 體積，推導出水井水位變化與拘限含水層 孔隙水壓變化之比，於本研究中皆假設為 1。本研究嘗試結合高密度的觀測站網紀錄 資料，配合簡化之黏彈性力學模式推導解 析拘限含水層中之液相孔隙水壓，描述沖 積扇地區因地震引起之地下水整體變化分 布情形。

前人研究

連體力學中單相物質之本質特性參數 Ψ之控制方程式可表示為

( )

(

)

f G t ρΨ +∇⋅ ρ Ψ −∇⋅ −ρ =ρ ∂ ∂ i v （1） 而在α 、β 不同相（氣、液、固三相） 的介面上存在如下的不連續之介面方程式

(

)

(

)

(

)

(

)

0 αβ α βα β ρψ ρψ − + ⋅ + − + ⋅ = w v i n w v i n （2）

而 Hassanizadeh [Hassanizadeh and Gray（1979a、1979b、1980）]藉由 Coleman [Coleman and Noll（1963）]之方法利用熱 力學第二定理熵不等式原理，將二相流之 動量控制方程式線性化，其基本假設為 1. 於不同相之介面上無相變之產生；2.液相巨 觀無摩擦性；3.溫度梯度僅與熱傳導有關； 4.忽略加速度中的傳導項即u_lu_k,_l項；5.毛細 效應可忽略。 則液相及固相之線性動量控制方程式 分別為

2 , 2 0 f f k f f f f f k f k f d km m u g p t R ε ρ ε ρ ε υ ∂ _{− +} ∂ − = （3a）

(

)

(

)

(

)

2 2 , , 1 1 1 0 s s k s s f f k f s f d f k kll km m u g t p R ε ρ ε ρ ε σ υ ∂ − − − + ∂ − − + = （3b） 唯參照隨後相關研究[Biot（1956a、

1956b）、Berrymman et al.（1988）、Jeng and

Lee（2001）]指出，控制方程式式（3）於 二 相 流 震 波 傳 遞 應 用 時 ， 即 模 擬 動 態 （dynamic）反應時應存在額外的控制項， 該項之存在應與液相與固相之相對加速度 有關，則式（3）可表示為 2 , 2 0 f f k f f f f f k f k f d s d km m km m u g p t R A a ε ρ ε ρ ε υ ∂ _{− + −} ∂ − = （4a）

(

)

(

)

(

)

2 2 , , 1 1 1 0 s s k s s f f k f s f d s d f k kll km m km m u g t p R A a ε ρ ε ρ ε σ υ ∂ − − − + ∂ − − + + = （4b） 相 對 於 Biot 的推導[Biot（1956a、 1956b）]，其利用流體力學中 Poiseuille 流 體的理論並結合波動力學之基礎，推導出 當飽和孔隙介質於波動中可傳遞兩種壓力 波，其波速為一快一慢，同時由於液相不 傳遞剪力波故僅有一種剪力波能通過孔隙 介質，其耦合方程式可表示為

(

)

(

)

, , 2 11 12 2 1 s f kll f k s f d k k m p u u b t σ ε ρ ρ υ − − = ∂ _{+ −} ∂ （5a）

(

)

d m f k s k f k f u u b t p ρ ρ υ ε + + ∂ ∂ = − 2 12 22 2 , （5b） 比較式（4）及式（5）之結果可知： 式（4）在忽略重力項之影響後，ρ₁₁+ρ₁₂為 單 位 孔 隙 介 質 體 積 中 固 相 之 密 度 ， 即 s sρ ε ；ρ₁₂ +ρ₂₂為單位孔隙介質體積中液 相 之 密 度 ， 即 f fρ ε ； 則 s km A = 12 ρ ， s km s s A − =ρ ε ρ11 ， s km f f A − =ρ ε ρ22 ， f km f k R b=µε 2 = 。 若假設式（5）為等向性， s km A 的存在 除 Biot 的推導[Biot（1956a、1956b）]可證 明 外 ， 根 據 Represent Theorem[Wang

（ 1970a 、 1970b ） 、 Drew and Passman （1999）]，當一方程式所描述的物理量具 等向性時，則其變數應不受座標轉換之影 響（Objective），同時由於速度差及加速 度差皆不受座標轉換影響之向量，因故為 完整描述控制方程式，應於方程式中加入 d m s kma A 項。除此之外ρ 可視為當固相顆粒₁₂ 於液相中震盪時所引起液相動能增加之質 量，相關理論可以於流體力學非旋性勢能 流中獲得應証[Lamb（1945）]。

Berryman [Berryman et al.（1988）] 利用變分法及前述之方程式，推求不混合 多相系統（含固、液及氣相，即非飽和土 壤）之波動線性方程式，其型態近於式 （5a）、（5b）取 Divergence 及 Curl，僅 係數項較二相流控制方程式複雜，同時利 用 Helmholtz 轉換將壓力及剪力之波動方 程式分離成為 2 2 2 2 2 t w k_s ∂ ∂ = ∇ χ χ （ 一 組 剪 力 波 方 程 式 ） （6a） 2 2 2 2 2 t A w k A ∂ ∂ = ∇ ± ± ± （ 兩 組 壓 力 波 方 程 式 ） （6b）

控制方程式

假設等向性，將式（5）取 Divergence 可以表示為

(

)

(

)

(

)

2 2 2 11 12 2 1 s f f f f d p p e b e b t t σ ε ε ρ ρ ς ς υ ∇ − − ∇ + ∇ ⋅∇ = ∂ _{+ −} ∂ _{− − ∇ ⋅} ∂ ∂ （7a）

(

)

(

)

2 2 12 22 2 f f f f d p p e b e b t t ε ε ρ ρ ς ς υ − ∇ − ∇ ⋅∇ = ∂ _{+ +} ∂ _{− + ∇ ⋅} ∂ ∂ （7b） 其中 s k u e=∇⋅ ， f k u ⋅ ∇ = ς 。 式（7a）與式（7b）之和為

(

)

(

)

(

) (

)

(

)

2 2 11 12 22 12 2 s f f sp f p e t σ ε ε ρ ρ ρ ρ ς ∇ − − = ∂ _{+ + +} ∂ （8） 依 據 混 合 理 論 （ mixture theorem ） [Green（1965）]，考慮孔隙水壓及土壤有 效應力，可將式（8）分離成兩組耦合之壓 力波

(

)

(

(

)

)

( )

(

( )

)

2 2 11 12 2 2 2 2 2 s f s s s s s s p e t e u t t σ ε ρ ρ ρ ε ρ ε ∂ ∇ − = + = ∂ ∂ ₌ ∂ _{∇ ⋅} ∂ ∂ （9a）

(

)

(

(

)

)

( )

(

(

)

)

2 2 22 12 2 2 2 2 2 f f f f f f f p t u t t ε ρ ρ ς ρ ε ς ρ ε ∂ ∇ − = + = ∂ ∂ _≅ ∂ _{∇ ⋅} ∂ ∂ （9b） 於 此 處 假 設 f f p ε ∇ ⋅∇ 、

(

)

2 12 12 2 e t ρ ρ ς ∂ − + ∂ 和 d b υ ∇⋅ 於 同 震 時 期 相 對於其他項可忽略並視為零。也就是說， 將 b 和ε 視為常數（與 Biot 推導假設相同_f [Biot (1956a, 1956b)]），同時式（9）因忽 略相對加速度項故可視為 quasi-static 控制 方程式。其中，於此時先忽略 b

(

e

)

t ς ∂ ₋ ∂ 項，並於下一節之本構方程式中利用阻尼 係數來描述孔隙水壓因相對速度而消散之 現象，亦即利用黏彈性模式取代b

(

e

)

t ς ∂ ₋ ∂ 項之效應。 根據研究指出於假設所在研究區域 於震動時為不排水狀況下，可僅考慮液相 行為[Prevost（1982）]，假設拘限含水層中 砂石沈積物所在研究區域於震動時為液相 可壓密且孔隙介質不排水狀況下（研究區 域含水層分層概念圖詳圖一），將液相及 固相耦合的波動方程式分離，僅考慮液相 波動方程式，即式（9b）。

本構方程式

有關飽和孔隙彈性土壤的動力行 為可用黏彈性力學模型來簡化[Eringen （1980）、Bardet（1992）]，可得出相 近的模擬結果，有關數值模式的模擬亦 有類似的處理[Zienkiewicz（1984）， scheme B]，於液相不可壓縮的假設下將 penalty 項平移至阻尼項。本研究依前述 僅考慮液相波動方程式（9b）同時假設液 相於孔隙中為黏彈性流中之 Maxwell-Fluid （其概念模型詳見圖二）。 則壓應力與應變之關係可表示如式 （10）

(

)

(

)

(

( )

f

)

f f f f f u t p t C p η ε η ε ε ∇⋅ ∂ ∂ = − ∂ ∂ + − （10） 利用阻尼係數η 來描述孔隙水壓因固 相及液相相對速度導致摩擦消散之情況。 這樣的作法其優點在於不需要完整描述固 相顆粒之運動情形，而僅需利用參數化的 控制方程式去描述所觀測之現地資料，然 後利用現有已知之物理係數（例如：水利 傳導係數）去解釋現地資料套配所求得之 參數之物理意義。 同時考慮一維水平流況，設此單位孔 隙 介 質 體 積 中 液 相 之 水 平 位 移 為 (lx wt) i x f fu e e − − = α ε ，則單位孔隙介質體積中 液相之應變速率可表示為

( )

_iw

_{( )}

_u

₍

_li

₎

t u _f f f f _{=− − +} ∂ ⋅ ∂∇ α ε ε （11） 假設單位孔隙介質體積中液相之壓應 力與單位孔隙介質體積中液相之應變速率 之關為

( )

t u p f f f f _∂ ⋅ ∂∇ Λ = −ε ε ，帶入式（10） 得 C i wη η − = Λ 1 （12） 將式（11）、（12）代入式（9b）得

(

)

(

)

(

)

(

)

2 2 2 1 f f f f f w i li u w i C li u t η _{α ε} η ρ α ε − _{− + ∇ =} − ∂ − + ∂ （13） 比較式（6b）

2 2 2 2 2 t A w k A ∂ ∂ = ∇ ± ± ± （兩組壓力波方程式） 若 假 設 k 為 液 相 波 之 波 向 量 及₋ f fu A₋=ε ，則

(

Qi

)

v w Q i Q i v w C i w C i w C w k f f f + = − = − = − 1 1 1 2 2 2 2 2 2 η η ρ （14） 其中 Q C wη ₌ 1 ，根據式（14）及式（10） 可將式（9b）之控制方程式表示為 0 1 2 2 2 2 2 = ∂ ∂ − ∂ ∂ − ∂ ∂ _{− − −} t A v t A x A f f η ρ （15）

拉普拉司轉換求解初始條件問題

假設濁水溪沖積扇地下水分區拘限含 水層砂石沈積物所在研究區域，於震動時 之概念模型為一端受力之黏彈性體（概念 模型詳見圖二）。 式（15）的初始條件為

( )

,0 =

( )

,0 =0 • − − x A x A （16） 即震動前之單位孔隙介質體積中液相 之位移及速度皆為零，將式（15）作拉普 拉 司 轉 換 將 t 轉 換 為 s 、 A₋

( )

x,t 轉 換 為

( )

x s A₋ , ，並將初始條件式（16）代入式（15） 則可表示為

( )

,

( )

, 1 2

( )

, 0 2 = − − ″ − − − s A x s v s x A s s x A f f η ρ （17a） 設

( )

x e s B A₋ = λ （17b） 式 （ 17b ） 代 入 式 （ 17a ） 得     ₊ ± =       + ± = s C s v s C s f f η η η ρ λ 1 1 為 符 合 物 理 現 象 ， 及 當 x→∞ 時 0 → − A ，故λ 取負值，同時由本構方程式 知，於x=0及t=0突 然施加應力 Pδ

( )

t 時 ， 其 邊 界 端 之 應 變 為

( )

( )

η η δ _t C t P t A + = ′ − 0, ，作拉普拉司轉換將 t 轉換為 s ，得

( )

B

( )

s C P s A₋′ 0, = 1 =λ ，並代 入式（17b）可將通解表示為

( )

−  +  −     ₊ − = s C s v x f f e s C s C Pv s x A η η , （18）

反拉普拉司轉換

經由查表知[Farrell, O. J. and Ross, B. （1971）]將 f

( )

t =I

(

a t2 −b2

)

H

( )

t−b 0 作拉 普 拉 司 轉 換 將 t 轉 換 為 s 可 得

( )

₍

₎₍

₎

b (s a)(s a) e a s a s s f − − + + − = 1 ，利用平移 定 理 可 求 將 f

( )

t =I0

(

a t2−b2

)

e−atH

( )

t−b 作 拉 普 拉 司 轉 換 將 t 轉 換 為 s 可 得

( )

₍

₎

b s(s a) e a s s s f 2 2 1 _{− +} + = ，同時將 η 2 C a= 及 f v x b= 代入可求解式（18）之反拉普拉 司轉換得

( )

2 2 2 0 , 2 C t f f f A x t Pv _{C x x} I t e H t C v v η η − − =  _{ }  _{ }   − _ −_ _{ } _ − _         （19） 依據式（12）

( )

(

∇⋅ −

)

∂ ∂ − = − A t C i w pf f η η ε 1 可求得則單位 孔隙介質體積液相壓力落後單位孔隙介質 體積液相應變速率之相位差為     = − Q 1 tan 1 θ ，同時將單位孔隙介質中液相 之壓力表示為

2 1 ₂ 1 ₂ 2 2 1 1 2 3 2 _{2 2} 1 2 1 1 2 2 2 0 2 2 1 2 ₂ 1 1 2 2 2 ₄ 2 2 8 f f f C C t t f _f C t f f p Pv C Q C C txI e xI e v _v C _C C C tx I I e v C x H t Z v η η η ε η η η η _η η η η − − − = +  _{   }  ∆ ∆  _{   }       ₊   _{ }  ∆  _{ } ∆                _{∆ + ∆}        _{    }  ₋   _{ }  ∆             − +       其中 2 2         − = ∆ f v x t （20） 式（20）之待定常數 Z 之求解應代入 含水層尖滅端應變速率為 0 之邊界條件（此 假設之合理化可由現地土讓液化現象加以 說明），即

(

∇⋅

( )

,

)

=0 ∂ ∂ − Lt A t ，得單位震源 壓力於拘限含水層所砂石沈積物所在區產 生之液相壓力分布為Unitf

( )

xt , [無因次式]

( )

2 2 2 1 1 2 3 2 _{2 2} 1 2 2 0 2 2 1 2 ₂ 2 1 2 , 1 1 2 2 2 ₄ 2 2 8 2 2 f f f C C t t f _f C t f C t L f Unit x t v t C Q C C txI e xI e v _v C _C C C tx I I e v C C tLI e v C η η η η η δ ε η η η _η η η η η η − − − − = +  _{   }  ∆ ∆  _{   }       ₊   _{ }  ∆  _{ } ∆       _{ }       _{  ∆ +  ∆}   _{    }  −   _{ }  ∆            _∆      ∆ − 2 1 2 3 ₁ 2 2 ₂ 2 0 2 2 1 2 ₂ 2 4 2 2 8 C t L L _{f L} C t L L f L f C LI e v C C C tL I I e v C x H t v η η η η η η η − −                                   _∆    _{ }         ₊     _{ }   ∆   _{ }           _{    }     _{  ∆ +  ∆ }     _{    }  −         _∆       _{   }        −       當 f v x t> 時 （21） 其中 2 2         − = ∆ f v x t 、 2 2         − = ∆ f L v L t 、 tδ 為震動時間。 相關研究[Hsieh et al.（1987）]指出孔 隙介質的儲水係數對水井水位變化影響不 大，由式（21）可知當孔隙介質孔隙率越 大時，或η 越小時，單位震源壓力於含水 層所產生之液相壓力越小。

交界面上孔隙水壓消散之解析

拘限含水層中礫石與砂石交界處因地 震導致急遽升高孔隙水壓之作用力為 F ， 當地震停止後，依據理論[Biot（1941）]，

f bc f p ε 即介面上單位孔隙介質體積中液相 之水壓，其消散行為符合擴散方程式，可 描述為：

(

)

(

)

_{( ) (}

₎

2 2 f f bc f f bc p t p kb F x t t S x ε ε δ δ ∂ − = ∂ ∂ − ′ + − ∂ （22） 初始條件p_bcf

( )

x,0 =0 令 f bc f p h=−ε ，利用一次拉普拉司及 一 次 傅 利 葉 轉 換 求 解 [ Abramowitz and Stegun（1972）]。 式（22）利用一次拉普拉司轉換將 t 轉 換成 s 、 h 轉換成 h ，並代入初始條件得

( )

st e x F x h S kb h s + − ′ ∂ ∂ = 2 δ 2 （23） 式（23）利用一次傅利葉轉換將 x 轉 換成ξ 、 h 轉換成 h 得 t s e F h S kb h s =− + − ′ π ξ 2 1 2 （24） 可解得                   +             = + = − ′ ′ − S kb Fe S kb s s S kb s S kb s e S F h t s t s π ξ ξ ₂ 2 2 2 （25） 式（25）利用反傅利葉轉換將ξ→x得       = − ′       − S kb Fe s e h t s x S kb s π π 2 2 （26） 式（25）利用反拉普拉司轉換將s→t 得

( )

₍

( )

₎

      ′ − − = ′ −       − S kb F t t e t x p f t t S kb x f bc π ε 2 , 4 2 （27） 當脈衝施於粗顆粒及細顆粒沈積物介 面瞬間，液相位移沿施力方向瞬間移動， 並以波傳方式將能量傳遞至扇尾部分，因 液相與固相顆粒摩擦導致震波能量之損 失，而震幅之平方正比於震波能量，故液 相位移隨震波傳遞距離及時間增長而逐漸 衰減（詳圖三），當震波通過以後液相會 因壓力不平衡及其流動之特性，而趨向回 復初始狀態及其速度沿施力方向之負方向 移動，唯因介面上液體孔隙水壓之消散隨 時 間 之 變 化 為

( )

₍

₎

      ′ − = S kb S F t t t p f f bc π ε 2 1 , 0 ，故細顆粒 沈積物所處之拘限含水層之液壓 f res p 應為

( )

,

( )

0,

(

, τ d

)

τ 0

∫

− = t _{f f} bc f res xt p tUnit xt p （28）

實例應用

背景說明 台灣於 88 年 9 月 21 日凌晨 1 點 47 分 12.6 秒，於南投集集附近發生強烈地 震，震央位於北緯 23.87 度、東經 120.75 度，芮氏地震規模達到 7.3 ML[Ma et al. （1999）、Yu et al.（2001）]，及於 88 年 9 月 26 日上午 7 點 52 分 50.0 秒於北緯 23.86 度、東經 121.00 度，發生芮氏地震規模達 6.9 ML之地震，位於台灣地區濁水區沖積 扇地下水分區之地下水觀測井皆完整記錄 此兩次地震前後之地下水位。 本研究區域位於台灣西部平原地帶， 其面積廣達 1700 平方公里，其東面為一連 串 的 西 向 逆 推 斷 層 （ west-vergent thrust faults），其造成鄰近沖積扇東面的八卦山 脈及斗六丘陵等背斜構造（anticlines）， 該區共有 70 個地下水文觀測站計 188 口水 井，每口水井僅於單一含水層設置濾網， 並以壓力式紀錄儀記錄每小時的水位[Hsu （1998）]。（詳圖四） 經濟部中央地質調查所[經濟部中央 地質調查所（1994、1999）]於研究區域共

進行十二個剖面的地質井鑽探調查，依據 地質調查井採樣繪製水文地質頗面分析顯 示，濁水溪沖積扇地下水分區可分為三個 含水層，各含水層位於扇頂位置為非拘限 含水層，而位於扇央及扇尾位置則屬拘限 含水層，各含水層於拘限含水層中又存在 粗顆粒礫石沈積物與細顆粒砂石沈積物交 界面（詳圖一），礫石沈積物靠近扇央部 分而砂石沈積物則位於扇尾，本研究假設 各含水層間之阻水層無滲流（leakage）產 生，亦即假設震波由扇頂傳遞至扇尾時於 含水層之拘限含水層部分為不排水狀態。 本研究為配合長時距模擬（50 日）， 故採用日平均水位，同時於此時間尺度之 下，可將地震為期數十秒的震動，以相對 應震動產生所產生之應力，簡化為單一脈 衝輸入項進行模擬。 集集大地震中車籠埔斷層地表破裂主 斷層長約 85 公里呈南北走向，約略與各含 水層之等水位線平行，根據日平均水位記 錄顯示，地震引起之水井水位變化最大處 皆位於各含水層中礫石與砂石之交界面處 （詳圖五），假設車籠埔斷層破裂帶對含 水層施予由東向西之衝擊力。本研究模式 應用於現地資料之套配，以一維模式模擬 東西向於礫石及砂石之交界處孔隙水壓消 散情形（詳圖二），利用位於鄰近交界處 東方之觀測井水位紀錄套配式（27），利 用該東方水井之L 、 k 、 S 反推求礫石與_bc 砂石交界面上孔隙水壓消散之狀況，再利 用前述推求出之結果視為邊界條件 f bc p 並 代入式（28）中，套配位於交界處西方之 觀測井水位紀錄，利用該觀測井之L 推求_bc η 值。 本研究共套配四組東西向之水井，其 編號如 3W1、3W2，其意義為，第一個數 字代表該觀測井濾網所在之含水層位置，3 代表含水層 3；第二個英文字表示該觀測井 位於交界面處之東西方，E 代表東方；第 三個數字表示組別，每一組有東西方各一 口水井（水井所在位置及特性參數詳表 一）。 套配結果 本研究於濁水溪沖積扇北側，含水層 三中礫石與砂石之交界面所在位置東西各 取一水井水位紀錄，分別為 3E1、3W1 及 3E2、 3W2 進 行 模 擬 （ 水 井 位 置 詳圖五 （c）），將 3E1 及 3E2 水井水位紀錄套配 式（27），3W1 及 3W2 水井水位紀錄套配 式（28），由於交界面位於拘限含水層中， 故孔隙水壓於交界面之消散符合式（27） 可視為一維之消散變化，同時本模式模擬 東西方水井水位結果皆符合地震後之水井 水位消散之變化趨勢（模擬結果詳圖六、 七），其差異除因本模式以簡化之一維模 式模擬其物理現象外，同時必須考慮於震 動過程中含水層之砂石沈積物所在區域其 阻水層並非完全無滲流，導致各含水層間 互相影響。本研究著重於孔隙水壓之效 散，故 quasi-static 模式足以描述其變化情 形。 本研究另於沖積扇南側，含水層二及 三中礫石與砂石之交界面所在位置東西各 取一水井水位紀錄，分別為 3E3、3W3 及 2E1、2W1 進行模擬，其中 3W3 及 2W1 位於同一位置，結果顯示（模擬結果詳圖 九，水井位置詳圖五（b）、五（c））， 水井 2E1 於地震後之水位消散狀況較本模 式預測快，亦即於交界面東方之孔隙水壓 消散不僅為一維之變化，根據中央地調所 調查結果得知，含水層一地下水流出區位 於水井 X 之位置（詳圖五（a）），9 月 21 日之地震造成水井 X 下游河道流量突升 （詳圖八），亦即地下水大量湧出，其過 程導致含水層一於水井 X 處水井水位受出 流地下水大量湧至而導致水位突升，含水 層二於水井 X 處水井水位受垂直流況影響 增加導致水井所量測之水位突降，而 2E1 受水井 X 流況之影響，為補充含水層二水 井 X 處垂直流向出流量之損失，導致後期

退水較理論值明顯增快。位於鄰近位置之 3E3，卻因含水層三因與含水層二間於水井 X 處存在阻水層，故 3E3 之水位變化趨勢 符合式（27）一維之假設。相同之現象亦 發生於水井 Y 處，含水層二及含水層三於 此處無阻水層相隔，含水層三水井水水位 受垂直流況影響增加導致水井所量測之水 位增加量減低，即下層含水層水位變化受 其上層阻水層位置之影響（詳圖五（b）、 五（c）），由上述兩個例子可以顯示當地 震發生時於非拘限含水層中，垂直之流況 為一重要的現象。另一方面，配合 3E3 及 2E1 之水井水位變化以式（28）套配 3W3 及 2W1 之水井水位變化獲致良好之結果 （模擬結果詳圖九）。 同時配合 921 集集大地震後，於濁水 溪沖積扇部分地區產生土壤液化現象，依 據發生之區域可分為兩類（詳圖五（a））。 發生於扇頂處非拘限含水層的土壤液化現 象（詳圖五（a）中 A、B1、B2、B3）， 明顯受到非拘限含水層中垂直流況之影 響，同時土壤液化之區域處於含水層中礫 石及砂石交界面處，亦即交界面處產生的 局部高水位透過非拘限含水層中垂直流況 之影響，導致所在區域之土壤液化。另一 種型態為發生於沖積扇北側靠近海邊地區 之土壤液化現象（詳圖五（a）中 C、D）， 沖積扇北側靠海部分處於沖積扇扇央的位 置，亦即扇尾位於離海岸線約 25 公理處， 同時由於靠近海岸邊，含水層頂部無黏土 不透水層，以本模式理論，拘限含水層中 水平流況之壓力於此處得到宣洩，導致該 區土壤液化。同時位於沖積善南側靠近海 邊地區，因處於含水層尖滅之扇尾端，雖 同樣無含水層頂部黏土層，並無發生土壤 液化現象，合理假設為含水層尖滅端應變 速率為 0。根據觀測資料顯示，當震波由含 水層中礫石沈積物所在區傳遞到砂石沈積 物所在區域，交界面處產生急遽升高之孔 隙水壓，若含水層其上方無阻水層則易導 致噴砂或土壤液化之現象，因此於評估區 域土壤液化潛能時需考慮該地區水文地質 之影響。 有關簡化模式中阻尼係數η 之物理意 義，根據前述四組八口東西方向之井，套 配式（28）求得       C Ln η 及kL 之關係為（詳_bc 圖十） 943 . 2 6401 . 1 1000 +     − = C Ln kL_bc η （29） 相關係數達 0.9831，式（29）表示當 水井距交界面越遠，受交界面之影響越 小，亦即η 值越小，另當孔隙介質之水力 傳導係數越大時，液相與固相越易產生相 對位移，因故其孔隙水壓消散能力越佳， 亦即η 值越小。同時地震所引起之水井水 位變化最大處皆位於含水層中礫石與砂石 交界面處，此關係可由式（21）得知，當 處於交界面之鄰近兩端，若砂石孔隙率ε_f 較小，及其因水力傳導係數較小導致阻尼 係數η 較大時，會產生較大之孔隙水壓， 其孔隙水壓消散也較慢。

結論

本研究利用阻尼係數來描述當地震應 力移除後，因液相與固相之相對位移導致 液相壓力消散的情況，建立一維模式模擬 此次 921 集集大地震後水井水位消散之現 象，獲致良好的結果。透過觀察地震後研 究區域之水井水位變化分布，利用本模式 於拘限含水層中，若已知地震作用力之方 向而不須知道地震作用力之大小，沿作用 力之方向上選取兩觀測井分別位於粒石及 礫石交界面處兩端，則可利用觀測井於地 震後之水位紀錄，配合水井試驗求得之水 力傳導係數及儲水係數，判別粒石及礫石 交界面處所在之位置。

符號表

± A ：勢能值[L ]2 d m a ：液相與固相之 m 方向加速度差 [L T2]

b ：含水層或水平流厚度[ L ] C ：為液相之統體壓密係數[

( )

2 LT M ] f ：物質之外來直接供給率[Ψ T] G ：物質的淨產生率[Ψ T]

( )

H ：Heaviside 函數

( )

n I ： n 階的 Bessel Function，即

( )

x J

( )

ix I_n = _n i ：物質於邊界之流通量 [

(

M L3

)

×

( ) ( )

L T × Ψ ] k ：液相之水力傳導係數[L T] ± k k_s, ：波向量（wave vector）[ L1 ] bc L ：位於礫石與砂石沈積物交界面兩側 之水井與交界面之水平東西向距離[ L ] l ：波向量[ L1 ] αβ n ：介面上由α 相指向 β 相垂直於 介面的單位向量 f k p, ：液相於 k 方向之壓力梯度 [

(

3

)

( )

2 T 1 L M × ] f km R ：液相與固相因速度差產生之動量交 換係數[

(

M L3

)

×

( )

1T ]，若假設參考座標 平行主軸則其為一純量即等向性 S ：儲水係數（coefficient of storage） α k u ：α 相於 k 方向之位移[ L ] v ：物質於空間上一點之速度[L T] w ：介面移動速度[L T] w ：震動頻率[ T1 ]

希臘字

Ψ ：本質特性參數 α ：震波之衰減係數[ L1 ] α ε ：α 相所佔體積比例（α = f 表示液 體、α =s表示固體） d m υ ：液相與固相之 m 方向速度差[ TL ] µ ：液相之黏性係數[M

( )

LT ] f v ：液相波傳遞之速度[L T]，為 C ρf ρ ：物質密度[ 3 L M ] α ρ ：α 相之密度[ 3 L M ] l s kl , σ ：固相所受表面應力於l 方向之應力 梯度[

(

3

)

( )

2 T 1 L M × ] χ ：勢能值[ 2 L ] η ：阻尼項係數，即壓力與應變速率之 比值[M

( )

LT ] 五、參考文獻 1. 經濟部中央地質調查所，1999，濁水溪沖積扇 水文地質調查研究總報告。 2. 經濟部中央地質調查所，1994，濁水溪沖積扇 水文地質調查研究報告。 3. 經濟部水利處，1999a，台灣地區地下水觀測站 網地下水位年報表。 4. 經濟部水利處，1999b，台灣地區日雨量年報 表。 5. 經濟部水利處，1999c，台灣地區日流量年報 表。

6. Abamowitz M and Stegun IA (1972) Handbook of mathematical functions with formulas, graphs, and mathematical tables. Dover Publications, New York, 1043 pp

7. Bardet JP (1992) A viscoelastic model for the dynamics behavior of saturated porelastic solids. Journal of Applied Mechanics. 59:128-135 8. Berryman JG, Thigpen L, and Chin CY (1988)

Bulk elastic wave propagation in partially saturated porous solids. The Journal of the Acoustical Society of America. 84:360-373 9. Biot MA (1941) General theory of

tree-dimensional consolidation. Journal of Applied Mechanics. 12:155-164

10. Biot MA (1956a) Theory of propagation of elastic waves in a fluid-saturated porous solid, I.

low-frequency range. The Journal of the Acoustical Society of America. 28:168-178 11. Biot MA (1956b) Theory of propagation of elastic

waves in a fluid-saturated porous solid, II. high-frequency range. The Journal of the Acoustical Society of America. 28:179-191 12. Coleman BD and Noll W (1963) The

thermodynamics of elastic materials with heat conduction and viscosity. Archive for Rational Mechanics and Analysis. 13:167-178

13. Cooper HH, Bredehoeft DJ, Papadopulos IS, and Bennett RR (1965) The response of well aquifer systems to seismic waves. Journal of Geophysical Research. 70:3915-3926

14. Drew DA and Passman SL (1999) Theory of multicomponent fluids. Springer Verlag New York Inc., New York, 308 pp

15. Eringen AC (1980) Mechanics of continua. R. E. Krieger Pub. Co., Huntington, New York, 592 pp 16. Farrell OJ and Ross B (1971) Solved problems:

Gamma and Beta functions, Lendre polynomials, Bessel functions. Dover Publications, New York, 410 pp

17. Grecksch G, Roth F, and Kumpel HJ (1999) Cosesmic well level changes due to 1992 Roermond earthquake comparing to static deformation of half space solutions. Geophysical Journal International. 138:470-478

18. Green and Naghdi PM (1965) A dynamical theorem of interacting continua. International Journal of Engineering Science. 3:241

19. Hassanizadeh M and Gray WG (1979a) General conservation equations for multi-phase systems: 1. averaging procedure. Advances in Water

Resources. 2:131-144

20. Hassanizadeh M and Gray WG (1979b) General conservation equations for multi-phase systems: 2. mass, momenta, energy, and entropy equations. Advances in Water Resources. 2:191-203 21. Hassanizadeh M and Gray WG (1980) General

conservation equations for multi-phase systems: 3. constitutive theory for porous media flow.

Advances in Water Resources. 3:25-40 22. Hsieh PA, Bredehoeft JD, and Farr JM (1987)

Determination of aquifer transmissivity from earth-tide analysis. Water Resources Research. 78:1824-1832

23. Hsu SK (1998) Plan for a groundwater monitoring network in Taiwan. Hydrogeology Research. 6:405-415

24. Jeng DS and Lee TL (2001) Dynamics response of porous seabed to ocean waves. Computers and Geotechnics. 28:99-128

25. Lamb SH (1945) Hydrodynamics. Dover Publications, New York, 738 pp

26. Ma KF, Lee CT, and Tsai YB (1999) The Chi-Chi, Taiwan earthquake: large surface displacement on an inland thrust fault. Eos (Transactions,

American Geophysical Union). 80: 605-611 27. Ohno M and Wakita H (1997) A water well

sensitive to seismic waves. Geophysical Research Letters. 24:691-694

28. Prevost JH (1982) Nonlinear transient phenomena in saturated porous media. Computer Methods in Applied Mechanics and Engineering. 20:3-18 29. Roeloffs EA (1998) Persistent water level changes

in a well near ParkField, California, due to local and distant earthquake. Journal of Geophysical Research. 103:869-889

30. Roeloffs EA (1996) Poroelastic techniques in the study of earthquake-related hydrologic

phenomena. In: Dmowska R and Saltzman B (eds) Advances in Geophysics, Academic Press, New York, Vol. 37, pp 135-195

31. Wang CC (1970a) A new presentation theorem for isotropic functions, part I. Archive for Rational Mechanics and Analysis. 36:166-197 32. Wang CC (1970a) A new presentation theorem

for isotropic functions, part II. Archive for Rational Mechanics and Analysis. 36:198-233 33. Yu SB, Kuo LC, Hsu YJ, Su HH, Liu CC, Hou

CS, Lee JF, Lai TC, Liu CC, Liu CL, Tseng TF, Tsai CS, and Shin TC (2001) Preseismic deformation and coseismic displacement associated with the Chi-Chi, Taiwan, earthquake. Seismological Society of America Bulletin.(in press).

34. Zienkiewicz OC (1984) Dynamic behaviour of saturated porous media: the generalized Biot formulation and its numerical solution.

International Journal for Numerical and Analytical Methods in Geomechanics. 8:71-91

表一 套配模式水井位置及特性參數表 站號 距礫石與砂 石交界處東 西長度L_bc （公尺） 水井所在經 度（二度分 帶表示法） 水井所在緯 度（二度分 帶表示法） 儲水係 數 S ** 水力傳導係數 k （公尺/天） ** 東西兩站 相距距離 （公尺） C η （天） 芳草 2（2E1）* 0 185350 2624184 0.00080 12.355 田洋 2（2W1） 6600 178717 2624989 30.499 6682 0.19 好修 3（3E1） 450 194052 2656100 0.00075 26.870 文昌 3（3W1） 3500 190120 2656250 56.572 3934 0.64 溪湖 3（3E2） 2000 196133 2649778 0.00114 52.445 趙甲 3（3W2） 6000 187624 2648441 67.133 8613 0.00025 九隆 3（3E3） 680 191168 2627781 0.00025 29.894 田洋 3（3W3） 12000 178717 2624989 11.568 12746 7.2 *受非拘限含水層垂直流況影響，退水趨勢不符合式（27）一維之假設，唯可將其視為 式（28）之邊界條件。 **水井抽水試驗[經濟部水利處（1999a）]。

︹ ︹ ︹ ︹ ︹ ︹ ︹ ︹ ︹ ︹ ︹ ︹ ︹ ︹ ︹ ︹ ︹ ︹ ︹ ︹ ︹ ︹ ︺ ︺ ︺ ︺ ︺ ︺ ︺ ︺ ︺ ︺ ︺ ︺ ︺ ︺ ︺ ︺ ︺ ︺ ︺ ︺ ︺ ︺ ← ← ← ←← ←_← ← ← ← ← ← ←← ← ← ← 27,500±300 6,500 90± 1,356±50﹖9,050±100 10,900±50 9,390±60 6,100±70 34,260±390 46,580±1,500 ＞50,000 ＞50,000 ＞50,000 ＞50,000 27,600±160 17,950±13016,000±370 6,290±90 6,240±70﹖ 10,700±1308,870±130 B3 B2 F3 T2 B1 F2 T1 F1 F1 T1 F2 T2 F3 T3 B3 B2 B1 T3 A9 A8 A7 A6 A5 A4 A3 A2 A1 100 100 50 50 0 0 -50 -50 -100 -100 -150 -150 -200 -200 -250 -250 -300 TF -300 --B---350 -350 0 5 公 里 10 15 礫石層 含水層（Aquifer） 水文地質概念分層界線 鑽探井 濾料深度範圍 濾管深度範圍 阻水層（Aquitard） 砂石層 定年結果（年前） 細砂層 粉砂、泥及黏土層 圖 例 公 尺 公 尺 ← 37,400±700 ← ← ← ←2,360±40 5,250±50 1,690±40 ← ←＞50,000 ← ←5,470±505,420±50 ← 圖一 濁水溪沖積扇分層概念圖 圖二 拘限含水層細顆粒砂石沈積物所在區域受震動時之概念模型 A1 A2 A3 A4 A5 A6 A7 A8 A9

C

含水層尖滅端視為 應變速率為 0 震波傳遞方向 x 座標 x=L x=0 η 拘限含水層礫石沈積物 阻水層 阻水層 拘限含水層砂石沈積物

圖三 拘限含水層中砂石沈積物液相位移量隨時間變化圖，式（19） 假設 L =30000 公尺、v =1500 公尺每秒、_f η C=20 天、 C Fv_f =1 公尺 0 0.2 0.4 0.6 0.8 1 0 5000 10000 15000 20000 25000 30000 X（公尺） 位移量（公尺） T=05 T=10 T=15 T=20 T=30 T=50 T=100 時間單位：秒

A B 0m 10000m20000m ★ 震央 斷層位置 河川 海岸線 研究區域 觀測井 彰化斷層 車籠埔斷層 大茅埔 - 雙冬斷層 台灣海峽 9 月 2 6 日地震震央 9 月 2 1 日地震震央 水裡坑斷層 八卦山脈 斗六丘陵 A B 圖四 研究區域地理環境圖 Quaternary Sequence Pliocene-Miocene Sequence Miocene Slate Paleogene Slate

+ 觀測井 2.50 9 月 21 日（地震後）日平均 水 位與 9 月 20 日（地震前）日 平均水位差（公尺）[水利處 （1999）] 含水層上方絕水層位置分布 含水層中粒石與砂石沈積物 介面位置分布[中央地調所 （1999）] 含水層尖滅位置 土壤液化區域 2W1 套配水井位置 3E1 3W1 3E2 3W2 3E3 3W3 Y X Y 2E1 2W1 X Y 0.00 15.00 30.00 45.00 60.00 75.00 90.00 105.00 120.00 135.00 150.00 0m 20000m 40000m （a）含水層一 （b）含水層二 （c）含水層三 k值分布（公尺/天） 圖五 含水層水位變化與水文地質關係圖 A B1 B2 B3 C D

-5.5 -4.5 -3.5 -2.5 -1.5 -0.5 0.5 1.5 0 10 20 30 40 50 時間（日） 水位高程（公尺） 圖六 3E1、3W1 水井水位模擬與現地資料比較圖 -+2.0 +0.0 +2.0 +4.0 +6.0 +8.0 +10.0 0 10 20 30 40 50 時間（天） 水位高程（公尺） 圖七 3E2、3W2 水井水位模擬與現地資料比較圖 + 觀測值

－

模式模擬值 3E1 3W1 （9 月 20 日） （11 月 9 日） 9 月 26 日，ML =6.9 9 月 21 日，ML =7.3 － + 觀測值

－

模式模擬值 3E2 3W2 （9 月 20 日） （11 月 9 日） 9 月 26 日，ML =6.9 9 月 21 日，ML =7.3

0 5000 10000 15000 20000 25000 30000 35000 1 91 時間（天） 自強大橋站日流量 （立方公尺/秒） 0 10 20 30 40 50 西螺站日雨量（公釐） 圖八 水井 X 下游河道流量突升圖 9 月 26 日，ML =6.9 9 月 21 日，ML =7.3 (9 月 1 日) (11 月 30 日) 雨量站 p流量站

-5 0 5 10 15 20 0 10 20 30 40 50 時間（天） 水位高程（公尺） 圖九 3E3、3W3、2E1 及 2W1 水井水位模擬與現地資料比較圖       C Ln η （天） 圖十 _      C Ln η 與 1000 bc kL 關係回歸圖 y = -1.6401x + 2.9431 R2 = 0.9831 -4 -3 -2 -1 0 1 2 1 2 3 4 5 1000 bc kL ︵平方公尺 \ 天︶ 3W1 2W1 3W2 3W3 + 觀測值

－

模式模擬值 3E3 3W3 （9 月 20 日） （11 月 9 日） 9 月 26 日，ML =6.9 9 月 21 日，ML =7.3 2E1 2W2