一、建立以資料同化為基礎之區域地下水井抽水量推估模式 (一) 地下水模擬模式

第肆章區域地下水抽水量推估模式建立

一、建立以資料同化為基礎之區域地下水井抽水量推估模式

(一) 地下水模擬模式

一個三維、飽和、拘限的地下水模擬模式可以下式表示(Bear, 1988; Willis and Yeh, 1987)：

[ ] ( , ) of three-dimensional spatial coordinates)， t 是時間 (T) (time)， K 是水力傳導系數張量(the hydraulic conductivity tensor) (LT ), and ^q 是控制方程式中的源項(T ) (the source term)或沉項(T ) (the sink term)。地下水井的抽水或補注在地下含水層所造成的效應可以用各井(wells)或各井的集合(well clusters)來表示。因此，

在控制方程式中，各井或各井的集合所造成的沉項可表示成：

在本計畫採用美國地質調查研究所(U. S. Geological Survey) 所發展的 MODFLOW (Harbaugh et al., 2000)做為地下水模擬模式，用以模擬地下水流況、地下水水位分布、地下水補注及抽水所造成的影響。

因為MODFLOW 採用有限差分法(finite difference)進行空間

上的離散化，採用一階背向式尤拉有限差分(first order backward Euler finite difference)進行時間上的離散化，因此在求解地下水方程式(4-1)時，其方程式可改寫為：時所採用的單位時間長度，而單位受力時段(stress period)，其維度為(T)，其定義為在一受力時段中源項或沉項^q^k1為定值的單位式中為外加之沉項(sink term)，這個沉項在納進法中被視為造成模式結果與觀測結果之間落差的原因，即為未知的地下水抽水量。

一個觀測井的水位，可能被週邊多個抽水井的抽水行為所影響，因此，一個觀測水位數據可能含有多個抽水井的抽水訊息。

所以，在推估地下水抽水量時，需將觀測水位數據所含有各抽水井的抽水訊息拆解出來，並還原、分配回到各抽水井。一般來說，

觀測井的數量比抽水井的數量少，因此，傳統的最佳化 (optimization)技術可能無法提供最佳的抽水量推估，因為資料不足會造成不適定的問題(ill-posed problem) (Yeh, 1986)。

納進法直接將觀測資料轉換成抽水量並導入控制方程式中，

因此，地下水控制方程式可寫為(Auroux and Blum, 2008)：

Ph_k1 h_k

t  Lh_k1 q_k_1 G h



^o Ch_k1



^(4-4)

式中 G (L /T) 為一(ⁿ^ⁿ^o)矩陣，亦稱為穫得矩陣 (the gain matrix)； ⁿ^o 為觀測資料的數量，儲存於具有ⁿ^o維度的向量^h^o之中；^C (無單位)是一個可將(ⁿ^ⁿ^o)的模式水位分布矩陣^h^k^¹投射至觀測空間(observation space)的矩陣。Hoke and Anthes (1976)指出穫得矩陣 G 為一純量(scalar)矩陣，此純量可用鬆弛的技術將模式輸出之狀態變數導向觀測數據^h^o。此外，Stauffer and Seaman (1990) 指出穫得矩陣 G 亦代表著不同數學過程的相對強度，以數值平衡的方式將不同的變數同化在一起。若穫得矩陣 G 具有較大的值，代表模式對觀測值放入較大的權重，因此，模擬及預測結果會以較快的速度接近觀測值。Vidard et al. (2003)以參數辨識的角度，嘗試用最佳化的技術得到最佳的穫得矩陣 G，這個結果顯示納進法實際上與卡門濾波法有相關，若假設 L 為一個對稱的正定矩陣，則^G ^^C^T^R^¹，其中 R 為觀測殘差的協方差矩陣 (covariance matrix)，因此方程式(4-4)則可被視為尤拉-拉格朗治方程式(Euler-Lagrange equation)，而且因為 L 為正定矩陣，所以具有充分而且必要的條件來解下列能量最小方程式(Auroux and Blum, 2008)：

min

h 1



h h_k



^T^{P h}



^{ h}_k



^^t₂ ^h^T^(Lh^{ 2q}k1) ^t₂ (h^o Ch)^TR^1(h^o Ch)

 

^(4-5)

因此，Quarteroni and Valli (1997)指出方程式(4-5)的解等同於方程式(4-4)的解，並用數值近似的方法得到將方程式(4-1)加入納進法沉項^{G h}



^o^{ Ch}^k1



_{的弱解。此外，}Li and Navon (2001)及 Vidard et al. (2003) 指出納進法亦可被視為四維變分法的一種 (four-dimensional variational data assimilation scheme)或是特殊條件下的卡門濾波。此，這些時間點亦稱為更新時間(update time)。

此處，以Ide et al., (1997)的著作為基礎，介紹資料同化常用 (innovation residual)並可表示如下：

o f

u u u

d h Ch . (4-6)

在地下水推估的問題中，運算子 C 是一個矩陣，矩陣內各

原素由 0 及 1 所組成，此矩陣乃將模式輸出的水位分布結果 (model output)投射至水位觀測空間(observation space)。最後，採用納進法可將控制方程式改寫成：

h_u1^a  F_u(h_u)G h



u1o Ch_u1^f



^{ F}^u^(h^u⁾^^Gd^u1 ^(4-7)

假設在方程式(4-4)中，時間^t^u 與時間^t^u1之間有 m 個計算時間長度，可能跨越數個受力時段(stress period)，運算子^F^u 及無因次矩陣 G (-)可由方程式(4-4)中推導而得到，並且包含了已透過資料同化技術更新後之沉項向量^q^u。在方程式(4-7)的右邊為預測 (forecast)，而左邊為分析(analysis)，雙邊皆屬於同一時間點^t^u1。因此，在實際執行上，方程式(4-7)是以時間序列的方式進行兩個步驟，第一個步驟為求解方程式(4-4)以得到水位分布預測值^h^u1^f ，此步驟不考慮資料同化技術納進法所得到在控制方程中的新增項。第二個步驟為根據水位(或水壓)觀測數據，得到殘差(residual)，

並以此進行地下水抽水量推估及求解方程式(4-7)以得到系統更新後的分析值(analysis) ^h^u1^a 。

在地下水推估的問題上，因觀測資料的特性，需採用觀測資料納進法(observation nudging)，在此方法之中納進法將觀測資料所得到之殘差(residual)以權重方式，分配到臨近的網格中，此分配可為空間上及時間上的分配，分配的範圍以其影響半徑(radius of influence)來定義。影響半徑亦可稱為影響範圍，影響範圍需透過敏感度分析，得到抽水行為對地下水水位變化造成的影響及其影響範圍，此影響範圍與地質條件及地下水模式單位時間長度 (stress period)有關。由納進法造成的向量^Gd^u1是一個在控制方程式中新增的沉項，因此，殘差可被視為因未估算的地下水抽水量而造成的累計地下水水位洩降(cumulative drawdown)，因此，資料同化技術納進法在此的功能即為將此洩降轉換成未知的地下

水抽水量，當資料同化技術納進法所推估的地下水抽水量放入模權重的方式(distance-weighted function)將洩降分配給位於觀測井的影響範圍內的各抽水井，因此^h^{i, j} ^{ }^w^{i, j}^d^u1,i，其中^w^^{i, j}為分配率或稱為分配係數(distribution ratio or interpolation coefficient)，

分配率的計算方式為以各觀測井及抽水井的座標，依距離計

需注意的是，各分配率均為非負之數值(non-negative value) 且總合需為1，亦即： method)(亦稱為反應係數法 response coefficient method)轉換成地下水抽水量( Becker and Yeh, 1972； Yeh, 1986)： coefficient)或 Jacobian 敏感係數(Jacobian sensitivity coefficient)，

此係數可透過敏感性分析而得到(Yeh, 1986)。若地下水系統反應為線性關係(或是非線性呈度不高)，則獲得矩陣（增益矩陣）在時間長度不變、抽水井(或井群)位置不變、觀測井位置不變、地質條件不變等各項條件相同下的前提下，是為定值(或可合理假設為定值)，而影響係數則用於把各抽水井所造成之水位殘差轉換成抽水量。^q^{i, j}即為透過觀測井的殘差^h^{i, j}推估的抽水井 ^j的抽水量增量(pumping rate increment)或抽水量減量(pumping rate decrement)。方程式(4-11)代表著每一觀測井水位觀測值與模式輸

量推估值，這些抽水量推估值需透過權重^w^{i, j}^* 將多個抽水量推估值整合成一個最終之抽水量推估值^^q^'^u^^1,^j，此整合過程可寫成下式：

1, , ,

'_u _j _{i j} _{i j}

q _ w q



 





(4-12) 抽水量推估^^q^'^u^^1,^j是一個有ⁿ^w維度的向量，需被投射回原地下水系統的狀態變數空間(state space)，而得到方程式(4-4)中的納進項的向量

1 ^T ' 1

u u

Gd _  gC q _ (4-13)

式中的純量g 即為前述之鬆弛倍數(gain relaxation time)。

抽水量^^q^{i j}^, 可由方程式(4-11)所計算。

(nudging term)僅於更新時間點有效。

i i

圖 4 -2 以序列方法採用資料同化技術進行地下水抽水量推估之流程示意圖

以資料同化的納進法進行更新的流程如圖 4-2 所示，此例子假設每一受力時段(stress period)結束時可得到地下水水位觀測資料。圖中顯示於長方形圖塊中的符號 T 為位於兩個連續更新時間點(update time)^t^u及^t^u1的受力時段(stress period)的編號。更新時間^t^u1位於受力時段(stress period) T 結束的時間點，舉例來說，更新時間^t^u ^²剛好是受力時段^{T 2}^ 結束的時候。需特別說明的是，

每一個受力時段週期可能含有多個數值積分上的時段(time step)。

此處假設每一受力時段結束時可得到水位觀測資料(或是每一更新時間點均有水位觀測資料)。以資料同化的納進法進行地下水抽水量推估並更新狀態變數(地下水水位分布)及預測值的過程包含兩個步驟，第一個步驟是以資料同化的納進法以在每一更新時間點^t^u1取得之水位觀測資料，推估前一時段^T之地下水抽水量^Q^'^T 。第二個步驟是以推估之地下水抽水量^Q^'^T (即為控制方程式(4-4) 中新增加之沉項(sink term))，再執行一次地下水模擬模式，以得到分析值(analysis)。此分析值即為利用資料同化的納進法對更新

T=1 T=2 T=3 T=4 T=5

Q’0 Q’1 Q’2 Q’3 Q’4

Update at tu=1

Update at tu=2

Update at tu=3

Update at tu=4

Q’1 Q’2 Q’3 Q’4

(Q’₀ is the initial guess)

時間點^t^u1及對下個時段^{T 1}^ 之預測，進行更新之地下水水頭分布。

在實務上，上述步驟是於每一時段結束時以離線方式(offline)執行。資料同化的納進法亦可以即時(real time)的方式在水位觀測資料可得的時間點來執行系統更新及預測更新。以下是資料同化納進法完整的流程：

步驟 0：計算影響係數、分配係數及權重係數

分配係數及權重係數可依給定的影響範圍、觀測井座標及抽水井座標，分別用方程式(4-30) 及 (4-31) 計算。影響係數 (^r^{i j}^, ⁽^{ }^hⁱ^/^^q^j⁾)可用敏感度分析取得，如果系統為一非線性系統，

則影響係數需要在每一更新時間點同時進行影響係數的更新，如果系統可視為一線性系統，則影響係數為不變之常數。在拘限含水層中，影響係數可視為常數，且可以表示為：

, , 1/ , 1 , 2 / , 2 ... , / ,

u u u

i j i t j T i t j T i t e j T e

r  h _ q _  h _ q _   h _ q _ ( 4-32) 式中^e為規劃總時間長度內之受力時段之數量，而影響係數之定義為於一時段內，一單位抽水量於抽水井^j所造成於觀測井ⁱ 之水頭變化。影響係數反應出由抽水行為產生洩降錐所造成之影響大小及範圍。

步驟 1：預測觀測井位置下一更新時間之地下水水位(或水頭) (hydraulic head)

利用模擬模式、已知及推估的地下水抽水量進行地下水水位預測。一開始，若不知地下水抽水量，則可假設地下水抽水量為零。本計畫因無抽水量資料，故一開始假設抽水量為零。

步驟 2：計算位於觀測井的殘差(innovation residual)

可利用方程式(4-8)計算殘差，殘差即為在更新時間點位於觀

測井的位置，模式輸出之地下水水位與水位觀測得到之地下水水位的差值。

步驟 3：推估地下水抽水量

利用前述步驟所得之資訊及方程式(4-13)計算納進項，即可計算得到推估之地下水抽水量。

步驟 4：於更新時間點進行地下水水位的更新

此步驟以推估之抽水量透過模擬模式於更新時間點進行地下水水位的更新，亦即計算分析值(analysis)。

重覆執行步驟 1 至步驟 4，直到所有規劃時間長度內的時段都完成為止，若是即時(real time)的應用，則重覆執行步驟 1 至步驟 4，不需停止。以本計畫為例，推估 88 年至 99 年之月抽水量，則每一時段長度為 1 個月，規劃時間長度為 12 年。時間精度依需求及模擬模式之時間模擬精度而定。

在文檔中應用資料同化方法推估區域地下水利用之研究 (頁 66-79)

一、建立以資料同化為基礎之區域地下水井抽水量推估模式 (一) 地下水模擬模式

第肆章 區域地下水抽水量推估模式建立