解析法和數值法在地下水溶質運移預測評價中的對比研究

吳 瑾 賈冬梅 李學寧

(1.天津市地質環境監測總站，天津 300191；2.天津市地質調查研究院，天津 300191)

2011年 6 月，環境保護部頒布了《環境影響評價技術導則—地下水環境》(HJ 610-2011)，對建設項目環境影響評價中的地下水環評工作起到了指導與規范作用。導則中明確要求一級評價應采用數值法，二級評價中水文地質條件復雜且適宜采用數值法時建議優先采用數值法，三級評價可采用解析法或類比分析法。因此，數值法和解析法在地下水環境影響評價的預測分析中都被廣泛應用。

解析法能夠給出控制微分方程的精確解；數值法采用一組代數方程近似處理微分方程[1]。當水文地質條件比較簡單而且己經取得足夠的水文地質資料和彌散參數資料時，可采用解析法。當含水層的邊界條件、結構和水文地球化學條件相對復雜、含水層之間存在越流時，數值法能夠靈活處理復雜的邊界條件、含水層的非均質性和各向異性等，此時宜選用數值法進行預測[2]。

數值法和解析法都有各自的優缺點。數值法的優勢是可以容易的解決水流和運移參數(如水力傳導系數、孔隙度、彌散度、陽離子交換量等)的分布異質性問題[3]。數值法比解析法靈活性更大，所以數值法可直接用于解決野外場地具體問題，模擬運算結果更能貼近實際。但數值法獲取專業的水文地質參數較為困難，且模型建立對數據的數量質量要求嚴格、操作繁瑣、周期較長，另外，它會導致所得的結果擴展超過與建模主體示蹤劑的彌散無關的溶質鋒面或污染羽[4]，這些都阻礙了數值模擬應用的范圍及廣度。解析法雖然對參數要求不高，操作起來也更加簡便易行，但解析法在使用過程中受地質條件限制較為明顯，鑒于區域因素、人為因素、環境因素的影響，在預測結果準確性上會有一定影響[5]。另外解析法僅能對一個投源點加以解析，在需要對面源進行分析時，具有一定局限性，預測結果精度不大[6]。

通過對數值法和解析法優缺點的對比，本文以某油田井場為例，分別用上述兩種方法進行溶質運移預測分析，通過預測結果的比對，探索兩種預測結果的差異性及對地下水環境影響評價預測分析的適宜性。

1 研究區水文地質概況

研究區位于天津市濱海新區某油田井場內，地勢西高東低，多洼淀或坑塘水域，處在我國典型的淤泥質海岸岸段北部渤海灣西岸，地貌類型屬于海積平原地貌。研究區場地范圍主要由濱海瀉湖洼地構成，地表以粘性土為主，土壤鹽漬化嚴重，并保留有貝殼堤，為距今200～5000年海岸變遷的遺跡。

第I含水組分為潛水和微承壓水，底界埋深約為70～80 m，含水層以粉細砂為主，一般4～5層，累計厚度10～20 m，東部最厚可達40 m。含水組的富水性弱，涌水量東部100～500 m3/d，西部多小于100 m3/d。咸水礦化度一般10～20 g/L。水化學類型為Cl-Na型。淺層多為咸水或咸鹵水，水質差，大部分地區均為不開采。研究區地下水水位埋深較淺，場地內潛水主要靠大氣降水入滲補給、地表水體入滲為主。地下徑流主要是自東南向西北方向。場地內地下水排泄方式為潛水蒸發、側向流出。

第II含水組底界埋深約為175～185 m，含水層以粉細砂為主，偶見粗砂，一般8～9層，單層厚度2～5 m，最厚約10 m。累計厚度北部40～50 m，中、南部27～36 m。其富水性由北向南變差，北部永定新河以北涌水量2 000～3 000 m3/d，向南至塘沽區中北部一帶，涌水量在1 000～2 000 m3/d，導水系數100～300 m2/d。塘沽區東部和南部廣大地區涌水量小于500 m3/d，導水系數50～100 m2/d。咸水的底界埋深在海河以北70～110 m，向南由110 m漸增至210 m，南部第II含水組為咸水。第II含水組總體上為淡水，北部礦化度0.4～0.9 mg/L，化學類型為HCO3-Na型，向南過渡為HCO3·Cl-Na和Cl·HCO3-Na型，礦化度0.7～1.0 mg/L，局部集中開采區地下水礦化度增高，有水質惡化趨勢，礦化度增高到2.21 mg/L。本含水組是塘沽區主要開采層之一。

第III含水組底界深度約為275～285 m，含水層以細砂、粉細砂為主，偶見中砂，一般6～8層，單層厚度3～6 m，累計厚度36～43 m，向南變薄。其富水性由北向南變差。東北部涌水量在2 000～3 000 m3/d和1 000～2 000 m3/d，導水系數100～300 m3/d，向南至海河以北變為500～1 000 m3/d，海河以南多小于500 m3/d。礦化度由北向南由0.6 g/L增至1 g/L左右，水化學類型由HCO3-Na過渡為HCO3·Cl-Na型和Cl·HCO3-Na型。本含水組也是塘沽區主要開采層之一。

第IV含水組底界深度約為405～415 m，下部包括部分新近系含水層。含水層巖性以粉砂、細砂為主，偶見中砂。北部單層厚度4～6m，累計厚度40～50 m，向南變薄為30～40 m。本組富水性較差，除西部涌水量大于2 000 m3/d 外，其余大部分地區在500～1 000 m3/d，向南部富水性更差，多小于500 m3/d。礦化度在0.4～0.7 g/L，以HCO3-Na和HCO3·Cl-Na型為主[7]。

2 數值法對溶質運移預測的應用

2.1 數值模型的建立

研究區采用對流—彌散方程來描述石油類污染物在三維地下水流中的運移，溶解于地下水中NAPLs運移的數學模型可表示為：

式中：n—含水層介質有效孔隙度；c—污染質濃度(mg/L)；Dx—縱向彌散系數 (m2/s )；Vx—縱向滲流速度(m/s)；t—時間(s)；M—生物降解項；φ—模型修正系數，根據前人研究成果，取模型修正系數為0.92。研究區水文地質結構模型和參數選取分別如圖1和表1所示。

圖1 水文地質結構模型圖

表1 參數匯總表

研究區溶質運移模型采用三維非穩定流數學模型：

(x,y,z)εG,t>0

C(x,y,z,t)lt=0=C0(x,y,z)∈G

C(x,y,z,t)lΓ1=C1(x,y,z,t),(x,y,z)∈Γ1

式中：C為濃度；Dij為水動力彌散系數張量的分量；W為源匯項；Cw為隨W注入的污染物濃度；Rd為滯留因子；n為孔隙度；ui為水流實際速度分量；G為研究區域;Γ為研究區域的邊界。

研究模型在垂向上概化為兩層。第一層為潛水含水層(-17～0 m)，巖性主要為粉質粘土、淤泥質粉質粘土、粉土；第二層為弱透水層(-22～-17 m)，巖性主要為粉質粘土。

研究區以井場為例，采用有限差分法，在平面上對井場進行矩形剖分，共剖分了230行200列，共計46 000個單元格，其中活動單元格15 187個，單元格長10 m，寬10 m，見圖2。

圖2 模擬區平面網格剖分圖

由于水文地質概念模型、模型邊界條件概化、源匯項及水文地質參數處理、模型識別驗證及調參、污染源強設定不是本文主要研究對象，這里不再對上述過程進行詳細論述。

2.2 預測結果

本次模擬選取的預測因子為石油類，石油類的標準限值參照《地表水環境質量標準》(GB3838—2002)中的Ⅲ類標準。當污染物濃度大于標準限值時，表示地下水受到污染，以此計算污染超標距離；當污染物濃度小于標準限值并大于檢出限時，表示地下水受到污染的影響，但不超標，以此計算污染影響距離；當預測污染物濃度小于檢出限時視同對地下水環境基本沒有影響。根據預測結果，地下水中石油類在100天、1000天、30年、50年時向地下水下游方向最大超標距離見表2和圖3。

100天石油類污染物運移范圍圖

表2 地下水中石油類遷移距離預測結果表(數值法)

3 解析法對溶質運移預測的應用

3.1 解析法原理

假設石油類泄露排放形式簡化為點源，在投源點瞬時投入質量為M的示蹤劑，地下水沿x軸正方向均勻等速流動，不考慮污染物在含水層中的吸附、揮發、生物化學反應等，且模型中所賦各項參數予以保守性考慮，污染物在潛水含水層中的遷移，可概化為瞬時注入示蹤劑(平面瞬時點源)的一維穩定流動二維水動力彌散問題，污染物濃度分布模型如下：

式中：x，y：計算點處的位置坐標；t：時間(d)；C(x，y，t)：t時刻點x，y處的示蹤劑濃度(g/L)；M：含水層的厚度(m)；mM：瞬時注入的示蹤劑質量(kg)；u：水流速度(m/d)；n：有效孔隙度(無量綱)；DL：縱向x方向的彌散系數(m2/d)；DT：橫向y方向的彌散系數(m2/d)；π：圓周率。

模型需要的主要參數包括：含水層厚度M=17.8 m；外泄污染物質量mM=680.4 g；巖層的有效孔隙度n=0.07；滲透系數K=0.13 m/d；地下水水力坡度I=1‰；水流速度u=0.001 86 m/d；污染物縱向彌散系數DL=0.011 5 m2/d；污染物橫向彌散系數DT=0.004 6 m2/d，這些參數可以由研究區水文地質勘察及類比區域收集成果資料來獲得，計算經過不在此處贅述。

3.2 預測結果

同樣選取石油類作為預測因子，石油類的標準限值參照《地表水環境質量標準》(GB3838—2002)中的Ⅲ類標準。預測結果及向地下水下游方向最大超標距離圖見表3及圖4。

表3 地下水中石油類遷移距離預測結果表(解析法)

圖4 地下水中石油類在不同時間污染暈運移圖(解析法)

4 不同預測方法的預測結果對比分析

根據預測結果，對于解析法，污染物向地下水下游擴散最大超標范圍呈規則的橢圓形，而數值法則呈不規則橢圓形，主要是因為解析法中假定地下水為均勻穩定流，而數值法可通過實際水位和計算水位的擬合分析，使二者的擬合效果基本趨于一致，從而使地下水流場基本符合地下水的實際運動特征[8]。

通過兩種預測方法的結果比較，發現兩個特點：(1)短時間內解析法預測結果小于數值法預測結果，達到最長預測年限50年時解析法預測結果超過數值法。(2)數值法在不同時間段內結果變化程度較解析法小。

針對特點(1)主要原因為：數值法模擬了污染物從地面下滲到含水層中然后再不斷擴散這一整個過程。而解析法在等厚、均值、各項同性、無限長邊界、穩定流邊界的理想化條件下，在預測過程中假定污染物直接在含水層中注入，不考慮下滲過程，也不考慮土壤及微生物對部分污染物的吸附、降解作用，只考慮污染物在地下水流場作用下的對流彌散作用，故污染物在向下游運移的過程中，較數值模擬方法在縱向擴散的速率更快，污染物的運移距離更遠[9]。但在短時間內污染物傳播距離遠的特點難以體現。

針對特點(2)主要原因為：在數值模擬過程中設定了實際的邊界條件并將預測區域按照單元格進行拆分，模擬情況更接近污染物實際運移，污染物的運移速率整體緩慢，隨著時間推移垂向上深度越大預測結果較解析法相比越小，不同時間段內污染物運移距離變化程度不如解析法的變化程度大。

鑒于本項目預測區域內水力參數各項異性、水流流場在廠區內也不能保持完全均一性，導致解析法對實際反映的準確性稍差。同時本項目工程運營時間較長(最長時間50年)，預測時間跨度較大也導致解析法運移距離變化較大。所以針對本項目解析法計算精度較數值法來看仍不太理想。但是各種預測方法都有其各自的適用范圍。解析法作為地下水溶質運移預測方法的一種，為地下水污染預測提供了一種途徑，為相關類似建設項目在地下水環境影響評價中的具體應用提供了科學依據，豐富了地下水環境影響評價方式，對及時做好地下水污染防治工作提供了一定依據。