曹紅勝 張琳

（1. 濱州市政務(wù)服務(wù)中心，山東濱州 256600；2. 濱州市引黃灌溉服務(wù)中心，山東濱州 256600）

1 引言

隨著我國(guó)城鎮(zhèn)化的不斷推進(jìn)，居住在城市周邊地區(qū)的人口數(shù)量越來越多［1-2］，如何有效保護(hù)水資源以及避免河流污染損害人類和生態(tài)健康成為一大挑戰(zhàn)［3-5］。目前，固體廢物和水污染是城市河流的主要污染源，并且人口的快速增長(zhǎng)與基礎(chǔ)設(shè)施建設(shè)之間的不平衡，也使得河流中存在大量的病原體和有機(jī)物質(zhì)。此外，大量文獻(xiàn)研究表明［6-8］，盡管城市地區(qū)的地表不易透水，但污染物仍會(huì)滲入地下污染地下水源，再加上過度抽水，各含水層間的相互作用增強(qiáng)，加劇了污染物在地下水系統(tǒng)不同層間的傳播和運(yùn)移，使問題進(jìn)一步惡化［9］。這種變化會(huì)影響河流的水利狀況，反過來也會(huì)加劇水質(zhì)污染問題，因此，需要針對(duì)河流含水層系統(tǒng)的管理采取更加綜合有效的措施。

確定產(chǎn)生河流和含水層系統(tǒng)流動(dòng)與運(yùn)移動(dòng)態(tài)的環(huán)境，對(duì)解決河流污染問題至關(guān)重要。然而，目前人們對(duì)此過程了解甚少，而且缺乏適當(dāng)?shù)哪Ｐ蛠砟M河流與含水層間的相互作用，無法獲得污染物的運(yùn)移情況。因此，本文提出了一種基于河流與地下水空間耦合分布的河流走廊污染物運(yùn)移模型，該模型克服了現(xiàn)有單一系統(tǒng)建模方法和集水區(qū)尺度水文模型的局限性，真實(shí)模擬了河流走廊中河流和含水層之間的相互作用及污染路徑，有利于調(diào)查現(xiàn)有污染的影響，為提出相應(yīng)的管理策略提供依據(jù)。

2 模型的建立

2.1 典型模型分析

河流和地下水交換流量與質(zhì)量的速率是隨空間和時(shí)間變化的，取決于界面處的水力和水質(zhì)特征。由于流量和質(zhì)量的交換速率是動(dòng)態(tài)的，因此會(huì)出現(xiàn)非線性交換模式。

為了描述該現(xiàn)象，需要建立耦合地表水（SW）和地下水（GW）成分的水文模型。已有的模型包括MODHMS，InHM，HydroGenoSphere，MIKE－SHE，tRIBS。但是，上述模型的比例尺雖然可以通過假設(shè)來簡(jiǎn)化陸上水流路徑，但對(duì)于河段和走廊比例尺研究而言，此類簡(jiǎn)化仍然是不合適的，尤其是不適合模擬洪水和湍流等現(xiàn)象以及相關(guān)的輸運(yùn)過程，因?yàn)樗鼈冃枰敿?xì)的信息。

隨著計(jì)算機(jī)的進(jìn)步，耦合地表非穩(wěn)定水動(dòng)力和地下水流動(dòng)解算器的模型逐漸成熟，包括SWIFT2D－SEWAT，MODFLOW LGR VSF newCFL，MODFLOW－OWHM 和2dMb MODFLOW。

但是，二維地表水流模型和三維地下水流模型的完全耦合仍處于起步階段，均沒有考慮洪水和濕地研究的移動(dòng)邊界問題，仍未在輸運(yùn)和水質(zhì)耦合方面取得進(jìn)展。

目前有多種典型的地表水和地下水模型，通常只適用于模擬河流或含水層的動(dòng)態(tài)情況。但是，由于以下原因，它們不適用于涉及瞬態(tài)條件、不規(guī)則形狀河床的復(fù)雜情況下的河流含水層流量和質(zhì)量估算：

（1）大多數(shù)是為解決典型河流或含水層問題而開發(fā)的單一系統(tǒng)模型。雖然針對(duì)影響系統(tǒng)的過程進(jìn)行了相應(yīng)的建模，與其他非建模系統(tǒng)相互作用產(chǎn)生的過程仍需要一組邊界條件，但是邊界條件的定義在多數(shù)情況下被簡(jiǎn)化并且不符合實(shí)際情況。

（2）河道和漫灘中河流與含水層的相互作用是一個(gè)復(fù)雜的動(dòng)態(tài)過程。相互作用的方向取決于地下水位和河流水位的相對(duì)位置。但是，河流與含水層之間的流量和質(zhì)量交換通常被忽略或簡(jiǎn)化為均勻或線性變化的預(yù)定義速率，這與實(shí)際情況是不符的。

（3）可能發(fā)生界面通量方向的反轉(zhuǎn)，從而導(dǎo)致每個(gè)系統(tǒng)交替充當(dāng)源或匯。因此，河流－含水層界面的物質(zhì)交換的模擬，受不斷變化的壓力梯度和作為源頭的整個(gè)系統(tǒng)的濃度分布的影響，需要并行求解河流與地下水流動(dòng)和輸運(yùn)方程，同時(shí)根據(jù)相互作用的變化不斷更新地表水和地下水系統(tǒng)狀態(tài)。對(duì)于外部耦合河流和地下水模型，一般無法計(jì)算由河流－含水層相互依存的動(dòng)力學(xué)引起的流量振蕩。

2.2 新模型構(gòu)建

2.2.1 模型架構(gòu)

本文所建立的模型主要將兩種運(yùn)移模型（QeS2用于地表水［10］，MT3DMS 用 于 地 下 水［11］）明 確 地耦合到一個(gè)整合的河流－含水層流動(dòng)模型（2dMb＋MODFLOW）中。該模型不僅規(guī)定了流動(dòng)和運(yùn)移方程的聯(lián)立解，還計(jì)算了水和示蹤劑在河流與地下水系統(tǒng)連續(xù)體中的交換通量［12-15］。模型架構(gòu)及模擬過程的簡(jiǎn)化表示如圖1 所示。

圖1 模型架構(gòu)及模擬過程的簡(jiǎn)化

地表水動(dòng)力與運(yùn)移模型部分分別為2dMb 和QeS2 模型的修改版本，主要用于求解二維歐氏空間中的二維深度平均淺水和運(yùn)移偏微分方程。由于地形的離散化用來獲取河流走廊區(qū)域的地表流和運(yùn)移動(dòng)態(tài)，因此無需指定河流的位置。2 個(gè)地下水模型部分分別為MODFLOW 和MT3DMS 的修改版本，主要用于求解二維/三維地下水流動(dòng)和運(yùn)移偏微分方程。每個(gè)流動(dòng)模型組件中所使用的時(shí)間步長(zhǎng)需要獨(dú)立計(jì)算，并且在整個(gè)模擬過程中要不斷更新，以符合Courant（CFL）條件。一般來說，河流模型的時(shí)間步長(zhǎng)通常為秒或1/10 秒，具體取決于網(wǎng)格分辨率和流量條件。而地下水模型所需的時(shí)間步長(zhǎng)通常為小時(shí)或天，比河流模型所需的時(shí)間步長(zhǎng)高出104倍。

2.2.2 污染負(fù)荷方案

河流的污染源類型主要分為點(diǎn)式和分布式兩種［16］。其中，第一種類型主要包括模擬溢出型負(fù)荷的瞬時(shí)點(diǎn)源和模擬出口型負(fù)荷的連續(xù)點(diǎn)源。后者主要包含線性斜率算法方案，因?yàn)樗梢酝ㄟ^逐次線性逼近而適用于幾乎任何負(fù)荷模式［17］。瞬時(shí)點(diǎn)源和連續(xù)點(diǎn)源分別由以下方程式描述：

式中，Km為負(fù)荷率，%，主要取決于移動(dòng)河床特征或自由水面面積。

2.2.3 模型設(shè)置與參數(shù)化

用于測(cè)試（瞬態(tài)）傳輸部分的解析解需要穩(wěn)態(tài)和均勻的流場(chǎng)以及給定的色散系數(shù)，然而，SW 和GW流量部分產(chǎn)生此類條件的邊界條件與參數(shù)的定義并不直接，需要進(jìn)一步考慮。在GW 流量分量的情況下，縱向壓力梯度和水力傳導(dǎo)率被定義為產(chǎn)生0.33 m/s的滲流速度，該速度使用達(dá)西定律估算，并可直接用于解析解。將橫向壓力梯度和電導(dǎo)率分別設(shè)置為0和Kx×10－5，以使流場(chǎng)為單向且沿縱向方向。在基于電導(dǎo)的模型部分中，參數(shù)kc/b（電導(dǎo)層厚度上的電導(dǎo)）值范圍在1～0.001 m/d 之間，在本模型中，使用的值為0.043 2 m/d。為了確保受地表水體滲透影響的含水層中流場(chǎng)的準(zhǔn)確估計(jì)，地下水網(wǎng)格單元分辨率不應(yīng)超過特征長(zhǎng)度，因此，其值設(shè)置為10 m。

初始和邊界條件以及所用污染源的特征見表1。為了便于驗(yàn)證模型，定義了2 種不同的場(chǎng)景“S1”和“S2”，分別測(cè)試SW 和GW 污染物傳輸部分以及SW－GW 相互作用模型部分。場(chǎng)景“S1”包括2 個(gè)連續(xù)的污染負(fù)荷點(diǎn)，1 個(gè)位于河流中，另1 個(gè)位于地下水中，分別測(cè)試SW 和GW 運(yùn)輸模型部分的性能。場(chǎng)景“S2”旨在測(cè)試SW－GW 交互模型，僅使用邊界和初始條件進(jìn)行強(qiáng)制。

表1 初始和邊界條件以及所用污染源的特征

3 結(jié)果與分析

用于測(cè)試傳輸組件的分析解決方案要求在整個(gè)區(qū)域內(nèi)具有穩(wěn)態(tài)和均勻流動(dòng)條件。為了確保基準(zhǔn)問題滿足上述條件，首先測(cè)試模型是否能夠再現(xiàn)所需的流場(chǎng)。

3.1 地下水運(yùn)移

通過模擬結(jié)果可知，形成了均勻、穩(wěn)定和單向的壓力梯度，并如前面所設(shè)計(jì)的一樣形成了0.33 m/s的穩(wěn)定滲流速度。將GW 傳輸部分的數(shù)值結(jié)果與2個(gè)不同時(shí)刻的相應(yīng)解析解進(jìn)行對(duì)比，結(jié)果表明，模型結(jié)果與解析解非常吻合，其中具有微小的差異，主要是由數(shù)值色散引起的，誤差在可接受的范圍之內(nèi)。

3.2 地表水運(yùn)移

基于分布場(chǎng)景“S1”中達(dá)到穩(wěn)態(tài)流動(dòng)條件后的水位、速度、剪切速度、湍流黏度和擴(kuò)散率分布結(jié)果，穩(wěn)態(tài)條件非常接近，只有微小的偏差。由于區(qū)域大小有限及摩擦，河岸附近的速度較小。此外，由于北部和南部邊界條件的影響，導(dǎo)致了渦流擴(kuò)散率和渦流黏度在這些區(qū)域分布也不均勻。但是，基準(zhǔn)問題的范圍涵蓋了較廣泛的區(qū)域，該區(qū)域不受非均勻行為的影響，可根據(jù)分析解決方案對(duì)SW 傳輸組件的性能進(jìn)行測(cè)試。所選區(qū)域的穩(wěn)態(tài)證據(jù)中，所有與流量相關(guān)的變量在該區(qū)域中都是恒定的，并用陰影交叉線突出顯示。只有渦流擴(kuò)散率沿流動(dòng)方向略有變化，但與平均值的解析解中使用的值相對(duì)應(yīng)。將數(shù)值解與解析解進(jìn)行比較，結(jié)果表明，數(shù)值模型很好地捕捉了解析解在空間和時(shí)間上預(yù)測(cè)的總體濃度分布。微小差異是由SW 流分量模擬的驅(qū)動(dòng)流變量的可變性以及數(shù)值解中遇到河岸的羽流引起的，其影響隨著時(shí)間的推移變得更加明顯。與地下水模型部分的誤差類似，可比區(qū)域的絕對(duì)誤差小于1%。

3.3 河流－地下水相互作用：流動(dòng)和運(yùn)輸

基于電導(dǎo)率模型中的輸運(yùn)分量沒有解析解，因此，當(dāng)外部耦合時(shí)，采用耦合模型得到的場(chǎng)景“S2”的結(jié)果與MODFLOW 和MT3DMS 獨(dú)立版本預(yù)測(cè)的結(jié)果進(jìn)行比較。由結(jié)果可知，2 種模型方法模擬的滲透通量相似，表明所建立的新型模型準(zhǔn)確性較高，在邊界附近的預(yù)測(cè)不匹配主要是由于表面流和運(yùn)移模型中邊界條件的影響。

4 結(jié)論

（1）基于河流與地下水空間耦合分布的河流走廊污染物運(yùn)移特征，本文建立了相應(yīng)的數(shù)學(xué)模型，將目前已有的多個(gè)獨(dú)立模型完全耦合到河流的運(yùn)移過程中。建立該模型是為了解決現(xiàn)有單一系統(tǒng)建模方法和現(xiàn)有尺度水文模型的局限性，提供一種分析復(fù)雜現(xiàn)象的方法。

（2）模擬結(jié)果表明，與現(xiàn)有單一系統(tǒng)模型結(jié)果相比，新型河流污染物運(yùn)移模型能夠更準(zhǔn)確地模擬相連河流含水層系統(tǒng)中的污染路徑，并且在當(dāng)前城市高速發(fā)展的背景下，有利于調(diào)查現(xiàn)有污染的影響，并為提出相應(yīng)的管理策略提供依據(jù)。

（3）雖然該模型仍存在一些局限性，但對(duì)城市河流污染的預(yù)防和恢復(fù)有著重要的指導(dǎo)意義，可在進(jìn)一步了解水文系統(tǒng)之間的復(fù)雜相互作用及發(fā)生在河流含水層界面上的水力和生物化學(xué)過程基礎(chǔ)上，不斷更新及優(yōu)化該模型。