國立臺灣大學工學院土木工程學系 博士論文
Department of Civil Engineering College of Engineering National Taiwan University
Doctoral Dissertation
利用伴隨狀態法耦合水流與傳輸 逆推地下水參數及未知條件
Coupling flow and transport for groundwater parameter and unknown condition identification using adjoint state method
劉宏仁 Liu, Hung-Jen
指導教授：李天浩、徐年盛 教授
Major Professor: Lee, Tim Hau, Hsu, Nien-Sheng
中華民國 98 年 6 月
Jun. 2009
誌謝
余致力於學術研究七載有餘，初學懵懂，幸賴恩師李天浩教授與徐年盛教授 辛勤指導，每遇壅塞之際給予學生明燈指引，於無涯學海中覓得一份天地，教授 們用心至深，萬分感謝。安孺、呈懋、文生、孝忠、德霖、智勇、國展諸位學長 姐，於研究過程中，無私的經驗分享與學思交流，一路引領與提攜著我踏上研究 的旅途，學長姐用意至深，萬分感謝。
摯友欣懋、宗翰、嘉興、俊宏，一路相伴、談天歡笑、悲喜與共，真誠的關 懷給予我無限的溫暖，使苦澀的青春染上繽紛色彩，摯友們用情至深，萬分感謝。
文明、佑誠、玉豪、柏凱、耀霖、威名、佑蓉等學弟妹，一同在學術的殿堂下，
教學相長、互相扶持，走過了這段研究所的歲月，學弟妹用誠至深，萬分感謝。
每一日，家人、表哥與玉樺給我最大的勉勵與鼓舞，伴我度過一個個春夏秋冬，
無論起落、無論得失，永遠無條件的支持，使本人永遠陽光進取、奮發向上，家 人們用愛至深，萬分感謝。
特別感謝陳主惠老師，對於繁雜的數學不厭其煩地幫助學生一一推導與檢 核，使得論文更加完備。最後感謝這些年來所有的研究伙伴們，能與你們相識是 我的榮幸。僅以此篇文章獻給各位。
摘要
本研究提出一個一般化逆向問題求解方法，結合了地下水模擬模式、伴隨狀 態變數法、梯度法與近似牛頓法，在最小誤差平方和的架構下，最佳化估計未知 參數、初始條件或邊界條件。應用伴隨狀態法推導所得之伴隨問題是一種對於參 數或條件估計錯誤所造成之模擬誤差的描述，藉由狀態變數與伴隨狀態變數的積 分，目標函數對應於所有未知數的梯度值可以快速獲得。
本研究以水平二維拘限含水層為考量對象，首先以一個設計過的示蹤劑試驗 探討水力傳導係數(K)檢定問題。除了水頭與一階動差觀測外，再加入了地下水 流速與示蹤劑零階動差的觀測。貢獻度指標分析結果指出流速僅包含觀測位置 K 的資訊，而水頭與動差則包含了大範圍內 K 值分佈的訊息。
本研究提出不同示蹤劑濃度釋放策略，營造一個隨空間變化的零階動差場，
解決零階動差二元分佈問題，且所得之零階動差其對於檢定 K 的貢獻度還高於 水頭與一階動差。參數檢定結果發現，水頭觀測資料描述的是地下水等勢能線的 分佈，此分佈直接反應參數場大小；動差觀測描述的是流線的分佈，可以表現出 流線與流線間 K 值差異與的資訊。
在完整逆向問題的探討上，本研究以一個試驗設計過的抽水實驗配合最佳化 演算法，同時估計蓄水係數、導水係數、初始水位、邊界水位與邊界流量。相關 性分析結果顯示蓄水係數與初始水位呈高度相關；導水係數與邊界水位、邊界流 量皆高度相關，這樣的相關性導致逆向問題成為非唯一解。
觀測貢獻度分析發現抽水實驗初期的非穩態洩降對於蓄水係數與初始條件 檢定貢獻度較高，但單一個水位觀測的貢獻度仍小於一；而抽水晚期的穩態水位 觀測，則是對於檢定導水係數與邊界條件特別有效。靠近邊界的觀測井對於邊界 條件特別敏感，而抽水井的洩降主要貢獻至參數檢定。在觀測充足的狀況下，逆 向演算法可有效融合上述這些資訊，使得檢定結果收斂至未知數的真值。
關鍵字：地下水模擬、水流與傳輸、參數檢定、完整逆向問題、伴隨狀態變數法、
試驗設計、觀測貢獻度、資料融合
Abstract
In this paper, groundwater simulation models, adjoint state method, gradient search method, and least square algorithm are combined to formulate an efficient optimization approach to solve the groundwater inverse problem. The adjoint state method is used to evaluate effectively the gradient of objective function with respect to parameter or condition.
The horizontal two-dimensional confined aquifer is the research target. First, the parameter identification is discussed through a designed tracer test. Head, flux, zeroth and first moment observations are utilized to estimate hydraulic conductivity(K) field.
The moment observations at different locations contain some indirect trend
information of K field. Next to head observations, they provide additional knowledge useful to parameter identification in groundwater system. Using all three kinds together, the case study demonstrates to elevate the efficiency and accuracy of solution substantially.
Second, the complete inverse problem is taken into consideration. A designed pumping test is performed to simultaneously estimate aquifer parameters, initial condition, and boundary conditions in groundwater modeling. The correlation analysis shows high connection between storage coefficient and initial condition.
Besides, transmissivity and boundary conditions are also highly correlated. A time series of unsteady head is needed for storage coefficient and initial condition estimation. The observation near boundary is very effective on boundary condition estimation. The observation at pumping well mostly contributes to the estimation of transmissivity.
Keywords: groundwater modeling, flow and transport, parameter identification, complete inverse problem, adjoint state method, experimental design, contribution of observation, data assimilation
目錄
口試委員會審定書 … … … ...… … … ..… … … … ... I 誌謝 … … … ...… … … ..… … … … ... II 中文摘要 … … … ...… … … ..… … … ... Ⅲ 英文摘要 … … … ...… … … … ..… … … ... IV 目錄 … … … ...… … … . V 表目錄 … … … ..… … … ... VII 圖目錄 … … … ..… … … … ..… ...… ... VⅢ
第一章 緒論 … … … ..… … … .… … 1-1 1.1 問題概述和研究動機 … … … ..… … … .... 1-1 1.2 文獻回顧 … … … .… .… … 1-3 1.2.1 地下水觀測和試驗設計 … … … .… … … .… . 1-4 1.2.2 參數檢定 ...… … … .… .. 1-8 1.2.3 完整逆向問題 ...… … … ...… ..… .. 1-10 1.3 研究目標 … … … .… … … .… … 1-12 1.4 論文架構 … … … ..… … … 1-13
第二章 研究方法論 … … … ...… … … … 2-1 2.1 最佳化演算法 … … … ..… … 2-1 2.2 地下水參數檢定 … … … ..… … … 2-3 2.2.1 地下水水流與示蹤劑動差控制方程式 … … … .… … .… . 2-4 2.2.2 目標函數 ...… … … .… .. 2-9 2.2.3 伴隨問題 ...… … … ..… .. 2-10 2.2.4 參數更新方法 ...… … … .… … … ..… .. 2-13 2.3 完整逆向問題 … … … ...… … … ... 2-15 2.3.1 非穩態地下水水流控制方程式 … … … ...… … … .… ..… … … . 2-15 2.3.2 目標函數 ...… … … . 2-16 2.3.3 伴隨問題 ...… … … . 2-17
第三章 參數檢定之案例研究 … .… … … .. 3-1 3.1 試驗設計 … … ..… … … … .… … … .… … … .… . 3-1 3.1.1 案例說明 … … … .… … .… . 3-2 3.1.2 狀態變數場 ...… … … .… .. 3-4 3.1.3 初始估計誤差 ...… … … ..… .. 3-9 3.2 貢獻度分析 … … … ..… 3-12 3.3 伴隨狀態變數場 … … … … .… … … . 3-20 3.4 檢定結果分析… … … … .… … … ..… … . 3-24 3.4.1 參數最佳化估計 … … … .… … .… . 3-24 3.4.2 水位、流速與動差觀測的價值 ...… … … .… .. 3-27
第四章 完整逆向問題之案例研究 … ... 4-1 4.1 試驗設計 … … … ...… … … .… … … ... 4-1 4.2 敏感度分析 … … … .… … … … .… … … 4-3 4.3 相關性分析 … … … ...… … .… … … .. 4-4 4.4 貢獻度分析 … … … .… … … .… … … 4-6 4.5 可檢定性分析 ...… … .… … … .… … … .… 4-8
第五章 總結與建議 … … … ... 5-1 5.1 伴隨狀態法 … … … .… … … . 5-1 5.2 試驗設計 … … … .… … … . 5-2 5.3 參數檢定 … … … ...… … … .… … … .. 5-3 5.4 完整逆向問題 … … … ...… … … .… … … .. 5-5
參考文獻 ..… … … ..… … … . 參-1 附錄 A 地下水參數檢定之伴隨問題推導 … … … . A-1 附錄 B 參數檢定問題之實際求解過程 … … … ...… … … ..… … B-1 附錄 C 完整逆向問題之伴隨問題推導 … … … ...… … … ..… … C-1 簡歷
表目錄
表 3.1 含水層真實參數場 … … … … .… … … ..… … … .. 3-3 表 3.2 OSSE 觀測資料 … … … ...… ..… … … .… … … … . 3-12 表 3.3 觀測精度與要求之參數精度 … … … ..… … ...… 3-14 表 3.4 水位、流速與動差觀測對於檢定 K1的貢獻度指標 … … … ...… .3-18 表 3.5 水位、流速與動差觀測對於檢定 K2的貢獻度指標 … … … ...… .3-19 表 3.6 不同觀測資料組合下之最佳參數估計結果 … … … ... 3-25 表 3.7 單一觀測資料下之參數最佳估計值結果 … … … ..… … … ..… . 3-28 表 4.1 參數、初始與邊界條件之相關性矩陣 … … … ...… … … ... 4-5 表 4.2 觀測精度與要求之參數、初始與邊界條件精度 … … … ...… .. 4-6 表 4.3 不同觀測組合下之最佳估計結果 … … … ..… … … .... 4-8 表 4.4 不同未知數設定下之最佳估計結果 … … … .… … … ...… . 4-9
圖目錄
圖 2.1 地下水一般化逆向問題求解流程圖 .… … … .… … … ...… ... 2-2 圖 2.2 參數最佳化搜尋過程圖 .… … .… … … ...… … … ... 2-15 圖 3.1 虛擬拘限含水層俯視圖 … … … .… … … 3-3 圖 3.2 真實地下水水位與流速分布場 ... 3-5 圖 3.3 真實零階動差分布場 ..… … … ..… . 3-6 圖 3.4 真實一階動差分布場 … … … .… … … ... 3-8 圖 3.5 初始參數估計所造成之水頭與流速誤差分布圖 … ..… … … .… .. 3-9 圖 3.6 初始參數估計所造成之零階動差誤差分布圖 … … … ..… ... 3-10 圖 3.7 初始參數估計所造成之一階動差誤差分布圖 … … … ... 3-11 圖 3.8 觀測 H(60,15)對於檢定各個計算節點 K 值的貢獻度分布圖 … … … . 3-14 圖 3.9 觀測 q(60,45)對於檢定各個計算節點 K 值的貢獻度分布圖 … … … .. 3-15 圖 3.10 觀測 m0(30,45)對於檢定各個計算節點 K 值的貢獻度分布圖 … … 3-16 圖 3.11 觀測 m1(90,15)對於檢定各個計算節點 K 值的貢獻度分布圖 … … .. 3-16 圖 3.12 伴隨一階動差分布圖, *E10 … … … ...… … … . 3-213 圖 3.13 伴隨零階動差分布圖, … … … ...… … … . 3-222 圖 3.14 伴隨水頭分布圖, … … … ...… … … . 3-231 圖 3.15 伴隨流速分布圖, … … … .… ...… … … .3-23 圖 3.16 參數最佳化過程中不同觀測資料組合下誤差平方和的遞減圖 … ... 3-25 圖 3.17 參數最佳化過程中單一觀測資料下誤差平方和的遞減圖 … … ..… . 3-27 圖 4.1 虛擬二維拘限含水層示意圖 … … … ...… … … .. 4-2 圖 4.2 抽水實驗過程之抽水井與監測井之洩降曲線圖 … … ...… … … ..… … 4-3 圖 4.3 水位 h(30,20)對於參數、初始與邊界條件的敏感度隨時間變化圖 ...4-3 圖 4.4 水位 h(10,20)對於參數、初始與邊界條件檢定之貢獻度隨時間變化圖
… … … ... 4-6 圖 4.5 水位 h(30,20)對於參數、初始與邊界條件檢定之貢獻度隨時間變化圖
… … … ... 4-7 圖 4.6 水位 h(50,20)對於參數、初始與邊界條件檢定之貢獻度隨時間變化圖
… … … ... 4-7
第一章 緒論 1.1 問題概述和研究動機
土壤與地下水受到有機物污染的問題近年來在台灣受到關注，如美國無線電 公司（RCA）有機溶劑污染地下水，中油高雄煉油廠大型油槽及台南縣永華等多 個加油站儲油槽洩漏污染事件等。汽油、柴油等石化油品以及三氯乙烯（TCE）、
四氯乙烯（PCE）等有機溶劑均屬於微溶於水的非水相液體（non-aqueous phase liquid，簡稱 NAPL），進入到地下含水層中，主要是以未溶解的殘留相（residual）、 自由相（free product）與污染池（NAPL pool）存在含水層中。由於 NAPL 在地 下水系統中的移動，受兩液相的介面力主控，抽出殘留相或自由相，可能需要垂 直:水平達 100:1 的水力梯度，現場實際上不可能實施。過去採用抽出處理（pump and treat）NAPL 污染場址的經驗，停抽後地下水中 NAPL 濃度會回彈（rebound）； 近年許多案例改採圍堵和自然衰減（Natural Attenuation）的方法控制。
地下水污染整治或控制的前提，是透過場址特徵調查（site characterization）， 掌握水文地質與污染分布。場址特徵包括：
一、 水文地質參數：地質特徵參數包括含水層地質剖面，如土壤組成、地質分 層等；以及水力傳導係數（hydraulic conductivity）、蓄水係數（storage coefficient）與延散係數（dispersion coefficient）等。
二、 水流特徵：地下水流動特徵包括地下水水位分布、流向、流速與流量、水 文初始條件和邊界條件等。
三、 污染特徵：污染物分布特徵包括液態污染物位置、污染總量、溶解相的污 染範圍，污染物的物理化學參數，如吸附係數（adsorption coefficient）、衰 減係數（decay coefficient）等。
常用的場址特徵調查方法有鑽井取樣（core sample）、抽水實驗（pumping test）
與示蹤劑試驗（tracer test）等。近年不同的井中流速量測儀（borehole flow meter）
技術趨於成熟，可以在水井內直接量測「井中流速」（未必等於地下水流速）。示 蹤劑試驗是在水井中直接注入示蹤劑，或在注入井中隨水注入示蹤劑，使其隨地
下水流傳輸；在觀測井處觀察其濃度隨時間變化情形，此濃度歷線又稱為突破曲 線（breakthrough curve，簡稱 BTC）。示蹤劑試驗過去多用來估計地下水的延散 係數，進年來井際分溶示蹤劑試驗（partitioning inter-well tracer test，簡稱 PITT，
Jin et al. 1995）則應用於調查污染場址的 NAPL 總量與分布。試驗方法是利用抽 水井與注水井營造一個局部穩態的地下水流場，於注水井注入多種守恆示蹤劑
（conservative tracer）與不同分溶係數的示蹤劑（partitioning tracer），取得觀測 井和抽水井的突破曲線資料，對於開始注入時間（t0）積分為各階動差
（moments），可提供場址污染物分布資訊。
由於 NAPL 分布不均、地文參數空間異質性、水流特徵因地文參數而變；
以及上述場址特徵調查方法，或為在平面上的點觀測（如鑽井取樣），或為參數 的區域平均值（如抽水試驗），無法透過取樣直接獲得全域、高解析度的特徵值。
所以場址特徵調查，除了必須從基本背景資料著手，透過試驗設計（experimental design）取得富含資訊的（informative）的調查、觀測資料外；還須配合地下水 模式（groundwater model），逆向反演各種場址特徵參數。反演場址特徵參數問 題可以分為三類：
一、 參數檢定（parameter identification）：「水文地質參數」和「汙染特徵參數」
不確定，透過遞迴調整參數，使得模式反應逼近調查數據的方式，反演估 計未知參數。
二、 資料同化（data assimilation）：「水文地質參數」已知，利用地下水模式反 演估計「水流特徵」中不確定的水文初始條件與邊界條件。
三、 完整逆向問題（complete inverse problem）：「水文地質參數」、「水文初始 條件」與「邊界條件」同時具有不確定性，利用地下水模式同步遞迴反演 推估。
傳統上對於「試驗設計」的認知，是在有限的預算內、對於反演最有助益的 空間位置、最短的試驗觀測時間，取得最大化場址特徵值的調查資訊。唯示蹤劑 試驗和井中流速等新類型調查資料和配置方法，對於反演場址特徵參數的效益，
並不明朗。James et al. (2000)在四周打入鋼板樁形成不透水條件的 PITT 現場實 驗場址，場址大小為 4.4*3.5*1.7m3，為小範圍的污染場址，比較平行和輻射兩 種注入井與抽水井的配置，發現所得的示蹤劑動差觀測對參數之敏感度不同。除 了井配置外，不同掃瞄範圍、邊界條件和回收條件，都有可能使示蹤劑動差對反 演的效益不同，但文獻中可供「試驗設計」參考的方法和案例研究極少；配合反 演方法論，完整呈現調查資訊效益的研究付之闕如。
伴隨狀態法（adjoint state method，簡稱 ASM）自 1990 年由 Sun 和 Yeh 提 出以後，陸續被應用於水文地質參數檢定問題的求解過程，取代「微小擾動法」
（perturbation method）計算賈克比矩陣中目標函數對於參數的梯度。相同論文 亦提出利用狀態變數與伴隨狀態變數，計算某個別觀測對於某檢定參數敏感度 的方法；以及利用參數自然變異的標準偏差，和該觀測誤差標準偏差作標準化 後的「貢獻度」（contribution of observation）量化指標，來評估該觀測對於檢定 不同參數的資訊有效程度。唯文獻中並無呈現伴隨狀態變數對於不同觀測值「貢 獻度」影響，並將其與試驗設計結合的研究。
最後，在參數檢定問題中，通常假設只有參數場為未知。實際狀況下，場址 的邊界條件未必能完全掌握；若場址範圍較大，初始條件與範圍內的補注量（或 滲漏量）和抽水量亦可能不確定，這些水文條件不確定性會影響參數檢定結果，
理論上參數檢定時必須同時考量參數和水文條件的不確定性，但目前仍未見到求 解完整逆向問題的研究。
1.2 文獻回顧
本節對於前人研究進行文獻回顧，1.2.1 小節說明現有地下水觀測與試驗設 計，1.2.2 小節描述在地下水模擬中對於參數檢定問題，已有的研究方法與成果，
由回顧中發現更進一步的問題；1.2.3 小節描述地下水完整逆向問題求解至目前 為止的成果，以及後續發展的方向。
1.2.1 地下水觀測和試驗設計
透過對於地下水系統的觀測與試驗設計，可以得到場址特徵資訊。早期的現 場實驗方法主要有抽水實驗（pumping test）與微水實驗（slug test），實驗方法是 透過抽水使得場址產生一水力梯度（hydraulic gradient），此水位差異造成了地下 水流動，於抽水井或觀測井觀測水位的的變化，利用地下水流模式與達西方程式 可估計地層的水文地質參數，這樣的求解方法在 Yeh (1986)的文章中有一個完整 的回顧，所有研究結果皆顯示了水位觀測包含了參數分布場資訊，所以可用於正 確檢定水力傳導係數（hydraulic conductivity，K）。
Sun and Yeh (1990)除了使用水頭觀測外，同時引入了示蹤劑試驗（tracer test），從上游注入井注入示蹤劑後於下游觀測井觀測流經示蹤劑之濃度
（concentration），將此濃度歷線作為觀測資料，並結合了地下水水流與示蹤劑 傳輸兩個偏微分方程，求解逆向問題，估計水力傳導係數的分布。在他們的兩個 K 分區的虛擬案例中發現，濃度觀測資料對於水力傳導係數的敏感度高於水頭觀 測，此結果顯示了濃度觀測資料在應用於求解地下水逆向問題時，提供了很有意 義的資訊。Wagner (1992)、Cheng and Yeh (1992)、Barlebo et al. (1998)、以及 Mayer and Huang (1999)的研究結果也支持了這一項發現，據此他們建議求解逆向問題 時同時使用水頭與濃度的觀測資料，可以達到較佳的檢定結果。
但是濃度的觀測資料是一個時間序列（time series），在數值模擬與逆向演算 上較為耗時，且數值正確性（accuracy）易受到離散處理方法的影響，所以後續 有人提出將分穩態濃度歷線積分為穩態動差（moment）的概念。此概念最早由 Jin et al. (1995)所提出，但一直發展至 Cirpka and Kitanidis (2000)才正式將示蹤劑 動差應用於水力傳導係數的檢定問題中。Jin et al. (1995)針對受 NAPL 污染之場 址提出了井際分溶示蹤劑試驗（partitioning inter-well tracer test，簡稱 PITT），PITT 是一種非破壞性間接量測方法，其試驗方法是利用抽水、注水井形成穩態地下水 流場，由注水井處釋放至少一種守恆與一種分溶示蹤劑以掃瞄井間的含水層，於 觀測井與抽水井處觀測並記錄示蹤劑濃度變化歷線或又稱為突破曲線
（Breakthrough Curve，簡稱 BTC）。守恆示蹤劑不與 NAPL 反應，而分溶示蹤劑 則會分溶進入 NAPL 中造成傳輸延遲（retardation），Jin et al. (1995)將突破曲線 對時間取零階及一階動差（zeroth and first moments）後，計算示蹤劑的平均到達 時間（mean arrival time），據以估計場址中之平均 NAPL 殘餘量多寡，到達時間 越晚，表示含水層中的污染量越大。
而 Harvey and Gorelick (1995)進一步直接將非穩態的溶質傳輸方程式
（unsteady-state solute transport equation）對整個時間域做積分，成為穩態的示蹤 劑動差方程式（steady-state tracer moment equation）。依其對時間積分的數量級大 小，對濃度積分所得為零階動差；對濃度乘上時間的一次方積分為一階動差；對 濃度乘上時間的二次方積分後為二階動差，依此類推。如此積分所得之動差僅在 空間中進行變化，而無時間上的變化問題，此舉將示蹤劑傳輸問題減低一個變化 維度，更容易求解。他們同時也說明一條完整的突破曲線上無限多個的時間觀測 資料可以由其前四個動差，即零階、一階、二階與三階動差（0thto 3rdmoment）， 代表了大部分的資訊。零階動差代表了整個突破曲線下方所包含的面積大小，即 為累積通過的示蹤劑質量除以觀測位置的流通量；一階動差代表了一個濃度權重 過後的示蹤劑到達時間，將其除以零階動差即等於抵達觀測井位的示蹤劑平均到 達時間；二階動差與三階動差則分別代表了突破曲線的開展（spread）與偏斜
（skewness）程度，此四個動差適切地描述了整條突破曲線。Harvey and Gorelick (1995)推導出示蹤劑動差方程式，使得後來的研究可以直接使用穩態的動差，解 決了非穩態的濃度時間歷線難以使用於逆向問題求解的缺點。
Rubin et al. (1997)推導了一個動差方程式的半解析解（quasi-analytical solution），可以快速計算瞬間吸脫附的溶質其在異質性（heterogeneous）地層中 的傳輸與動差分布情形。他們的研究結果發現，平均的示蹤劑傳輸速度約等於全 域流場流速的幾何平均值（geometric mean）。由於一階動差除上零階動差即代表 了示蹤劑的平均到達時間，也代表了示蹤劑傳輸速度的快慢，據此他們建議取一 階動差的觀測資料來估計整個水力傳導係數場的平均值，這是第一篇說明了濃度
積分後的動差包含了有意義的水力傳導係數場資訊的文章。
James et al. (2000)採用 Jin et al. (1995)的 PITT 試驗方法，場址四周打入鋼板 樁形成不透水封閉流場，於場址中央注入井注入多種示蹤劑，於四周圍空間均勻 分布的抽水井處完整回收（recover）所有的示蹤劑，記錄其突破曲線後取積分成 為動差觀測資料。利用穩態水流與 Harvey and Gorelick (1995)所提出之動差方程 式作為模擬模式，取水頭與一階動差作為觀測資料，逆向檢定三維含水層之延遲 係數（retardation factor，R）分布場。延遲係數與 NAPL 殘餘飽和度為一對一函 數，藉此推估得 NAPL 殘留量的空間分布。該研究宣稱其方法也可推估含水層 水力傳導係數，唯文中並未呈現其成果。James et al. (2000)也探討了抽水、注水 井的設置對於參數檢定的影響，發現 PITT 最先以上游注入下游抽出示蹤劑，對 含水層進行水平掃瞄，得到的動差觀測資料對於延遲係數的敏感度比較低。而若 改以中央注入四周圍抽出示蹤劑的輻射狀掃瞄，可以得到較為敏感的觀測。所以 他們由敏感度分析（sensitivity analysis）結果，建議採用輻射狀的抽水、注水井 配置。
Cirpka and Kitanidis (2000)是第一個將地下水水流與示蹤劑動差方程式公式 化（formularize）為一個地下水反向問題，以最佳化求解水力傳導係數場的研究。
其實驗方法採用傳統的示蹤劑試驗，在穩態流場中，於上游一水井中注入示蹤 劑，於下游多口觀測井處觀測。他們利用地下水水流模式結合示蹤劑動差模式正 向模擬示蹤劑於含水層中的傳輸情形，逆向問題則與 James et al. (2000)一樣同時 採用了水頭與一階動差作為觀測資料，反演水力傳導係數場的分布。結果發現一 階動差對於觀測位置上游一狹窄區段內的水力傳導係數場很敏感，此區段內參數 分布的細部特徵可以藉由引入一階動差觀測的方式而加以辨識。但是如果僅僅使 用動差觀測而無水頭觀測資料，逆向求解的結果將不穩定（unstable），這表示 水頭觀測仍是必要的限制條件（constraint），而示蹤劑的觀測資料為補充
（complemental）條件。研究結果也顯示當兩種不一樣的觀測資料同時使用時，
由於資訊的增加使得參數檢定的結果比單獨使用者為佳。
在 James et al. (2000)以及 Cirpka and Kitanidis (2000)求解參數檢定問題時，
發現其試驗設計使得在示蹤劑流經區域零階動差為固定常數，在示蹤劑未流經區 域零階動差為零的二元（binary）狀態；這樣斷然的、不連貫的（abrupt）分布，
使得後續在計算目標函數對應於未知參數的梯度值時呈劇烈變化，導致零階動差 觀測無法應用在逆向問題求解中，Cirpka and Kitanidis (2000)的論文文末即建議 針對此點繼續進行研究探討。
另外 Cirpka and Kitanidis (2000)之試驗設計採用傳統的示蹤劑試驗方法，這 樣的試驗設計會導致示蹤劑僅掃瞄了部分含水層，在非完全掃瞄（non-fully scanned）的狀況下，動差的觀測僅僅包含了掃瞄區域範圍內的資訊，將有限的 局部資訊用於推估全域的參數場時，形成了條件不足（ill-posed）的問題，這也 是另一個導致反演結果不穩定的原因。
Barth and Hill (2005)的研究則是指出突破曲線的動差是一種整合性的資 料，相較於單一的濃度觀測來說提供了一種更為強健的（robust）資訊，特別適 用於敏感度分析與參數檢定的問題。
除了上述普遍使用的水頭與濃度或動差的觀測外，對於地下水系統的參數分 布，我們還可透過地下水流速的測量予以瞭解，在固定的水力梯度下，地下水流 速越快，表示場址的水力傳導係數越高。傳統上流速觀測相較於水位與濃度觀測 而言較為困難，但隨著量測技術的發展，目前已有兩種主要可行的現地量測方 法：電磁流速儀（electromagnetic borehole flowmeter），可參考 Dinwiddie et al.
(1999)與被動流量計（passive flux meter），可參考 Basu et al. (2006)。實驗是以分 層抽水的方式，量測得到的各個不同深度的流速資料，以反演水力傳導係數的垂 直分布，反演過程假設每個深度的含水層在水平方向上為均質（homogeneous），
K 值僅在垂直方向上有所變化，每個不同的深度量測到的流速資料即代表該深度 的平均 K 值資訊。
Zijl (2004)首先嘗試同時使用水頭與流速觀測於逆向估計三維的導水係數分 布場，其反演方法採用直接求解參數法（direct approach），結果發現新加入的流 速觀測資料，使得檢定過程中融合（assimlate）了額外的資訊，參數估計結果明
顯變得穩定。然而其所給之虛擬案例依然假設整個地下水流場為均勻流（uniform flow），流速到處相同，所以他的研究結果與前人研究相同，都是將觀測流速置 於邊界條件上作為限制式使用，而非將場址內觀測得到的流速資料視為觀測資料 放入逆向演算邏輯中，所以地下水流速觀測資料中所包含的有用資訊仍然未被引 用到地下水逆向問題的求解中。
Wu et al. (2008)進行沙箱實驗與現地實驗，結果發現地下水流速儀量測所得 井中水流向與實際地下水流向相符，但量測所得之井中流速則為實際土壤中平均 流速的 3 至 6 倍，相關的研究也顯示了相同的結果，凸顯了地下水流速儀於現地 試驗應用的限制性，引用流速資料於參數檢定時需注意此問題。
綜上所述，地下水水頭、示蹤劑動差與地下水流速等這些可以觀測得到的資 料，都直接或間接包含了地下水含水層的資訊（information），應該同時被使用 於參數估計的問題中，以收最佳的反演結果。再者，Jawitz et al. (2003)發現零到 三階動差中又以零階動差與一階動差最具意義，但目前為止，在運移項主控、延 散程度很小的示蹤劑試驗中，僅有一階動差被應用於逆向反演水力傳導係數，而 零階動差則礙於其二元的分布型態，無法有效地應用，這是有待解決的問題。傳 統的示蹤劑試驗，只掃瞄了局部區域，產生了條件不足的問題，此點可以透過試 驗設計進行改善。井際分溶示蹤劑試驗即以抽水與注水井系統使得示蹤劑可以完 全掃瞄整個場址，並且採用空間均勻布置觀測井的方式，取得了足夠多的空間資 訊，此試驗方法值得做為參考。
1.2.2 參數檢定
對於地下水參數檢定問題的求解，除了給定適當的觀測資料外，還需有一套 有效率的演算邏輯，此邏輯中包含了水流與傳輸模擬以及參數的更新。Sun and Yeh (1990)回顧自 60 年代以來應用伴隨狀態法（adjoint state method，簡稱 ASM）
於地下水參數檢定之研究，歸納並提出了系統化的伴隨操作規則（adjoint operation rules），使得後續的研究可以依其規則很簡單且快速地由原始的模擬問
題（simulation problem）中推導出伴隨問題（adjoint problem），求解此伴隨問題 即可評估參數估計錯誤對於模擬所得狀態變數之誤差大小的影響。在應用上，只 要透過求解一次模擬問題與求解一次伴隨問題，目標函數對應於所有未知參數的 梯度值或敏感度就可以同時得到。相較於傳統的擾動法（perturbation method），
一次變化一個參數值的大小，以求取其梯度，當參數數量很多時，擾動法的計算 量會過於龐大，伴隨狀變數法解決了此一問題，大幅地提升了計算效率，James et al. (2000)以及 Cirpka and Kitanidis (2000)皆應用伴隨狀態法於其逆向演算程序 中。
James et al. (1997)利用序率方法（stochastic approach），在水文地質參數具有 定常性的假設下，利用擾動技巧，以貝氏推估法（Bayesian Estimator）計算突破 曲線之動差、地下水流速、延遲係數與水力傳導係數的共變異數（covariance）
及交互共變異數（cross-covariance）；再藉由空間中某些觀測點的突破曲線觀測 資料，推求三維含水層之延遲係數分布場。進一步，James et al. (2000)利用伴隨 狀態法改良其目標函數對各參數之敏感度矩陣計算，以空間點突破曲線動差和水 頭為狀態變數的方程式，配合高斯牛頓法求解三維地下水流與污染傳輸之反向問 題。其方法應用至一個模擬現場的 PITT 試驗，正確估計了的 NAPL 殘留的分布 趨勢，但其虛擬案例設定成四周為不透水的全封閉（closed）邊界，導致數學上 因為缺乏一個水頭參考值而無法計算。這問題在給定內部一參考點的水位之後，
獲得解決；但這個已知參考點於後續求解伴隨問題時，卻成為一個無限點源（point sink/source），誤差無法正確分布至伴隨狀態變數場，使得梯度值計算錯誤。這些 一方面是試驗設計不良的問題，另一方面也是對於伴隨狀態法原理的不瞭解所造 成的錯誤，需要由方法的原理以及伴隨狀態變數（adjoint state variables）的意義 出發進行探討與改善。
Cirpka and Kitanidis (2000)採用了伴隨狀態法及水流與示蹤劑動差模式於水 力傳導係數的估計問題上，利用 Kitanidis (1995)所提之半線性地質統計方法逆向 求解水力傳導係數場的分布，求解過程中不斷更新參數場，要求其模擬狀態變數
與實際觀測之水頭與一階動差相符。他們的參數場採用地質統計（geostatistics）
的觀念，呈對數常態（log-normal）分布，每一個計算格點上皆有一未知水力傳 導係數 K。逆向推求結果顯示，伴隨狀態法確實可應用於多參數求解問題中，快 速且正確地估計目標函數對應於大量未知參數的梯度值，反演所得之參數場可以 抓到整個真實分布趨勢以及部分的細部特徵。但由於其示蹤劑僅部分掃瞄了含水 層，非全域完整掃瞄，所以應用伴隨狀態法時，僅在局部幾個觀測點有誤差資訊 進入，導致估計所得之梯度值正確性降低。
Sun and Yeh (1990)提出伴隨狀態法的同時，也提出了一個「觀測貢獻度
（contribution of observation）」量化指標，此指標改良了敏感度的計算，先將敏 感度乘上參數解析度再除以觀測精度，等於將敏感度做了標準化的動作。若貢獻 度指標值大於一，表示在要求的參數解析度下，不同參數數值所造成的狀態變數 差異，可以被量測出來，如此該觀測即為有效觀測，可有效用於檢定參數場。此 貢獻度指標解決了敏感度分析易受參數本身數值大小影響，往往是參數數值越小 者敏感度越高的問題，提供了一個更客觀的評比工具，唯此一指標鮮少被引用於 參數檢定的問題中。
總結上述地下水參數檢定的研究成果發現，透過觀測資料的取得，配合地下 水模式與伴隨狀態法，確實可有效辨識含水層特徵與參數分布場。進一步的研究 應著重在嘗試將更多可觀測到的資料放入參數檢定流程中，配合最佳化演算法將 資料整合，以使得資訊最大化；最佳化演算法也必須探究其原理，以增加演算效 率，減少檢定的誤差。貢獻度這樣一個有意義的指標應該用來取代傳統的敏感 度，以辨識出真正有效的觀測資料，而且可以嘗試將貢獻度分析與上一小節之試 驗設計結合，使試驗設計更為完整。
1.2.3 完整逆向問題
完整的地下水逆向問題，除了參數場是未知需要檢定外，初始的地下水水位 分布與場址四周圍的邊界水位或流入、流出量也是未知的，同樣需要透過最佳化
程序尋找最適合的條件限制。少數幾篇文章探討了參數與邊界條件檢定問題，如 Rubin and Dagan (1988)利用序率方法（stochastic method）分析，發現靠近邊界的 水位受邊界條件影響很明顯。若邊界的位置或形式不確定時，他們建議採用開放 邊界的形式較為合理，也就是給定其為自由流通的邊界。LaVenue and Pickens (1992)先分析地下水模式對於何處的參數或邊界水位最為敏感，將這些點位指定 為控制點（pilot point），針對控制點進行參數與邊界水位檢定，而後其他未知參 數場與邊界水位再以克利金法（kriging approach）內插得到。Yeh and Mock (1996) 發展一個三階段的參數與邊界條件檢定方法，首先檢定一個全域平均參數值，然 後調整邊界條件以降低水頭模擬誤差，最後再微調局部的參數值，使其符合局部 地下水位觀測。
Gomez-Hernanez et al. (1997)以地質統計（geostatistic）結合最佳化方法估計 導水係數場與邊界條件。Bonilla and Cushman (2000)研究發現異質性含水層的邊 界條件對地下水位計算影響大。Hendricks Franssen and Gomez-Hernanez (2002) 以伴隨狀態法同時檢定了三維實際含水層之水力傳導係數與邊界水位。Brochu and Marcotte (2003)建議以觀測水位內插至邊界位置上，做為邊界水頭。Sun and Yeh (2007)提出了一個強健的（robust）參數檢定方法，結果顯示檢定所得之最佳 參數會隨著模式結構、事前資訊、邊界條件與模式應用目標不同而異。
除了邊界條件外，地下水模式的初始條件也是一個未知的分布場，在一般模 式模擬時，都是使用觀測水位內插後即假設為已知給定。Chen and Zhang (2006) 使用了卡門濾波（Kalman filter），在逐一時刻演算地下水位的過程中，不斷更新 上一時刻的水位。
回顧地下水完整逆向問題的文獻發現，額外考量了邊界條件與初始條件的不 確定性後，確實可讓檢定得到的模式比單純檢定參數者為佳。然而，初始條件的 影響，初始條件、邊界條件與參數間的關係（correlation）以及初始條件的檢定 問題皆鮮少有研究探討；更進一步，未知的初始條件、邊界條件與參數場應該同 時被檢定，這樣一個完整的逆向問題（complete inverse problem）求解到目前為
止還未有人提出有效的方法，而且對於這樣一個複雜問題的試驗設計也未有探 討，預料其將比參數檢定之試驗設計更為複雜，需想辦法取得有效的時間與空間 之觀測資料。
1.3 研究目標
本研究的目標有五點，分點敘述如下：
首先，發展一套可以融合（assimilate）水頭、流速、示蹤劑動差，計算效率 高，適用於逆向估計地下含水層水力傳導係數場的參數檢定方法。選擇包括地下 水守恆方程式、達西定率，和不同階示蹤劑動差方程式為控制方程式的模擬模 式，狀態變數包含水頭、流速、和各階示蹤劑動差等變數。本研究首創將地下水 流速、示蹤劑零階動差此二觀測，與水頭觀測、示蹤劑一階動差觀測，共同納入 目標函數值（objective function），探討這四種不同觀測對於地下水參數檢定的貢 獻。
其次，選擇伴隨狀態法，在穩態流場狀況下，推導包含地下水流速，假設水 力傳導係數為未知地下水參數的地下水流，和示蹤劑傳輸動差等變數的伴隨狀態 方程式；以及，在非穩態流況下，假設導水係數、蓄水係數等地下水參數，同時 初始、邊界兩水文條件亦不確定，推導在這些假設條件下的伴隨方程式。導出目 標函數對於各地下水參數、初始條件和邊界條件等未知數的梯度計算式，再結合 梯度搜尋法（gradient search method）建立逆向反演各不確定參數和條件的最佳 化演算方法。
第三，嘗試透過守恆示蹤劑注入、抽取、觀測井和邊界條件的空間配置，進 行示蹤劑試驗設計，使系統能夠達成五項目標：1.避免所得之示蹤劑零階動差為 二元化觀測，期望使零階動差能對檢定參數產生實質資訊；2.示蹤劑能完整掃瞄 整個場址，動差觀測資料可以表現全域所有參數的資訊，同時掃描範圍不隨地下 水參數改變；3.示蹤劑能完全面回收，避免從注入井流出示蹤劑的部份流管
（stream tube）掃描資訊流失；4.部分邊界條件為定水頭邊界條件，以滿足地下 水控制方程式的數學穩定性要求；5.部分示蹤劑邊界條件為定動差值條件，以滿 足示蹤劑動差控制方程式的數學穩定性要求。
第四，設計參數檢定的「觀測系統模擬試驗」（Observing System Simulation Experiment，簡稱 OSSE）案例，透過案例呈現伴隨狀態變數所隱含的意義與其 如何用於計算「貢獻度」指標，以利設計者理解不同場址調查井位安排下，水頭、
流速、示蹤劑動差等觀測，對於不同空間位置檢定參數的「貢獻度」。再利用不 同水頭、流速、示蹤劑動差組合，逆向反演地下水參數，驗證不同位置、各種觀 測資料對於參數檢定的實際影響。
第五，設計完整逆向問題的 OSSE 案例，透過案例分析參數、初始條件與邊 界條件間的相關性（correlation）；利用貢獻度分析尋找可有效用於檢定這些不同 型態未知數的時間與空間觀測資料，然後將這些觀測以不同的組合應用在完整的 地下水逆向問題求解中，對檢定結果進行分析，討論各個未知數的可檢定性以及 方法論的適用性。最後將所有可得的資訊同時使用於參數檢定、參數與邊界條件 檢定以及完整的逆向問題求解中，點出參數檢定將不確定的初始條件與邊界條件 假設為已知對檢定結果的影響，以及完整逆向問題的重要性。
1.4 論文架構
本論文架構與各章節重點如下：第二章針對地下水逆向問題的求解提出本研 究的最佳化方法論，結合模擬模式以及伴隨狀態法，將逆向問題公式化為一目標 函數最小化問題。後續將其應用在參數檢定問題以及完整逆向問題中，對所使用 的地下水模擬模式、目標函數、伴隨問題、梯度積分式以及參數更新方法做一完 整的描述。第三章是參數檢定問題的案例研究，藉由設計好的試驗，驗證最佳化 方法的可行性。首先詳細描述試驗設計的整體構想與作法，接著依反向求解過程 一一呈現狀態變數場分布、貢獻度指標分析、初始估計誤差、伴隨狀態變數場分 布及參數最佳化估計結果，詳細探討最佳化過程中每一個單元環節，並由最佳化
結果揭示各種觀測的價值。第四章是完整逆向問題的初探，以一個簡單的抽水實 驗配合上試驗設計，有效收集地下水洩降曲線資料。檢視未知參數、初始條件與 邊界條件對於逆向問題求解的影響，並探討這些未知數間彼此的關連性。應用貢 獻度指標於辨識有效觀測的時空位置，最後就完整逆向問題的求解結果，顯示出 本最佳化方法的適用性。第五章提出本研究過程中主要的發現和貢獻。
第二章 研究方法論
本章描述地下水逆向問題求解的方法論，2.1 節先清楚呈現逆向求解的程 序，包括未知數的給定、模擬模式與目標函數的計算，以及未知數的迭代更新。
2.2 節是將所提之最佳化方法應用於地下水參數檢定問題，詳細描述所使用的控 制方程式、邊界條件、目標函數、伴隨問題、梯度計算式與參數更新方法。2.3 節則是將研究方法擴展至完整的逆向問題求解，不僅檢定參數，同時檢定初始條 件與邊界條件，此節詳細描述控制方程式、初始條件與邊界條件、其對應之伴隨 問題與未知數之梯度計算式以及未知數之更新方法。
2.1 最佳化演算法（Optimization Algorithm）
地下水參數檢定問題的定義是模式參數未知，而有部分的狀態變數觀測資 料。反向問題的求解即在尋找最佳的參數場，使得模擬所得狀態變數場與觀測資 料相符。此搜尋過程是一個不斷更新（update）參數場的過程，更新的方法可分 為搜尋法（search method）、梯度法（gradient method）與擬牛頓法（quasi-Newton method），更新參數的目的在使得目標函數值，即模擬與觀測之誤差平方和最小。
由數學上來看，此三種更新方法分別為目標函數的零階、一階與二階近似方程 式，最佳化過程即在求使得目標函數為零的解，此解就是真實參數場。
而完整的地下水逆向問題除了上述的模式參數未知以外，初始條件與邊界條 件同樣也是未知數，需要以一個有效率的逆向求解程序估計之。於是針對地下水 逆向問題之求解，本研究提出一個一般化的（generalized）逆向問題求解方法，
此方法使用地下水模擬模式進行正向模擬，以模擬與觀測之誤差平方和為目標函 數，並應用伴隨狀態法推導所得之伴隨問題與梯度積分式計算目標函數對應於所 有未知參數與條件之梯度值，以梯度法或擬牛頓法更新參數場，不斷迭代演算以 達最佳未知數估計，此一般化逆向問題求解的流程如圖 2.1 所示。
圖 2.1 地下水一般化逆向問題求解流程圖 以下描述最佳化估計所有模式未知數的每個步驟：
步驟 1、給定初始猜值：一開始對於所有未知的參數、初始條件或邊界條件 給定其初始猜值，此猜值可依據現場初步的量測與實驗或者模式使用者本身的經 驗判斷，給定一個值與其合理範圍；此初始猜值只是一個初步估計值，此後會根 據最佳化演算法求解其最佳估計值。
步驟 2、正向模擬：利用控制方程式、初始條件與邊界條件，在給定的未知 數初始猜值下，正向模擬以求解地下水模擬問題（simulation problem），得到模 擬狀態變數分布場。
步驟 3、計算目標函數：計算全域目標函數，即模擬狀態變數與實際現場觀 測的誤差平方和。目標函數是一個具體的量化指標，藉由狀態變數場的誤差評估 未知參數與條件的估計誤差。最佳化的過程將不斷調整所有的未知數，使得模擬 與觀測越來越接近，目標函數達最小。
步驟 4、判斷是否收斂：判斷目標函數值是否小於一給定之容許值（given 1、給定初始猜值
2、正向模擬
3、計算目標函數
計算結束、獲得最佳化估計 是
4、是否 收斂？
否 5、伴隨問題模擬
7、更新猜值
6、計算梯度
tolerance），容許值由模式使用者本身決定，可以依據量測資料的多寡與精確度 給定。如果目標函數已小於容許值，表示模擬與觀測的誤差已降低至可接受範圍 內，參數與條件以達最佳化，逆向問題求解完成。
步驟 5、伴隨問題（adjoint problem）模擬：若目標函數大於容許值，表示 未知數估計偏離實際情形，導致模擬狀態變數與實際觀測仍有段差異，此時將這 些誤差放入伴隨狀態方程式的源匯項中，配合伴隨邊界條件與終點條件（final condition），求解伴隨問題，得伴隨狀態變數場。
步驟 6、計算梯度：利用伴隨狀態法推導所得之梯度積分式，此梯度式僅為 狀態變數與伴隨狀態變數的函數，將其依各個未知數所代表的時空區域進行積 分，即可同時計算在此估計未知數場下目標函數對應於所有未知數的梯度值。
步驟 7、更新猜值：將計算所得之梯度值以梯度法或擬牛頓法計算未知數的 修正方向，再輔以直線搜尋法（line search method）決定最佳的修正步幅，一次 更新所有未知的參數、初始條件與邊界條件，使其往真值場步步逼近。
重複步驟 2 至步驟 7，不斷迭代直到滿足收斂門檻，目標函數值小於容許值 為止，如此完成了地下水逆向問題的求解，得到最佳的參數、初始條件與邊界條 件估計。
上述的最佳化流程，已經將地下水逆向問題給一般化，此流程可適用於地下 水參數檢定以及初始條件與邊界條件的估計，最佳的未知數分布可以很有效率地 在幾次迭代（iteration）搜尋之後完成，使得檢定所得之地下水模擬模式最可以 適切的描述實際狀況。以下將就應用目標之不同，分別描述如何將本研究所提之 方法論應用於地下水參數檢定與完整的逆向問題求解上。
2.2 地下水參數檢定（Parameter Identification）
對污染場址特徵調查而言，井際分溶示蹤劑試驗以普遍使用於評估現場的污 染分布情形，其收集了富含訊息的水頭與動差觀測資料，應不只限定於估計污染 分布，而可往回一步檢定地下水模式，使用守恆示蹤劑的動差資料配合水頭觀
測，可以正確估計水利傳導係數。本節描述此一地下水參數檢定問題的求解過 程，包含地下水流與示蹤劑動差控制方程式、目標函數、伴隨問題、梯度積分式 與參數更新方法。本文中僅呈現重要結果，伴隨狀態法的理論與詳細推導過程置 於附錄 A 中。
2.2.1 地下水水流與示蹤劑動差控制方程式
關於地下水參數檢定問題，本研究以水平二維等向性拘限地下含水層
（horizontal two-dimensional isotropic and confined aquifer）為研究對象，探討其 水流與示蹤劑傳輸問題，所用到的方程式包括有：連續方程式、達西方程式、傳 輸方程式積分後之零階動差方程式與一階動差方程式，分述如下。其中穩態地下 水流之控制方程式為：
( ) Qsource Qsink 0
B B
q (2.1)
, (x,y) K H
q 0 (2.2)
方程式(2.1)為基本之連續方程式，代表了一個拘限含水層中內部水量的增減 與其邊界面上的流通量變化相等；方程式(2.2)為達西方程式，此式說明流通量的 大小與單位距離的水頭差成正比，其中乘上一個參數水力傳導係數 K，這是一個 地下含水層的水文地質參數，此參數越大表示地下水越容易在含水層中流量，反 之則流動不易。
而這兩個方程式有其附帶的邊界條件，其通式表示如下：
1
2
1
2
|
|
H f
f
q n (2.3)
其中
H =H(x,y)是壓力水頭（piezometric head），單位是長度[L]
( , ) ( , )
x y
q x y q x y
q i j 是達西流速（Darcy velocity），單位是速度[L/T]
K=K(x,y)即水力傳導係數場，也是速度單位[L/T]
B 是含水層厚度，此處假設含水層厚度均等，為一固定常數值，單位是長度[L]
Qsink=Qsink(x y, )是在抽水井位置才有的固定抽水量，為單位時間單位面積的抽 出水量，所以單位是長度除以時間[L/T]
Qsource= Qsource(x y, )是注入井處的注水量，單位是長度除以時間[L/T]
是整個水流的模擬區域，邊界位於與1 ，2 是一個給定水頭邊界；1 是一2 個給定流通量的邊界。
n = ( , ) x y
n 是邊界上的一個垂直於邊界的向外單位法向量。2 f1= f1(x,y)、 f2= f2(x,y)則是用於描述邊界條件的兩個函數。
在拘限含水層中，不考慮吸脫附（adsorption/desorption）的非穩態示蹤劑傳 輸方程式如下：
+ ( ) ( ) Qsink Qsource s
C C C C C
t B B
v D (2.4)
傳輸方程式(2.4)其附帶的初始條件與邊界條件如下，此處考量實際現場之示 蹤劑試驗，實驗前並無任何示蹤劑於含水層中，所以初始濃度為零。
1
2
3
4
0 on , 0
|
|
C t
C f
C f
D n
(2.5)
其中
C=C(x,y,t)是示蹤劑的濃度，其為時間與空間的函數，單位是質量除以體積[M/L3] Cs=Cs(t)是由注入井處所注入之示蹤劑濃度歷線，單位為質量除以體積[M/L3]
= ( , ) qx qy x y u v
v v i j i j是孔隙流速，等於達西流速除上孔隙率大小，單位
是長度除以時間[L/T]
f3= f3(x,y)、f4= f4(x,y)都是已知函數，用於描述邊界條件
是含水層的孔隙率，是一個比例所以沒有單位[ ]
2
2
( ) / ( ) /
( ) / ( ) /
xx xy T L T L T
yx yy L T T L T
D D V u V uv V
D D uv V V v V
D 是水力延散係數
（hydraulic dispersion coefficient），單位是長度二次方除以時間[L2/T]，此處假設 分子擴散（molecular diffusion）很小可以忽略。其中的速度V u2 、v2 是L 縱向延散度（longitudinal dispersivity）是長度單位[L]、 是橫向延散度（transverseT dispersivity）是長度單位[L]， 和L 皆假設其在研究區域中為一定值常數。T
將傳輸方程式以及邊界條件對時間域無窮積分（實際上僅需積分至試驗完成 時間即可），獲得零階動差方程式與其邊界條件。積分完成後，原來的非穩態傳 輸模擬模式變為一個穩態的動差模擬模式，其控制方程式與邊界條件如方程式 (2.6)、(2.7)所示。
0 0 0
( ) ( ) sink source s s
Q Q
m m m t C
B B
v D (2.6)
1
2
0 5 3
0 6 4
|
|
c
c
m f t f
m f t f
D n (2.7)
其中
0 0( , ) 0 ( , , )
m m x y
C x y t dt是零階動差，是將突破曲線對時間積分，所以單 位是濃度乘上時間[TM/L3]ts是示蹤劑注入時間[T]
tc是試驗完成時間[T]
f5= f5(x,y)、f6= f6(x,y)皆為已知函數，用於描述邊界條件
由動差的定義可以看出，零階動差代表了整個突破曲線下方所包含的面積大 小，即為累積通過的示蹤劑質量除以觀測位置的流通量；其積分範圍從時間零開 始至時間為無窮大，如此才能確定示蹤劑以被完全回收，突破曲線以包含完整的 場址訊息。實際現場觀測時，是取樣至濃度很小即可停止試驗，將時間離散的瞬 間觀測，利用數值積分，可估計得到零階動差的觀測。而濃度取樣頻率需依據突 破曲線之斜率變化進行取樣，在斜率變化大的濃度上升段以及濃度下降段，取樣 點要多；試驗出奇雨後其較為平緩的濃度變化取樣之時間間隔可以較長，如此才 可利用有限的量測資料代表完整的突破曲線，再由數值積分所得之動差才會正
確。
零階動差方程式是由傳輸方程式積分而來，所以其基本形式與傳輸方程式相 同，唯其時間項已經被積分掉了。等號右手邊的源匯項代表了注入示蹤劑的多 寡，此項會直接影響模擬所得之零階動差大小。零階動差是突破曲線與時間軸包 圍而成的面積，在一個運移主控（advection-dominant）的穩態地下水流場中（即 擴散作用相對而言較小），零階動差值在示蹤劑流經範圍（plume）內，等於注入 示蹤劑質量除以注入速率，也等於注入時間乘以注入濃度（m0 t Cs s），是一個 常數值；相反地，在示蹤劑沒有掃瞄過的地方，其零階動差值恆等於零，形成了 一個二元的分布型態。Cirpka and Kitanidis (2000)首先發現這樣一個二元分布行 為，這樣的分布只適合用於辨識示蹤劑的流徑，以及劃分參數可檢定的區域，卻 無法對於參數檢定提供有意義的資訊。為了解決零階動差資訊不足的問題，本研 究提出了一個新的示蹤劑釋放策略（release strategy），以巧妙的操作零階動差分 布場，使之產生一個空間連續變化的分布形式，如此一方面由於此分布形式與水 利傳導係數場分布相關，可增加其所包含的參數場資訊；另一方面連續的零階動 差分布場，使得零階動差可以被引入反向問題的求解流程中，其目標函數對應於 未知參數的梯度值可以被正確的估計。
一階動差方程式同樣由傳輸方程式推導而來，是將傳輸方程式與邊界條件乘 上時間一次方後對無窮時間域積分；而二階動差方程式則是乘上時間二次方後對 時間積分，三次、四次動差等依此類推。推導方法可以參考 Harvey and Gorelick (1995)的文章，裡頭有詳細的方法論、推導流程與模擬結果。本研究直接引用其 成果，由於動差方程式已經是穩態方程式了，所以不再帶有初始條件，一階動差 方程式與其對應之邊界方程式如下所示：
2
1 1 1 0
( ) ( )
2
sink source s s
Q Q t C
m m m m
B B
v D (2.8)
1
2
1 7
1 8
|
|
m f
m f
D n (2.9)
其中
1 1( , ) 0
m m x y
Ctdt是一階動差，為一個二維空間的狀態變數，單位是時間二 次方乘上濃度[T2M/L3]f7= f7(x,y)、f8= f8(x,y)為已知函數，用於描述邊界條件
一階動差是示蹤劑從注入井到觀測井的抵達時間（arrival time）乘上了抵達 時的濃度大小作為權重，所以再將一階動差除以零階動差（所有濃度值總和）進 行正規化（normalize），得到了平均抵達時間（mean arrival time），這是一個很重 要的資訊，他代表示蹤劑於含水層中傳輸的平均速度快慢，平均抵達時間越久則 示蹤劑流速越慢，示蹤劑流速又與地下水流速相近，所以也包含了場址的水利傳 導係數的資訊。實際現場觀測時，是將觀測所得之濃度變化歷線，乘上其對應之 觀測時間與離散之時間間隔，然後累積加總得到一階動差的觀測。
由一階動差方程式(2.8)可以看出，積分完成後，注入的示蹤劑以及零階動差 皆以源匯項形式出現在等號的右手邊，實際上所有高階的動差方程式皆包含了比 自身低一階的動差於源匯項中，所以在求解高階動差方程式之前，一定要依序求 解比其低階的動差。以一階動差為例，首先求解零階動差方程式，再將零階動差 場放入一階動差方程式的源匯項中，求解得一階動差分布場。上述四個方程式於 實際使用時，是以有限插分法進行離散後，再以 FORTRAN 語言編碼為電腦程式。
動差是濃度對時間積分後的資料，包含了整合性、壓縮的（lumped）資訊，
比被時間稀釋的（temporally diluted）濃度觀測資料更可以反應出水力傳導係數 場的變化，所以更適合被使用於參數檢定的問題。而且無論幾階的動差方程式皆 為穩態方程式，這些方程式的形式也都一樣，所以相較於非穩態的傳輸方程式而 言，其求解容易且精度提高。基於上述的資訊充足性與計算的簡單性，本研究中 選擇以穩態的零階動差與一階動差模擬模式取代非穩態的傳輸模式，以更正確的 描述示蹤劑的傳輸行為，並提供更有效的觀測資料以逆向推估水力傳導係數場。
此外，在邊界條件的設定上，本研究同時採用了給定通量（given flux）與給
定值（given value）這兩種邊界條件，如此在正向模擬時才能正確求解。若邊界 條件均為給定通量，則無狀態變數值的限制，其值高低變成相對值而非絕對值；
相反地，若全部都是給定值邊界條件，則場址的水力梯度確定但通量卻不知道，
則任意的參數場皆可滿足水位觀測，所以在數學上的處理，邊界條件需包含兩種 形式。但是此限制可以透過內部計算條件給定的方式進行處理，如在四周為給定 水頭的邊界條件下，再額外給定內部一點的抽水或注水量，狀態變數場即可得唯 一解；同樣地，若四周都是不透水的邊界條件，也只需給定一個內部點的水頭值，
全域的狀態變數場也可解得唯一解。但這樣的處理方法在應用伴隨狀態法時，將 面臨發散的問題，而變得不適用，此點將在 2.2.3 小節說明。
2.2.2 目標函數
求解地下水逆向問題的最佳化方法，是由一個初始猜測的參數場出發，在此 參數場分布下配合水流與動差模式模擬水頭、流速、零階動差與一階動差的空間 分布。參數檢定的目標就是透過不斷調整參數分布場，使得模擬結果符合觀測資 料，如此才可以確定檢定得到的參數場符合實際現場狀況，為最佳參數。
估計過程中，會設立一個具體的量化指標稱為目標函數，以判斷參數估計之 正確與否。此指標乃模式計算值與現場實際觀測值的差異，差異越小表示模擬狀 態變數場越正確，也就表示了估計的參數場越正確，參數檢定的過程也就是目標 函數最小化的過程。一般選擇以計算值與觀測值之誤差平方和作為目標函數，是 一個可微分的二次方程式，分為單一觀測點誤差的局部目標函數 f 與所有研究區 域內觀測點誤差總和的全域目標函數 E，計算式如下。
2 2 2 2 2
0 0 1 1
2 2 2 2 2
0 0 1 1
,
1[( ) ( ) ( ) ( ) ( ) ]
2
1 [( ) ( ) ( ) ( ) ( ) ]
2
com obs xcom xobs ycom yobs com obs com obs
com obs xcom xobs ycom yobs com obs com obs
x y
f H H q q q q m m m m
E f d
H H q q q q m m m m
(2.10)
其中Hcom、qxcom、qycom、m0com、m1com分別是模式模擬的水頭、x 方向流速、y
方向流速、零階動差與一階動差；而Hobs、qxobs、qyobs、m0obs、m1obs則是其相對 應的現場觀測資料。
由於此目標函數整合了所有可得的觀測資料，但觀測資料本身數量級不同
（order），如果直接使用在後續參數檢定時會發現檢定過程由數值最大的一階動 差所主控，其他觀測的資訊無法被有效融合，所以上述的觀測資料在用於計算目 標函數前需事先經過正規化（normalized）使其數值大小相近，將位於零到一之 間，如此使得各個狀態變數的收斂速度相近，所有的觀測資訊可以有效的利用。
本研究中，水頭觀測資料以初始水頭正規化之；流速觀測以最大流速觀測正規化 之；零階動差是以示蹤劑平均注入濃度與試驗時間正規化之；一階動差是以示蹤 劑平均注入濃度乘上試驗時間的平方正規化之，正規化後所有觀測資料已成為無 因次。
2.2.3 伴隨問題
目前為止，模擬模式已架構完成，可以透過數值模擬計算流場與示蹤劑傳輸 歷程；目標函數也已設定完成，可以評估模擬結果的好壞；這時只需要有一個可 以有效修正參數的演算法即可完成求解地下水反向問題。這個演算法必須可以評 估參數估計錯誤對於模擬結果的影響，以正確的計算目標函數對應於未知參數的 梯度值（the gradients of objective function with respect to unknown parameters），梯 度之定義為任意一個參數值的單位變化所造成的目標函數變化大小，此一梯度值 指示了參數該如何更新。
一般的逆向演算法耗費最多的計算時間在求取這些梯度值，在參數眾多的情 況下，將大幅降低演算法的求解效率。為了增進反演的效率，本研究選用 Sun and Yeh (1990)所提出之伴隨狀態法以加快求解梯度值的速度，提升最佳化方法的效 率。伴隨狀態法基本作法是將原本的模擬問題（simulation problem）取變分，經
