基于转移概率随机模拟的DNAPL运移二维和三维数值模拟对比研究

米 东，樊 皓，罗小勇，郑 菲，施小清*，吴吉春

1.表生地球化学教育部重点实验室，南京大学地球科学与工程学院，南京210023 2.长江水资源保护科学研究所，武汉430051

基于转移概率随机模拟的DNAPL运移二维和三维数值模拟对比研究

米 东1，樊 皓2，罗小勇2，郑 菲1，施小清1*，吴吉春1

1.表生地球化学教育部重点实验室，南京大学地球科学与工程学院，南京210023 2.长江水资源保护科学研究所，武汉430051

非均质介质的空间维度变化对重非水相流体（DNAPL）的运移具有重要的影响。在充分考虑地质体的空间连续性、不对称性以及各向异性等特征的基础上，采用基于马尔可夫链的转移概率（transition probability）模型来构建非均质随机场。该文通过TMVOC-MP软件来模拟DNAPL在非均质介质中的运移规律，探讨非均质随机场的水平空间连续性、空间维度变化以及侧向运移过程对DNAPL运移的影响。结果表明，介质的水平空间连续性越好，DNAPL在水平方向的迁移范围越大，在垂向的迁移范围越小；相比于三维模型，二维模型中DNAPL在水平方向的展布更大、在透镜体上的蓄积量更多，在实际应用中以二维模型代替三维模型会加大模拟结果与实际污染情况之间的误差；侧向运移过程削弱了单个平面的非均质性对DNAPL运移的控制，当存在侧向运移时，DNAPL绕过透镜体所运移的距离以及在透镜体上的蓄积量会相应减小。

重非水相液体；非均质；水平空间连续性；空间维度变化；侧向运移

1 引言

随着有机化工产品的广泛使用，大量有机溶剂在炼制、使用和储存过程中会发生泄漏，加之人类的随意排放，使得地下介质中的有机物污染问题不断加剧（崔俊芳等,2003;李宏和Ranjith, 2008;杨明星等,2013）。有机污染物进入地下环境后多以非水相流体（NAPL,Non-aqueous phase liquid）的形式赋存于地下介质中，其中尤以密度大于水的重非水相流体（DNAPL,DenseNAPL）的污染最为严重（Majiand Sudicky,2008）。

DNAPL在地下环境中的运移和分布受多种因素的控制，除自身的物理化学性质外，介质的非均质性对其运移具有重要的影响（施小清等, 2012;张蔚等,2013）。相比传统的基于交叉变差函数（cross-variogram）的地质统计模型，基于马尔可夫链的转移概率（transition probability）模型能更准确的反应地质体的空间结构特征。该方法建立的模型有比较成熟的随机理论支撑，并且能够更加准确的模拟、预测DNAPL的运移和分布规律（何芳和吴吉春,2003;靳萍等,2009）。

目前，国内外对DNAPL的运移和分布规律已进行了大量的研究，并取得了诸多的理论成果（Dekker and Abriola,2000;Pantazidou and Liu，2008;Christ et al.,2010;施小清等,2012;郑菲等, 2015）。然而，相关研究多侧重于对DNAPL泄露速率以及介质非均质性的探讨，涉及水平空间连续性和空间维度变化对DNAPL运移影响的研究较少。忽略水平空间连续性和空间维度变化对污染物运移的影响，会导致模拟结果与实际污染情况之间存在误差（Christetal.,2005），从而难以真实模拟和预测实际场址中的污染情况，加大治理的难度。因此，本文采用TMVOC-MP软件来模拟DNAPL在具有不同水平空间连续性的非均质介质中的运移过程，对比分析DNAPL在二维和三维模型中的运移结果，以此来探讨介质的水平空间连续性和空间维度变化对污染物运移的影响，从而为污染物的修复治理提供科学依据。

2 研究方法

2.1 TMVOC-MP简介

TMVOC-MP是由美国劳伦斯-伯克利国家实验室开发的TMVOC软件的并行版本，主要用于地下介质中环境污染问题的数值模拟研究。该软件采用积分有限差分方法进行空间离散，可模拟多种NAPL（包括挥发性有机物如原油、汽油、柴油以及有机溶剂等流体）在多维非均质介质中的迁移、溶解、吸附以及生物降解等过程，也可以用于分析水-气-NAPL三相之间的相互转化（Pruess and Battistelli,2002;施小清等,2009;刘晓娜等, 2012）。针对多相流控制方程具有强非线性的特点，TMVOC-MP采用Newton-Raphson迭代法求解非线性方程组，并调用并行迭代解法器库AZTEC中的Krylov子空间方法来进行并行计算，显著提高了软件的计算效率（Zhang et al.,2007;米东等, 2016）。

2.2 二维模型和三维模型概述

三维模型是一个长（X）、宽（Y）和高（Z）分别为16m、10m和12m的长方体，其水平方向和垂向的离散间距分别为0.5m和0.15m，共分为51200（32×20×80）个网格。为了满足对比条件，二维模型的所有参数均与三维模型中对应的XZ剖面相同，如图1。采用基于相元的方法将地质体概化为粘土、粉土、粉砂和粉质粘土四种岩性，分别用岩性1、2、3和4来表示（其中岩性2为背景岩性），其所占比例分别为10%、60%、20%和10%。假定研究区内岩性的孔隙度均为0.3，固体颗粒的比重均为2650 kg/m3。模型顶、底部以及前后两个面（XZ剖面）均位于地下水面以下且为零通量边界，左右边界（YZ剖面）为定水头边界。模型先在水-气两相体系中达到平衡，使得研究区内的水相饱和度达到1.0，假定区域内的水力梯度为0.1，水流从左向右流动。选取四氯乙烯（PCE）为DNAPL污染物的代表，其在顶部以下0.975m处以面源方式（NX=4）×（NY=5）均匀注入，每个注入点的注入速率为5.54 kg/d。模型的其他主要参数设置见表1。

图1 概念模型示意图Fig.1 Schematic diagram showing the 2-d and 3-d injectionmodels

2.3 非均质随机场

非均质介质的水平空间连续性对DNAPL运移具有重要的影响。在充分考虑地质体的水平空间连续性的基础上，采用地质统计软件T-PROGS（Transition Probability Geostatistical Software）生成三维非均质随机场。该软件的主要方法是基于马尔可夫链的转移概率模型，其建立的三维非均质模型能够真实反映岩相的空间结构（Carle，1999；Majietal.，2006）。为了探究地质体的水平空间连续性对DNAPL运移的影响，本文构建了三组不同水平相关长度的非均质随机场，其具体参数见表2。随着水平相关长度的增大，非均质介质的水平空间连续性越好，如图2。

表1 模型计算的参数设置Table 1 Hydrogeologicaland other parametersused in the simulation

表2 非均质随机场的参数设置Table2 Parameters for stochastic permeability fields

图2 非均质随机场的一次实现（下图为上图的切片图）Fig.2 One realization ofstochastic permeability fieldwith differentcorrelation length(The chartsbelow are the slicegraphsof theupperdiagrams)

3 模拟结果

3.1 单个实现的模拟结果

从三组非均质随机场中分别选取一个具有代表性的实现，模拟二维模型和三维模型中PCE的运移过程和分布结果。在连续注入80 d后，对比分析PCE在不同相关长度的非均质随机场中的模拟结果，探讨水平空间连续性、空间维度变化以及侧向运移过程对PCE运移和分布的影响。

3.1.1 水平空间连续性的影响

为了探究非均质介质的水平空间连续性对PCE运移的影响，分别模拟三维模型中PCE在不同非均质随机场的运移过程。在连续注入80 d后，PCE在非均质介质中的运移和分布结果如图3所示。

从图3可以看出，在三维模型中，随着非均质介质水平空间连续性的增大，PCE在透镜体上蓄积的量就越多、蓄积长度（面积）就越大，并且其运移到模型底部的量就越小。这是因为水平空间连续性越好，非均质随机场中渗透性较差的透镜体的范围就越大，从而限制了PCE向下运移的过程。该结果表明，水平空间连续性对PCE的运移具有重要的影响。

3.1.2 空间维度变化的影响

在考虑空间维度变化对PCE运移的影响时，要保证除了模型的维度外，其他条件（例如二维模型与三维模型的网格剖分、二维模型与三维模型中对应平面注入的PCE总量、初始条件和边界条件等）均相同。在连续注入80 d后，PCE在非均质介质中的运移和分布结果如图4。

图3 PCE在不同非均质随机场中的模拟结果（Snapl=0.1）Fig.3 The simulation resultsof PCE in differentstochastic permeability fields(The figures show the 0.1 isosurface of PCE saturation)

图4 PCE在二维和三维模型中的模拟结果Fig.4 The simulation resultsofPCE in 2-D and 3-Dmodels

在三维模型的二维剖面上，PCE的浓度值低于二维模型中的浓度值、PCE的运移范围小于二维模型中的运移范围。这是因为PCE在三维模型中运移时存在侧向运移过程，从而造成该剖面上PCE的总量小于二维模型中PCE的总量。该结果表明，在进行实际污染场址的模拟和修复时，应采用三维模型来模拟和预测PCE在非均质介质中的运移过程和分布范围，以二维模型代替三维模型会导致模拟结果与真实污染情况不符。

相比于二维模型，三维模型中PCE在渗透性较差的透镜体上蓄积的长度和浓度均小于二维模型的模拟结果。造成该现象的主要原因是侧向运移过程削弱了单个平面的非均质性对DNAPL运移的控制：在二维模型中，PCE向下运移到渗透性较差的透镜体时只能绕过透镜体，从而造成PCE在透镜体上产生大量的蓄积，并导致其在透镜体上蓄积的长度与透镜体的长度大体相等；而在三维模型中，PCE向下运移到渗透性较差的透镜体时，当侧向运移的距离小于透镜体的长度时，其优先通过侧向运移过程绕过该透镜体，这样不仅缩短了PCE的运移路径，而且显著减小了PCE在透镜体上的蓄积量。

3.1.3 垂向运移随时间的变化

考虑水平空间连续性以及空间维度变化对PCE向下运移速率的影响，选取Z方向的质心位置为指标，分别计算其在PCE连续注入10 d、20 d、30 d、40 d、50 d、60 d、70 d和80 d后的值，从而定量描述PCE向下运移的过程，如图5。

从图5可以看出：

（1）二维模型中，在20 d以前，水平方向的相关长度越大，PCE向下运移的速率越大；20 d后，水平相关长度越小，PCE向下运移的速率越大。这是因为，随着水平方向相关长度增大，透镜体也随之增大，PCE在向下运移过程中遇到透镜体的概率逐渐减小。在水平方向相关长度较大的非均质随机场中，在20 d以前，PCE在向下运移的过程中并未遇到透镜体。而在20 d以后，PCE逐渐在透镜体上蓄积，从而影响了向下运移的速率。

图5 PCE在垂向的运移规律Fig.5 The PCEmigration distance in zdirection versus time

（2）三维模型中，PCE在水平方向相关长度中等的非均质随机场B中运移速度最快，其次为非均质随机场A，最后是非均质随机场C。这是因为在三维模型中，当水平方向的相关长度较小时，PCE多以离散状赋存在运移路径上；而当水平方向的相关长度较大时，PCE运移到较大的透镜体时会发生蓄积，从而影响其向下运移的速率；只有当水平方向的相关长度中等时，PCE运移的通道较为集中，加之侧向运移过程减小了透镜体对PCE运移的控制，使得PCE可以快速运移到模型底部。

（3）在20 d以前，PCE向下运移时没有遇到较大的透镜体，因此二维模型和三维模型的运移结果差别不大。而在20 d以后，随着PCE持续注入，其在向下运移时遇到较大的透镜体，三维模型由于存在侧向运移过程，从而使得单个剖面上PCE的总量少于二维模型PCE的总量，造成PCE在三维模型中向下运移的速率小于其在二维模型中向下运移的速率，并且该现象会随着时间的增加而不断加剧。

3.2 多个实现的模拟结果

由于非均质随机场单个实现具有随机性和不确定性的特点，采用蒙特卡洛方法对每组随机场分别生成50个实现，从统计的角度来研究非均质介质的水平空间连续性和空间维度变化对DNAP L运移和分布的影响。采用空间矩分析方法定量描述污染物的运移行为，选取质心、展布和GTP（定义为连续的离散状污染物的体积除以连续的池状污染物的体积）三个参数来描述PCE的污染源区结构特征，具体如下：空间一阶矩和二阶矩分别表示污染羽的质心位置和污染羽围绕质心的展布范围，采用GTP来探讨污染物在介质中的存在形态。

3.2.1 质心位置

当PCE在运移过程中遇到渗透性较差的透镜体时，会产生蓄积和绕流现象，随着水平方向相关长度的增大，透镜体的范围增大，PCE要绕过透镜体所运移的距离也随之增大。因此，随着水平方向相关长度的增大，PCE在X方向的运移距离逐渐增大、在Z方向的运移距离逐渐减小（图6）。

PCE在三维模型的二维剖面上的运移范围小于其在二维模型中的运移范围，其质心的离散程度也远小于其在二维模型中的离散程度。这是因为在三维模型中存在侧向运移过程，从而造成PCE在三维模型中单个剖面上的总量小于二维模型中的总量。同时，侧向运移过程会削弱单个平面上非均质性对PCE运移的影响，使得三维模型中质心的离散程度远小于二维模型中质心的离散程度。

图6 PCE在X、Z方向质心位置的箱形图Fig.6 The box-p lotsof PCEmass centroid in x direction and z direction

3.2.2 展布

随着水平方向相关长度的增大，PCE在水平方向的展布逐渐增大、在垂直方向的展布逐渐减小，这表明水平空间连续性对PCE在非均质随机场的运移具有重要的影响。同时，由于侧向运移过程不仅会减少三维模型中二维剖面上PCE的总量，而且会削弱单个平面的非均质性对PCE运移的影响，因此PCE在三维模型中二维剖面上的展布及其离散程度远小于其在二维模型中的展布及其离散程度（图7）。该结论与3.4.1的结论一致。3.2.3 GTP

在连续注入80 d后，PCE在水平相关长度中等的非均质随机场中的GTP值最大（图8）。这是因为，当水平相关长度较小时，非均质介质中的透镜体较小，PCE可以快速运移到模型底部并形成池状污染；而当水平相关长度较大时，非均质介质中的透镜体较大，PCE绕过透镜体的所需的距离也随之较大，从而在透镜体上会产生大量蓄积并形成污染池；只有当水平相关长度中等时，PCE即不会快速运移到模型底部，也没有在透镜体上产生大量蓄积，而是多以离散状赋存在非均质介质中。

对比PCE在二维模型与三维模型中的GTP值时可以发现，三维模型中二维剖面上GTP的值大于二维模型GTP的值，造成该现象的主要原因是三维模型中存在侧向运移过程。由于侧向运移过程会削弱单个平面的非均质性对PCE运移的影响，使得PCE在透镜体上的蓄积现象减弱，从而使得PCE更多的以离散状态赋存在非均质介质中。

图7 PCE在X、Z方向展布的箱形图Fig.7 Thebox-plotsofspatialdistribution of PCE in x direction and zdirection

图8 GTP的箱形图Fig.8 Thebox-plotsofGTP

4 结论

（1）非均质介质的水平空间连续性会影响透镜体的大小，从而改变DNAPL的运移路径和分布范围。非均质介质水平空间连续性越好，透镜体的范围就越大，DNAPL运移的距离以及在透镜体上的蓄积量也随之增大。

（2）空间维度变化对DNAPL的运移和分布具有重要的影响。DNAPL在二维模型中运移的范围会大于其在三维模型中对应剖面上的运移范围，所以使用二维模型来模拟和预测DNAPL在非均质介质中的运移和分布结果，将会增大数值模拟结果与真实污染情况之间误差。

（3）侧向运移过程削弱了单个平面的非均质性对DNAPL运移的控制。在模拟和预测DNAPL污染运移和分布结果的过程中，当存在侧向运移过程时，DNAPL绕过透镜体所运移的距离以及在透镜体上的蓄积量会相应减小。

崔俊芳，郑西来，林国庆.2003.地下水有机污染处理的渗透性反应墙技术[J].水科学进展，14(3):363-367.

何芳，吴吉春.2003.基于马尔可夫链的多元指示地质统计模型[J].水文地质工程地质，30(5):28-32.

靳萍，邵景力，李长青，等.2009.基于T-PROGS的地下水三维数值模拟及应用[J].水文地质工程地质，36(4):21-26.

李宏，Ranjith PG.2008.地下水中非水相流体污染运移的三维有限元法[J].岩石力学与工程学报，27(12):2509-2516.

刘晓娜，程莉蓉，张可霓，等.2012.地下水LNAPL层的原位曝气模拟研究[J].环境科学与技术，35(2):19-24.

米东，施小清，张可霓，等.2016.基于MPI的TOUGHREACT并行性能分析[J].水文地质工程地质，(1):34-40.

施小清，姜蓓蕾，吴吉春，等.2012.非均质介质中重非水相污染物运移受泄漏速率影响数值分析[J].水科学进展，23(3):376-382.

施小清，张可霓，吴吉春.2009.TOUGH2软件的发展及应用[J].工程勘察，(10):29-34.

杨明星，杨悦锁，杜新强，等.2013.石油污染地下水有机污染组分特征及其环境指示效应[J].中国环境科学，6:1025-1032.

张蔚，施小清，吴剑锋，等.2013.渗透率空间变异性对重非水相流体运移的影响[J].高校地质学报，(4):677-682.

郑菲，高燕维，施小清，等.2015.地下水流速及介质非均质性对重非水相流体运移的影响[J].水利学报,46(8)：925-933.

Carle SF.1999.T-PROGS:transition probability geostatisticalsoftware,usersmanual version 2.1[R].UniversityofCalifornia,Davis.

Christ J A,Lemke L D and Abriola L M.2005.Comparison of two-dimensional and three-dimensional simulations of dense non-aqueous phase liquids(DNAPLs):Migration and entrapment in a nonuniform permeability field[J].Water Resources Research,41(1): 205-230.

Christ JA,Ramsburg CA,PennellK D,etal.2010.Predicting DNAPLmass discharge from pool-dominated source zones[J].Journal of Contaminant Hydrology,114(1-4):18-34.

Dekker T Jand Abriola LM.2000.The influenceof field-scaleheterogeneity on the surfactant-enhanced remediation of entrapped non-aqueous phase liquids[J].JournalofContaminantHydrology,42(2):219-251.

MajiR and Sudicky EA.2008.Influence ofmass transfer characteristics for DNAPL sourcedepletion and contaminant flux inahighly characterized glaciofluvialaquifer [J].JournalofContaminantHydrology,102(1-2):105-119.

Maji R,Sudicky E A,Panday S,et al.2006.Transition probability/Markov Chain analyses of DNAPL source zones and plumes[J].Ground Water, 44(6):853-863.

Pantazidou M and Liu K.2008.DNAPL distribution in the source zone: Effect of soil structure and uncertainty reduction with increased sampling density[J].Journal of Contaminant Hydrology,96(1-4): 169-186.

Pruess K and Battistelli A.2002.TMVOC,a numerical simulator for three-phase non-isothermal flows of multicomponent hydrocarbon mixtures in saturated-unsaturated heterogeneousmedia[R].CA:Berkeley:Lawrence Berkeley NationalLaboratory.

Zhang K,Yamamoto H and Pruess K.2007.TMVOC-MP:a parallel numerical simulator for three-phase Non-isothermal flows ofmulticomponent hydrocarbon Mixtures inporous/fractured media[R].CA:Berkeley:Lawrence Berkeley National Laboratory.

Com parison of Two-dimensionaland Three-dimensional Simulations of DNAPL M igration in Saturated PorousM edia

MIDong1,FAN Hao2,LUO Xiaoyong2,ZHENG Fei1,SHIXiaoqing1*,WU Jichun1
1.Key LaboratoryofSurficialGeochemistry,MinistryofEducation;SchoolofEarth Sciencesand Engineering,Nanjing University, Nanjing210023,China; 2.ChangjiangWaterResourcesProtection Institute,Wuhan430051,China

The cChange in the spatial dimension of heterogeneousmedia significantly influences themigration and distribution of Dense Non-Aqueous Phase Liquid(DNAPL).Based on the Markov Chain transition probability model,we generate stochastic permeability fields and take a fully account of horizontal spatial continuity,asymmetry,and anisotropy of heterogeneousmedia.By simulating themigration and distribution of DNAPL in heterogeneousmediawith TMVOC-MP,we investigate the effects of horizontal spatial continuity,reduced dimensionality,and lateral migration on the migration and distribution of DNAPL.Numerical analysis reveals thata better horizontal spatial continuity results in a larger horizontalpollution scopeand a smaller verticalmigration extentof DNAPL.Comparedwith thebehavior in the three-dimensionalmodel,DNAPL distributes in a larger scope in the horizontal direction, accumulatesmore on the low-permeability lenses,and ismostly in a pool state in the two-dimensionalmodel.Lateralm igration processcan reduce the impact of heterogeneity of a single plane on the DNAPL migration.So it is necessary to use the threedimensionalmodelin practicalapplications instead of two-dimensionalmodel.

DNAPL;heterogeneity;horizontal spatial continuity;physical dimensionality;lateralm igration

SHIXiaoqing,Associate Professor;E-mail:shixq@nju.edu.cn

P641

A文献标识码：1006-7493（2016）04-0733-08

10.16108/j.issn1006-7493.2016061

2016-04-24；

2016-11-01

国家自然科学基金-新疆联合基金（U1503282）；国家自然科学基金（41030746；41172206）联合资助

米东，男，1990年生，硕士研究生，主要从事地下水数值模拟研究；E-mail:midong1990@163.com

*通讯作者：施小清，男，1979年生，副教授，博士，主要从事地下水数值模拟研究；E-mail:shixq@nju.edu.cn