具多內部源二維多物種溶質傳輸解析解模式發展及現地應用

國立臺灣大學生物資源暨農學院生物環境系統工程學研究所 碩士論文

Department of Bioenvironmental Systems Engineering College of Bioresources and Agriculture

National Taiwan University Master Thesis

具多內部源二維多物種溶質傳輸解析解模式發展及現地應用 Development and field application of analytical model for two-

dimensional multispecies transport involving multiple internal sources

許時凱 Shih-Kai Hsu

指導教授：劉振宇 博士 、陳瑞昇 博士

Advisor: Chen-Wuing Liu, Ph.D.、Jui-Sheng Chen, Ph.D.

中華民國 105 年 7 月

July, 2016

誌謝

兩年的碩士生涯即將畫下尾聲，在這段時間學習到很多知識，也認識了很多 人，並受到許多人幫助，使本論文能夠順利完成，並感謝所有幫助過我的貴人。

首先特別感謝指導教授劉振宇老師，提供豐富的資源與研究環境，並細心指導與 協助修改論文。感謝共同指導教授陳瑞昇老師的指導與幫助，使我能夠順利地在 中央大學完成我的論文研究。

此外也特別感謝口試委員陳世楷老師，大學時期即受到老師很多幫助，並給 予日後研究生活許多方向與建議。感謝口試委員張誠信老師提供論文上寶貴的建 議與指正。首先感謝台大研究室的大家，感謝學姊高雨瑄博士在台大研究生活的 所有幫助，熱心提供我很多研究方向與資源。感謝范致豪老師每次 Meeting 上的 指導，感謝卓姐在行政事務上的協助，感謝學長蔡政諺博士的關心與幫助，感謝 庭祐在研究生涯的鼓勵與陪伴，感謝學長宗男、學弟振佑。

感謝中央大學污染傳輸研究室的各位，讓我在中央的研究生活感受到像家人 一樣的溫暖，感謝學長賴庚辛博士在解析解模式上的指導與相關建議，感謝學姊 慶芳、秉叡、懷恩、思燕、爽德及虛子，謝謝你們豐富我的研究生活。感謝北科 研究室的大家，使我在研究生活遇到挫折能夠更加堅強，感謝學長誠斌在論文上 的提攜，感謝狗任與致恆。

感謝大學好友智豪、重諺、寶賢及翊群，感謝高中好友峻毅，研究生活總是 互相給對方打氣。最後我要感謝我阿婆、阿公、爸爸與媽媽的幫助與體諒，感謝 哥哥時倫研究生涯的關心與幫助。一路走來受過太多人的幫助與鼓勵，要感謝的 人太多了，不及備載，在此一併致上謝意。

摘要

在現地污染整治上常見場址內有多處污染源長期釋出的情況，然而前人發展 之移流－延散方程式多為邊界源(boundary source)溶質傳輸解析解，在模擬現地場 址污染問題時容易受到限制且較難廣泛應用。而考慮內部源(inner source)之溶質傳 輸解析解模式可模擬流場任意位置污染源注入情形。且過去發展之單物種傳輸模 式，無法模擬序列降解下多物種的衰減與傳輸情形。本研究發展具多內部源二維多 物種溶質傳輸解析解模式，並考慮移流、延散、線性平衡吸附和一階衰減情況，以 描述溶質在二維孔隙介質中的傳輸行為，模式亦可應用在模擬放射性核種衰變鏈、

氨氮及含氯有機物之降解過程。解析解推導依續使用 Laplace 轉換、finite Fourier cosine 轉換及廣義型積分轉換消去時間及空間微分項，將偏微分方程式轉換為代數 方程式進行求解，再利用一系列逆轉換求得原域之解析解。本研究發展之解析解與 有限差分(Finite difference method)數值方法進行驗證，結果顯示兩者相當吻合，確 定發展模式之正確性。模式應用在 Dover 空軍基地第六區，模擬含氯有機溶劑污 染場址 PCE 及 TCE 污染團分別在在 10 年、20 年、30 年及 40 年情境之分布與擴 散情形，並與現地場址數據做比較，闡述模式的適用性。本研究可了解含氯有機物 序列降解上多物種污染物傳輸行為，且作為現地場址初步評估污染整治的基礎。

關鍵字：內部源、邊界源、污染傳輸、二維多物種、解析解、現地應用

Abstract

The long term releases of multiple sources are common situation in remediation of contaminated site. However, most of the advective-dispersive equation from literatures are analytical model for solute transport developed only for boundary source, which could be limited to simulate the transport of contaminant for filed-scale site. Considering the analytical model for solute transport of inner source could simulate any position for the contaminant sources injection. The analytical model for single-species solute transport could not simulate the decay and transport of multispecies in the condition of sequential decay chain. This study presents an multiple sources analytical model for two- dimensional multi-species advective-dispersive equation which considered linear equilibrium sorption and first order decay reaction. The model could also be applied in predicting the behavior of decaying contaminant such as radionuclide, nitrogen transformation and dissolved chlorinated solvent. Generalized analytical solutions are derived through the sequential application of the Laplace, finite Fourier cosine and generalized integral transform to remove the temporal and spatial derivatives in the coupled partial differential equation. The system of coupled partial differential equation is tranform to a set of linear algebraic equations, and the solutions in the original domain are then obtained through consecutive integral transform inversions. The developed analytical solution is tested by comparing their results against the numerical solutions using a finite difference scheme. Results show perfect agreements between the analytical and numerical solutions. The developed model was applied to analyze field-scale transport and biodegradation processes occurring at the Area-6 site in Dover Air Force Base, and predict the distribution of PCE and TCE plumes in different scenario (10 years, 20 years, 30 years and 40 years). The developed model in this study is a useful tool for

preliminary assessment to the remediation of the field site.

Keyword: Inner Source, Boundary Source, Contaminant Transport, Two-dimensional Multi-species, Analytical Solution, Field Application

目錄

謝誌……….……....i

摘要………...ii

Abstract……….………....iii

目錄………..…..v

表目錄……….….vii

圖目錄………..…...viii

符號表……….…...ix

第一章 前言……….……….……1

1-1 研究背景……….………...1

1-2 研究目的………..….……...2

第二章 文獻回顧……….3

2-1 多物種解析解模式發展……….3

2-2 具內部源解析解模式發展………5

2-3 模式現地應用……….5

第三章 二維多物種溶質傳輸解析解數學模式……….8

3-1 數學模式建立……….……….8

3-2 解析解推導………..13

第四章 模式收斂性測試與驗證...………...24

4-1 收斂性測試………..24

4-2 模式驗証………34

第五章 模式現地場址應用………...………41

5-1 現地場址應用……….………...41

5-1-1 背景介紹………41

5-1-3 模式現地場址應用………42

5-2 結果與討論………42

5-3 各模式之比較………44

第六章 結論與建議………55

6-1 結論………..55

6-2 建議………56

參考文獻……….57

附錄...……….62

表目錄

表 3-1 源/匯項函數與對應之 Heaviside 函數...10

表 4-1 數值收斂性測試所輸入之參數設定...27

表 4-2 四物種在 Pe=10 時，不同累加項數 l 的收斂情況...28

表 4-3 四物種在 Pe=10 時，不同累加項數 n 的收斂情況...29

表 4-4 四物種在 Pe=25 時，不同累加項數 l 的收斂情況...30

表 4-5 四物種在 Pe=25 時，不同累加項數 n 的收斂情況...31

表 4-6 四物種在 Pe=50 時，不同累加項數 l 的收斂情況...32

表 4-7 四物種在 Pe=50 時，不同累加項數 n 的收斂情況...33

表 5-1 Dover 空軍基地現地參數(From: Clement et al. 2000)...48

表 5-2 Dover site 污染源預估質量釋出率，ṁ (kg/year)...48

圖目錄

圖 3-1 二維有限域溶質傳輸概念模型...11

圖 3-2 具多內部源二維多物種溶質傳輸解析解推導流程...12

圖 4-1 內部源收斂性測試概念化模型...26

圖 4-2 內部源模式驗證之概念模型...36

圖 4-3 物種一解析解與數值解在不同時間下濃度分布情形...36

圖 4-4 物種二解析解與數值解在不同時間下濃度分布情形...37

圖 4-5 物種三解析解與數值解在不同時間下濃度分布情形...37

圖 4-6 物種四解析解與數值解在不同時間下濃度分布情形...38

圖 4-7 PCE 之解析解與數值解濃度分布情形...38

圖 4-8 TCE 之解析解與數值解濃度分布情形...39

圖 4-9 DCE 之解析解與數值解濃度分布情形...39

圖 4-10 VC 之解析解與數值解濃度分布情形...40

圖 5-1 Dover site 場址示意圖(From: Clement et al. 2000) ...46

圖 5-2 Dover site 地下水位等高線圖(單位=ft)(From: Clement et al. 2000)...46

圖 5-3 Dover site 概念化流場模型...47

圖 5-4 本研究發展模式模擬 t=10 years 的 PCE 濃度分布...49

圖 5-5 本研究發展模式模擬 t=20 years 的 PCE 濃度分布...49

圖 5-6 本研究發展模式模擬 t=30 years 的 PCE 濃度分布...49

圖 5-7 本研究發展模式模擬 t =40 years 的 PCE 濃度分布...50

圖 5-8 現地數據模擬 t =40 years 的 PCE 濃度分布...50

圖 5-9 現地數據與本研究發展模式在 40 years PCE 濃度分布為 200 μg/L 比較....50

圖 5-10 本研究發展模式模擬 t =10 years 的 TCE 濃度分布...51

圖 5-11 本研究發展模式模擬 t =20 years 的 TCE 濃度分布...51

圖 5-13 本研究發展模式模擬 t =40 years 的 TCE 濃度分布...52

圖 5-14 現地數據模擬 t =40 years 的 TCE 濃度分布...52

圖 5-15 現地數據與本研究發展模式在 40 years TCE 濃度分布為 1000 μg/L 比較..52

圖 5-16 本研究發展之模式與現地數據之 PCE 濃度分佈比較...53

圖 5-16 RT3D 與現地數據之 PCE 濃度分佈比較...53

圖 5-16 MIKSS 與現地數據之 PCE 濃度分佈比較...53

圖 5-16 本研究發展之模式與現地數據之 TCE 濃度分佈比較...54

圖 5-16 RT3D 與現地數據之 TCE 濃度分佈比較...54

圖 5-16 MIKSS 與現地數據之 TCE 濃度分佈比較...54

符號說明

符號 說明 因次

Ci(x,y,t) 物種 i 溶質濃度 [ML⁻³]

Gi(X,Y,s) 物種 i Laplace 轉換後的濃度 [ML⁻³] Hi(X,n,s) 物種 i finite Fourier cosine 轉換後的濃度 [ML⁻³] Ui(X,n,s) 物種 i 變數變換後的濃度 [ML⁻³] Z_i(ξ_l,n,s) 物種 i 廣義型積分轉換後的濃度 [ML⁻³]

Bi 物種 i 連續注入之溶質濃度 [ML⁻³]

DL 縱向延散係數 [L²T⁻¹]

D_T 側向延散係數 [L²T⁻¹]

L 介質長度 [L]

W 介質寬度 [L]

v 地下水流速度 [LT⁻¹]

Pe Peclet number [−]

x 空間座標 [L]

y 空間座標 [L]

t

Ri 物種 i 的遲滯因子 [−]

ki 物種 i 的一階衰減常數 [T⁻¹]

Si 物種 i 源/匯項濃度強度參數 [−]

px,i(x) 物種 i 源/匯項空間變化函數 Py,i(x) 物種 i 源/匯項空間變化函數 qi,i(t) 物種 i 源/匯項時間變化函數

ξl 特徵函數

K(ξl,X) 正規化特徵函數

N(ξ_l) 範數

Qi(s) 物種 i Lapace 轉換後源/匯項時間變化函數

P_y,i(n) 物種 i finite Fourier cosine 轉換後源/匯項空間變化函數 Px,i(ξl) 物種 i 廣義型積分轉換後源/匯項空間變化函數

第一章 前言

1-1 研究背景

近年來台灣工業迅速發展，工業製程產生的廢棄物若未經妥善處理或非法棄 置容易造成環境污染，其中污染物若處置不當或隨雨水滲漏至地表下，將會對土 壤與地下水系統造成嚴重污染，並進一步威脅生態系統或人民生命安全，故污染 場址的管制與整治工作相當重要。我國於民國八十九年公告施行「土壤及地下水 污染整治法」針對有土壤、地下水污染之虞地區優先檢測，進行相關水文地質參 數調查，並依照調查結果訂定污染整治計畫。

為了預測污染物在地下水系統中的傳輸與宿命，我們常用污染物傳輸模式進 行模擬，而污染傳輸模式的建立通常採用移流－延散方程式，其中移流與延散傳 輸為污染物在地下水系統中主要機制，移流主要描述污染物隨地下水流移動，而 溶質在地下水系統中會受到孔隙異質性影響導致水流速度與水流路徑差異，故加 入延散觀念來描述溶質在地下水系統中的分散現象。

在處理現地污染整治問題時，常見場址內有多處多物種污染源釋出情況，而前 人發展之溶質傳輸解析解多為邊界源(boundary source)注入，在模擬現地場址污染 問題時容易受到限制且較難廣泛應用，考慮內部源(inner source)注入之溶質傳輸模 式可模擬流場任意位置污染源注入情形。且過去發展之單物種傳輸模式，無法模擬 序列降解下多物種的衰減與傳輸情形。為了模擬及預測現地場址多物種污染團的 傳輸與分布情形，以作為土壤與地下水污染整治之基礎，本研究發展多內部源二維 多物種溶質傳輸解析解。

本研究選擇 Clement et al. (2000)之 Dover 空軍基地污染場址作為模式發展與 應用之案例，場址由於過去各種軍事作業產生的廢棄物或污染物滲漏，造成土壤及 地下水遭受含氯有機溶劑嚴重污染。本研究應用發展的模式模擬污染場址內四氯 乙烯(PCE)及三氯乙烯(TCE)污染團之傳輸與分布情形，並與現地數據作比較，闡明

1-2 研究目的

現地污染場址通常為場址內有多處污染源注入情形，前人研究大多著重在發 展邊界源溶質傳輸模式，在處理現地場址污染問題時容易受到限制且較難廣泛應 用，考慮內部源注入之溶質傳輸模式，可模擬流場任意位置污染源注入情形。本研 究考慮序列降解下多物種溶質傳輸反應，發展多內部源二維多物種溶質傳輸解析 解模式，依序利用 Laplace 轉換、finite Fourier cosine 轉換及廣義型積分轉換移除 時間及空間之偏微分項，再進行一系列逆轉換及代數運算求得二維多物種解析解 通解，模式可模擬任意數目污染源注入情況，及序列降解鏈下各種污染物衰減過程，

如放射性核種衰變鏈、氨氮轉換及含氯有機物自然衰減過程，進一步了解多物種污 染團的衰減與傳輸情形。推導之解析解模式進行收斂性測試，並與數值模式進行驗 證，以確認解析解的正確性。為了測試發展之模式於現地場址的適用性，本研究利 用過去文獻含氯有機污染場址之相關參數與現地數據進行模擬，預測 PCE 及 TCE 污染團濃度分布情形，並與現地濃度分布數據作比較，希望所發展的解析解能夠應 用在現地場址的整治與調查上。

第二章 文獻回顧

為了評估孔隙介質中污染物的傳輸與宿命，各種溶質傳輸解析解模式和數值 模式相繼發展，而解析解模式可以了解土壤與地下水中污染物之行為與溶質傳輸 機制外，亦可作為數值模式評估的有效工具(Pérez, 2010 ; 田傑仁, 2013)，此外，

土壤與地下水污染整治的後續評估工作亦顯重要，故溶質傳輸模式在現地應用的 模擬與預測工作上扮演一個重要的角色。本研究旨在發展一套有效的多內部源二 維多物種溶質傳輸模式，並將模式應用在分析現地場址的整治評估上。

2-1 多物種解析解模式發展

過去數十年，溶質傳輸解析解問題已被廣泛研究與討論，Van Genuchten et al.

(1982)發展一維移流－延散方程式，其方程式描述一維空間座標下單物種溶質傳輸 問題，並考慮線性平衡吸附與一階降解反應，推導不同控制方程式與邊界條件下之 解析解。

Batu (1989)利用 Laplace 轉換與 Fourier analysis 推導二維單物種解析解，其解 考慮縱向與側向延散、放射性衰減反應與線性平衡吸附，並與有限元素法之數值解 比較驗證，結果十分吻合。

上述研究發展不同溶質傳輸問題的解析解，但僅限於模擬單物種溶質傳輸現 象。相對於單物種溶質傳輸模式，在討論放射性核種衰變鏈、氨氮與含氯有機物 之降解問題時，多物種溶質傳輸模式扮演極重要的角色，但其解析解發展較單物 種溶質傳輸問題來的困難且相對稀少。

Cho (1971)利用 Laplace transform 推導一維無窮遠下多物種溶質傳輸問題之解 析解，其解考慮依序一階降解反應下三物種傳輸模式，模擬地下水系統中氨氮的傳 輸與宿命。Van Genuchten (1985)利用 Laplace transform 推導一維無窮遠下四物種降 解反應，並模擬放射性核種衰變練與氨氮轉化過程。Lunn et al. (1996)利用有限 Fourier transform 推導 Cho (1971)之解析解。Chen et al. (2011)提出利用 Laplace 轉

換和廣義型積分轉換依序求解時間域和空間域偏微分方程式之新方法，發展有限 域下四物種解析解。這些利用 Laplace transform、finite Fourier transform 或廣義型 積分轉換推導多物種解析解方法大多受限於求解一維系統。

Sun and Clement (1999)與 Sun et al. (1999)發展通用方法求解多維多物種依序 降解反應，其方法可求解任意數目物種之解析解，但轉換方法只能適用於所有物種 遲滯係數相同時。Bauer et al. (2001)利用 Laplace transform 提出求解不同遲滯係數 下多物種移流－延散方程式方法，此方法利用遞迴公式將單物種基礎解擴展至更 複雜的多物種解，並可應用在二維及三維系統。

但上述方法需對各別物種依序求解其解析解，在處理不同數目物種問題時需 要重新推導解析解，在應用上缺乏適用性。Clement (2001)發展任意數目物種解析 解通解，此方法可應用在一階降解反應下一維、二維及三維系統，在求解多物種 溶質傳輸問題上有更大的適用性，但其在計算上受限於需假設所有物種遲滯係數 相同。

Quezada et al. (2004)延伸 Clement (2001)所提出的方法，利用 Laplace 轉換及線 性轉換法求解偏微分方程式，推導不同遲滯係數下多維多物種解析解通解，但其考 慮出流邊界為無窮遠。

然而上述研究多假設出流邊界為無窮遠，此假設較不符合實際狀況，Parlange and Starr (1984)、Parlange et al. (1985)、Parlange et al., (1992)提出在延散作用主導 的傳輸情形下，出流邊界條件為無窮遠情況對溶質傳輸有極大影響。Pérez Guerrero et al., (2009)亦提出污染傳輸模式在進行的數值模擬時，通常考慮為有限邊界，Chen et al., (2011)也發現在出流邊界為無限遠時會低估污染物濃度，因此本研究考慮出 流邊界條件為有限長度情況。

Chen et al. (2016)為了預測多物種污染團的傳輸與分布情形，延伸 Chen et al.

(2011)的研究，並依序利用 Laplace 轉換、finite Fourier cosine 轉換及廣義型積分轉

擬放射性核種衰變鏈與含氯有機污染團的傳輸。

2-2 具內部源解析解模式發展

Batu (1993)利用 Fourier analysis 和 Laplace 轉換求得均勻流場下多邊界源注入 模式，邊界源考慮隨時間及空間變換之函數，模式可模擬隨時間及位置變換之污染 源注入，但其考慮污染源為邊界源注入，無法模擬場址內任意位置污染源釋出情形。

V. Srinivasan et al. (2007)利用 Laplace 轉換和線性轉換方法導出一維多物種 closed-form 解析解，並各別推導隨空間變化之 Bateman-type 初始條件下邊界條件 分別為 Dirichlet 和 Cauchy 之解析解。然而其控制方程式不考慮源/匯項，故無法 模擬內部源注入情形。

上述研究多著重在發展邊界源溶質傳輸解析解，然而邊界源注入模式無法模 擬現地場址任意位置之污染源注入情形。考慮內部源注入之溶質傳輸模式，可模擬 現地多處污染源釋出情況，對現地場址的整治評估工作有很大的幫助。

Park et al. (2001)考慮有限厚度含水層，利用 Green’s function 方法求得具內部 源之三維溶質傳輸解析解，並分別考慮點、線及面注入情形，惟其發展之單物種模 式，無法模擬現地場址多物種溶質衰減情形。

Pérez Guerrero et al., (2009)考慮出流邊界為有限長度，利用廣義型積分轉換求 得具源/匯項之三維解析解，其源/匯項可模擬隨時間及空間變化的內部源，分別考 慮穩態及暫態下的溶質傳輸機制，但所發展的之單物種解析解模式亦無法模擬序 列降解下的傳輸行為。

上述研究在發展具內部源注入之單物種溶質傳輸模式，無法模擬多物種溶質 序列降解下的傳輸機制，在使用上有許多限制，本研究欲發展具多內部源之二維多 物種溶質傳輸模式，模擬現地場址多物種污染團衰減與傳輸行為。

2-3 模式現地應用

工業發展產生的廢棄物若處置不當容易造成土壤與地下水污染，而地表下污

染物的傳輸與宿命會受到孔隙介質中物理、化學及生物反應間交互作用影響，為了 模擬污染整治情形並建立污染場址相關風險評估，現地場址的模擬與相關模式發 展與應用非常重要。

Rifai et al. (1987)發展一個二維生物整治數值模式 BIOPLUME-II 模擬厭氧生 物降解反應下碳氫污染物的傳輸與宿命，且 Rifai et al. (1988); Chiang et al. (1989) 利用 BIOPLUME-II 模擬各種碳氫燃料場址的厭氧生物整治方法及其他污染物降解 機制。

Essaid et al. (1995)發展二維多物種反應傳輸模式 BIOMOC 模擬原油洩漏場址 自然生物衰減過程。BIOMOC 可模擬好氧反應下揮發性與非揮發性有機溶劑衰減、

錳還原、鐵還原及甲烷生成反應，並應用模式研究明尼蘇達州 Bemidji 原油洩漏場 址現地傳輸及衰減過程。

Clement (1997)與 Clement et al. (1998)發展多物種反應傳輸模式 RT3D 模擬地 下水系統中生物整治情境。不同於其他的生物整治模式通常只能解特定種類污染 傳輸問題，RT3D 為通用反應傳輸模式可模擬任意種類與任意數目在生物化學反應 下溶質傳輸方程式。例如 Lu et al. (1999)利用 RT3D 模擬好氧及連續厭氧環境下碳 氫化合物污染團降解及傳輸之反應方程式。Sun et al. (1998)利用 RT3D 模擬碳氫化 合物在細菌生長下降解的多物種傳輸模式。Clement et al. (1998)利用 RT3D 發展四 物種依序降解模式，並模擬含氯有機溶劑如 PCE 及 TCE 在厭氧反應下的自然衰減 過程。

針對現地場址污染整治，美國環保署亦推廣一些解析解模式，模擬現地污染物 自染衰減過程，例如 BIOSCREEN (Newell et al., 1996)及 BIOCHLOR (Aziz et al., 2000)，其中 BIOSCREEN 可模擬因石油洩漏而造成碳氫化合物污染場址中自然衰 減的過程，而 BIOCHLOR 可模擬受含氯有機物污染場址之自然衰減現象。Clement et al. (2002)利用 BIOCHLOR 模擬含氯有機污染物降解過程，並評估監測式自然衰

上述二種解析模式 BIOSCREEN 及 BIOCHLOR 均假設入流邊界為第一類邊界 條件。由於入流邊界為第一類邊界條件，相較於第三類邊界條件通常會高估濃度 (Parker and van Genuchten, 1984; Parlange et al., 1985)，故本研究採用第三類邊界條 件作為入流邊界發展解析解。

Clement et al. (2000)為了分析含氯有機污染場址自然降解過程，發展一個多維 多物種反應傳輸模式，該模式可模擬多物種同時進行的衰減過程，包含好氧及厭氧 生物降解過程。模式應用在分析德拉瓦州 Dover 空軍基地現地污染傳輸和生物降 解過程，並利用現地資料校正，校正後之模式不僅可驗證地下水流況，亦可重建 PCE、TCE、DCE、VC 及含氯污染團觀測分布。

Atteia et al. (2008)分別利用發展之半解析解模式 MIKSS 與 RT3D 模擬 Dover 空軍基地污染團濃度分布與相關降解反應，結果顯示兩模式的濃度分布圖形相當 吻合，並利用現地 PCE 及 TCE 同位素資料驗證模式。

本研究延伸 Chen et al. (2016)之方法發展具多內部源二維多物種溶質傳輸解析 解模式，為了模擬真實場址內各物種內部污染源情況，本研究改良並加入源/匯項 情形(凃佑霖, 2015)，考慮不同遲滯係數、線性吸附、一階降解反應。模式配套文獻 之含氯有機污染場址現地資料(Clement et al. 2000)，預測 PCE 及 TCE 污染團在各 種情境的濃度分布，並與現地數據比較，以說明模式的適用性。

第三章 二維多物種溶質傳輸解析解數學模式

3-1 數學模式建立

溶質在土壤與地下水系統中的濃度與分布，會受到溶質種類、地表下物理、化 學特性及生物種類間的交互作用影響。一般而言，會造成溶質分散的機制由延散和 擴散控制，但在地下水問題中由於擴散係數遠小於延散係數 （擴散係數介於 1~2×10⁻⁹m²/sec遠小於延散係數），因此在計算時分子擴散常被忽略。故討論污染 傳輸問題時須考慮以下機制：

(1)移流：溶質在孔隙介質中隨水流移動地行為，其中孔隙水流速度遵循達西定律 (Darcy’s law)。

(2)延散：由於孔隙介質在微觀尺度下的異質性，造成溶質在孔隙水流中的流徑及 流速不同，因此採用延散的觀念來描述達西定律的不足，並假設延散作用遵守 Fick’s law 的型式。

(3)吸附：溶質和孔隙介質表面結合的現象稱之為吸附。

(4)生物降解：土壤與地下水系統中的微生物會對特定污染物進行自然分解，將污 染物傳化成毒害較小的物質。

根據上述機制及假設，發展具內部源二維多物種溶質傳輸解析解。假設溶質傳 輸在二維等向、均質流場，且地下水流為平行 x 軸之穩態均勻(steady uniform)流場，

建立一個序列降解之聯立控制方程組，方程式加入源/匯項以模擬流場內部源注入 情形，並設定初始與邊界條件。其二維有限域溶質傳輸概念模型見圖 3-1。

為求出聯立控制方程組之解析解，首先將控制方程式及初始與邊界條件無因 次化，接著依序進行 Laplace 轉換、有限 Fourier cosine 轉換及廣義型積分轉換移除 時間及空間之偏微分項，最後透過一系列逆轉換求得原域之解析解，二維多物種溶 質傳輸解析解流程見圖 3-2，考慮線性平衡吸附及一階降解反應之序列降解二維控 制方程式及初始與邊界條件如下：

控制方程式：

t t y x R C

t q y p x p S

t y x C R k y

t y x D C

x t y x v C x

t y x D C

NS

m

m my

mx m

T L



 



 

 



 







) , , ) (

( ) ( ) (

) , , ) (

, , ( )

, , ( )

, , (

1 1 1

1 1 , 1 , 1

1 1 2 1

2 1 1

2 2 1

(3-1)

N i

t t y x R C

t q y p x p S t

y x C R k

t y x C R k y

t y x D C

x t y x v C x

t y x D C

i i NS

m

im mi m y

i m x i i

i i

i i i i

i T L i

...

2

) , , ) (

( ) ( ) ( )

, , (

) , , ) (

, , ( )

, , ( )

, , (

1

, ,

1 1 1

2 2 2

2





 





 

 



 











 (3-2)

此處

C

( x , y , t )

k

R



im mi m y

i m x

i p x p y q t S

1

,

, ( ) ( ) ( )為控制方程式之源/匯項[ML^3]，其中 N S 為內部源數

目；S_i^m為內部源 m 物種 i 濃度強度參數[]，p^m_x,i(x)與p^m_y,i(y)代表內部源 m 物種

i 隨空間變化的函數，q_i^m(t)代表內部源 m 物種 i 隨時間變化的函數。

初始條件：

假設實驗進行前，孔隙介質中無任何溶質存在，則初始狀態可表示為 0

) 0 , ,

(x y t  

C_i (3-3) 邊界條件：

入流處考慮為第三類邊界條件以滿足質量守恆，並假設 y 方向注入邊界為 y₁到

y2，其中B_i為物種 i 之注入濃度，H為 Heaviside 函數用以表示空間段注入情形。

則入流處之邊界條件可表示為

   

 

N 1 i

y y H y y H vB t y x x vC

t y x

D_L Cⁱ _{i i}

...

) , , 0 ) (

, , 0 (

2 1













 



 

(3-4)

考慮x方向出流處邊界濃度梯度變化為零，表示為 ) 0

, ,

( 







x t y L x C_i

i1...N (3-5) 考慮 y 方向上、下邊界濃度梯度變化為零，表示為

) 0 , 0 ,

( 







y t y x C_i

i1...N (3-6)

) 0 , ,

( 







y t W y x C_i

i1...N (3-7) 控制方程式(3-1)及(3-2)中，加入源/匯項S_i^mp^m_x,i(x)p^m_y,i(y)q_i^m(t)以模擬現地應 用情況，場址某空間座標上污染源隨時間持續或週期排放的現象，其中p^m_x,i(x)及

) , (y

p^m_yi 代 表 隨 空 間 方 向 x 及 y 變 化 之 函 數 ， 分 別 以

   

p H x x

x x

H  ₁   及

  



p H y y

y y

H  ₁   表示 x 方向上x^m_p1xx^m_p2及 y 方向上y^m_p1 y y^m_p2區域段所圍 成之矩形區域有源項之濃度釋放，q_i^m(t)表示隨時間 t 變化之函數，以H t



 

表 3-1 源/匯項函數與對應之 Heaviside 函數

源/匯項函數 對應之 Heaviside 函數

)

, (x

p_x^m_i

   

x

H  ₁  

) , (y

p^m_yi



 



p H y y

y y

H  ₁  

) (t

q_i^m H

 

 

圖 3-1 二維有限域溶質傳輸概念模型

圖 3-2 具多內部源二維多物種溶質傳輸解析解推導流程 STEP 1.設定控制方程式與初始、邊界條件

STEP 2.將控制方程式與初始、邊界條件無因次化

STEP 3.進行 Laplace 轉換，移除時間之微分項

STEP 4.進行 finite Fourier cosine 轉換，移除空間方向 y 之微分項

STEP 5.控制方程式變數變換與邊界條件齊次化

STEP 6.進行廣義型積分轉換，移除空間方向 x 之微分項

STEP 7.進行一系列逆轉換求得原域之解析解

3-2 解析解推導

STEP 1. 設定控制方程式與初始、邊界條件

本研究為發展二維多物種解析解模式，考慮控制方程式、初始及邊界條件如 下：

t t y x R C

t q y p x p S

t y x C R k y

t y x D C

x t y x v C x

t y x D C

NS

m

m my

xm m

T L



 



 

 



 







) , , ) (

( ) ( ) (

) , , ) (

, , ( )

, , ( )

, , (

1 1 1

1 1 , 1 , 1

1 1 2 1

2 1 1

2 2 1

(3-1)

N i

t t y x R C

t q y p x p S t

y x C R k

t y x C R k y

t y x D C

x t y x v C x

t y x D C

i i NS

m

im mi m y

i m x i i

i i

i i i i

i T L i

...

2

) , , ) (

( ) ( ) ( )

, , (

) , , ) (

, , ( )

, , ( )

, , (

1

, ,

1 1 1

2 2 2

2





 





 

 



 











 (3-2)

0 ) 0 , ,

( x y t  

C

   

 

N 1 i

y y H y y H vB t y x x vC

t y x

D_L Cⁱ _{i i}

...

) , , 0 ) (

, , 0 (

2 1













 



 

(3-4)

) 0 , , 0

( 







x t y x C_i

i1...N (3-5) ) 0

, 0 ,

( 







y t y x C_i

i1...N (3-6) ) 0

, ,

( 







y t W y x C_i

i1...N (3-7)

STEP 2. 將控制方程式與初始、邊界條件無因次化

發展數學模式時，首先將控制方程式與初始、邊界條件無因次化，讓所求出 的解能應用於不同的尺度大小。假設無因次參數

L X  x ，

W Y  y ，

W Y₁ y¹，

W Y₂  y² ，

L

T  vt代入方程式(3-1)-(3-7)進行無因次化，則方程式可改寫成

T T Y X R C

T q Y p X p

T Y X C Y

T Y X C Pe X

T Y X C X

T Y X C Pe

NS

m

m my

xm m

T L



 





 

 











) , , ) (

( ) ( ) (

) , , ) (

, , ( )

, , ( )

, , ( 1

1 1 1

1 1 , 1

, 1

1 2 1

2 1 1 2

2 2 1

 

(3 -8)

N i

T T Y X R C

T q Y p X p T

Y X C

T Y X C Y

T Y X C Pe X

T Y X C X

T Y X C Pe

i i NS

m

m i m

i y m

i x m i i

i

i i i

T i

i L

...

2

) , , ) (

( ) ( ) ( )

, , (

) , , ) (

, , ( )

, , ( )

, , ( 1

1

, ,

1 1

2 2 2 2

2





 







 

 













 

 

(3-9)

初始條件：

0 ) 0 , ,

(X Y T  

C_i (3-10) 邊界條件：

   





) , , 0 ) (

, , 0 1 (

Y Y H Y Y H B T Y X X C

T Y X C

Pe ^{i i}

i L











 



  i1...N (3-11)

) 0 , , 1

( 







X T Y X C_i

i1...N (3-12) ) 0

, 0 ,

( 







Y T Y X C_i

i1...N (3-13) ) 0

, 1 ,

( 







Y T Y X C_i

i1...N (3-14)

此處

L

L D

Pe  vL ，

T

T D

Pe  vL，

v L R k_{i i}

i 

 ，

W

 L

 ，

v S_i^m L

im 



針對(3-8)-(3-14)之控制方程式與初始及邊界條件，本研究連續使用 Laplace 轉換、有限 Fourier cosine 轉換及廣義型積分轉換將偏微分方程式轉化成代數方程 式求解，推導過程如下所示：

STEP 3. 進行 Laplace 轉換，移除時間之微分項

利用 Laplace 轉換移除方程式中的時間 t 之微分項，Laplace 轉換的定義如下：

dT T Y X C e s

Y X

G_i^{( , , )}





 0 ( )

)

(s e q T dT

Q_i^{m sT} _i^m (3-16) 此處G_i(X,Y,s)與Q_i(s)為 Laplace 轉換後之函數，s 為 t 之 Laplace 轉換參數。而 Laplace 逆轉換公式如下：

ds s Y X G i e

T Y X

C ⁱ _i

i

i sT ( , , )

2 ) 1 , ,

( 



 (3-17) 對方程式(3-8)-(3-14)進行 Laplace 轉換後可得

 











 





 





NS

m

m my

xm m

T L

s Q Y p X p

s Y X G s

X R s Y X G Y

s Y X G X Pe

s Y X G Pe

1

1 1 , 1

, 1

1 1 1 1

2 2 1 2 2

2 1

) ( ) ( ) (

) , , ) (

, , ( )

, , ( )

, , (

1  

(3-18)

 

N i

s Q Y p X p s

Y X G

s Y X G s

X R s Y X G Y

s Y X G X Pe

s Y X G Pe

NS

m

im myi

mi m x i i

i

i i i i

i T i

L

...

2

) ( ) ( ) ( )

, , (

) , , ) (

, , ( )

, , ( )

, , ( 1

1

, ,

1 1

2 2 2 2

2













 





 









 

 

(3-19)

   





) , , 0 ) (

, , 0 1 (

Y Y H Y Y s H s B Y X X G

s Y X G Pe

i i

i L











 



  (3-20)

) 0 , , 0

( 







X s Y X G_i

(3-21)

) 0 , 0 ,

( 







Y s Y X G_i

(3-22)

) 0 , 1 ,

( 







Y s Y X G_i

i1...N (3-23) N

i1...

N i1...

STEP 4. 進行 finite Fourier cosine 轉換，移除空間方向 y 之微分項

執行有限 Fourier cosine 轉換移除方程式中的 y 方向二階微分項。有限 Fourier cosine 轉換的定義如下：

dY Y n s

Y X G s

n X

H_i( , , )^



Y n Y

p n

P_y^m_i( ) ^{1 m}_yi( )cos( )

0 ,

, ^



) cos(

) , , ( 2 ) , 0 , ( ) , , (

1

Y n s n X H s

n X H s Y X G

n i i

i









 (3-26)

方程式(3-18)-(3-23)經有限 Fourier cosine 轉換後可得



















  





NS

m

m ym

xm m

T L

s Q n P X p

s n X Pe H

s n dX R

s n X dH dX

s n X H d Pe

1

1 1 , 1

, 1

1 2 2 2 1 1 1

2 2 1

) ( ) ( ) (

) , , ) (

, , ( )

, , (

1   

(3-27)

N i

s Q n P X p s

n X H

s n X Pe H

s n dX R

s n X dH dX

s n X H d Pe

NS

m

i i y i x i i

i

i T i

i i i

L

...

2

) ( ) ( ) ( )

, , (

) , , ) (

, , ( )

, , ( 1

1

, ,

1 1

2 2 2 2

2





















  









 



 

(3-28)

) ( )

, , 0 ) (

, , 0 1 (

s n s B n X dX H

s n X dH Pe

i i i

L







 

 i1...N (3-29)

) 0 , , 1

(  

dX s n X dH_i

i1...N (3-30)

其中

   

...

3 , 2 , 1 sin 0

) sin

( _{2 1}

1 2









 



 n

n n

Y n Y

n

Y Y n





 為 Heaviside 函數 H





 

有限 Fourier 轉換後之函數。

STEP 5. 控制方程式變數變換與邊界條件齊次化

進行廣義型積分轉換前，Pé rez Guerrero et al. (2009)建立先將邊界條件齊性化，

再利用變數變換消除 x 方向空間一階微分項之方法，本研究根據 Chen et al. (2011) 之方法將邊界條件齊性化及變數變換之函數定義為

) , , ( )

( )

, ,

( n e ² U X n s

s s B n X

H _i

e X P i i

L





 (3-31)

此處U_i(X,n,s)為齊性化及變數變換後之函數。

方程式(3-27)-(3-30)經由齊性化及變數變換後可得











 











  













   



NS

m

Pe X m

ym xm

m

Pe X i

T

L T

L

L L

e s Q n P X p

e s n

B Pe

s n R

s n X Pe U

Pe s n

R dX

s n X U d Pe

1

1 2 1 , 1

, 1

2 2

2 2 1 1

1 2

2 2 1 2 1

2 1

) ( ) ( ) (

) (

) , , 4 (

) , , ( 1



 



 

(3-32)

N i

e s Q n P X p

s n X U e

s n e B

s n B Pe

s n R

s n X Pe U

Pe s n

R dX

s n X U d Pe

NS

m

Pe X im

mi m y

i m x i

i Pe X

i i Pe X

i T i

i

L i T

i i i

L

L

L L

...

2

) ( ) ( ) (

) , , ( )

( )

(

) , , 4 (

) , , ( 1

1

, 2 ,

2 1 1 1

2 2

2 2

2 2 2 2

2





























 











  













   









 



 



 



 

(3-33)

0 ) , , 0 2 (

) , , 0

(    

 Pe U X n s

dX s n X U d

L i

i (3-34)

0 ) , , 1 2 (

) , , 1

(    

s n X Pe U

dX s n X U d

i L

i (3-35)

STEP 6. 進行廣義型積分轉換，移除空間方向 x 之微分項

本研究利用廣義型積分轉換移除 x 方向二階微分項(Cotta, R.M., 1993; Pérez Guerrero et al., 2009)，其求解過程包含以下步驟：

(1) 選 擇 合 適 的 特 徵 值 問 題 並 根 據 問 題 求 解 特 徵 值 (eigenvalue) 、 特 徵 方 程 式 (eigenfunctions)、範數(norm)及正交化特性(orthogonalization property)。

(2)發展廣義型積分轉換及廣義型積分逆轉換。

(3)將微分方程式轉化成代數方程式並求解。

(4)利用廣義型積分逆轉換求得未知的函數。

根據方程式(3-33)及邊界條件(3-34)、(3-35)定義特徵值問題，並求得特徵值_l、 正規化特徵函數K(_l,X)及範數N(_l)，其定義如下（詳細推導過程見附錄）：

) cos(

) 2 sin(

) ,

( Pe X X

X

K _l  ^L _l _l _l (3-36)

2 ) 4

(

2 2

l L L

l

Pe Pe N

  ^{ } (3-37)

廣義型積分轉換的定義如下：

dX s n X U X K s

n

Z_i(_l, , )^



e X K

Pe X l

l





 ¹

0

) 2

, ( )

( 

2 2

4 ^l

L l L

Pe Pe







 (3-39)

dX e

X p X K

P ^X

Pe i

x l

l i x

L



 ¹

0

, 2

, ( ) ( , ) ( ) (3-40) 而廣義型積分逆轉換之公式如下：







1 ( )

) , , ( ) , ) (

, , (

l l

l i i l

N

s n Z X s K

n X

U 



 (3-41)