0.   引言

近些年来中国的现代化工业发展迅速，也很大程度上带动了城市化进程。然而，快速的人口集聚效应给环境也带来了沉重的负担^[1-2]，自20世纪90年代起国内外广泛采用建造垃圾填埋场来对居民生活产生的垃圾进行填埋处理^[3]。填埋场中的污染对于地下水资源及周围环境而言是一个潜在的威胁，且多数填埋场由于服役时间久而已经出现渗漏等问题。因此，有必要对填埋场中的污染物迁移进行深入研究，为采取相应的措施减少其对环境危害提供参考。

在填埋场中，由于各种被填埋的生活垃圾长期受到雨水等介质的影响而发生降解，产生各种有机物以及重金属离子等有毒物质^[4-6]。为了避免这些有毒物质迁移至周围的土体或地下水中，在填埋场底部一般设置衬垫来进行阻隔，如高密度聚乙烯（HDPE）土工膜、膨润土防水毯（GCL）或压实黏土层（CCL）等组合形成的复合衬垫^{[3, 7-8]}。然而，在实际工程中，土工膜由于其厚度较薄而极容易受到外力或不均匀沉降的影响而破损，导致阻隔效果急剧下降^{[3, 8-10]}。

目前针对污染物迁移问题已有大量研究，其中相关理论解析模型大多是基于Van等^[11-13]提出的一维解析模型进行展开的。如污染物在通过具有破损土工膜衬垫并进入下层土体的模型，给出了相应的实验、数值与解析结果^{[8-10, 14-20]}。然而这些解答均局限在一维的情况下，不能准确地描述污染物在土中的二维迁移特性。另外，Chen等^[21]通过假设土层顶部污染源为圆形的情况下给出了极坐标下的解析解，但该模型本质上假设顶部的浓度边界为同一类边界，不能很好地模拟土工膜及破损处的浓度边界情况^[1]。事实上土工膜破损处的土体与污染物直接接触，则该边界应为与浓度相关的第一类边界；而土工膜与土体接触的范围应该对土层顶部产生了类似于Robin条件的第三类边界条件。这样土层顶部的边界应由第一类与第三类边界共同组成，形成了一种混合边界^[22-23]，且以往的研究均忽略了填埋场中由于土工膜破损形成的混合边界。

本文考虑了垃圾填埋场中广泛存在的土工膜破损问题，以及填埋场中边界处浓度随时间变化的情况，建立了污染物透过土工膜并进入下层土体的二维迁移模型。通过积分变换（Laplace变换、Fourier变换及相应的逆变换）及边界转换法对控制方程及边界条件进行半解析解求解，采用数值解方法与得到的解析解结果进行对比，验证了本文研究的可靠性。基于该半解析解，对土工膜及土层参数进行研究，对浓度函数参数、破损率、各向异性参数、厚宽比进行了影响性的研究，得到了不同情况下污染物在土中的迁移规律。该研究可为垃圾填埋场的设计、土中污染物浓度发展趋势及相关污染防治措施等方面提供理论参考。

1.   理论模型

建立的模型包含土工膜及下伏土层的污染物二维迁移，如图 1所示。其中图 1（a）为填埋场中的研究单元，图 1（b）为模型截面图。坐标原点在土工膜缺陷底部的中心处，z为深度方向，x为水平方向。该模型由污染物、土工膜及土层构成，污染物的水头高度为h_w，土工膜的厚度为L_g，土层厚度为L_s；模型的宽度为2L，缺陷宽度为2D。在模型中，土工膜顶部及缺陷内受到渗沥液的污染，且污染物可通过在土工膜中扩散及缺陷内对流与扩散至下方土层。污染物在土层中进行二维扩散，在深度方向由于土工膜缺陷的存在发生对流作用。

图  1  考虑土工膜破损的污染物二维迁移

Figure  1.  Two-dimensional migration model for contaminants through defected GM: (a) research unit, and (b) sectional scheme

下载: 全尺寸图片 幻灯片

1.1   基本假设

基于图 1中的污染物迁移模型，做了如下假设：

（1）忽略土中污染物的化学反应、衰变等过程。

（2）污染物的水头稳定且较低，认为其对流速的影响忽略不计且保持恒定。

（3）土工膜中仅存在扩散作用，污染物在破损处及土层中存在对流作用，在土层中进行二维扩散。

（4）污染物溶液在土中的对流符合达西定律。

1.2   控制方程

污染物在土工膜及土层中的迁移控制方程如下。

土工膜中的污染物迁移控制方程（$- {L_{\text{g}}} \leqslant z \leqslant 0$）：

土层中的污染物迁移控制方程（$0 < z \leqslant {L_{\text{s}}}$）：

式中：t为时间；C_g，C_s分别为土工膜与土层中的浓度；R_d为阻滞因子，R_d=1+ρ_dK_d/n_s，ρ_d为土的干密度，K_d为反映土中吸附能力的分布系数，n_s为土的孔隙率；D_g为土工膜的扩散系数；D_sx，D_sz分别为土体的水平与竖直方向的有效扩散系数；v_s为污染物溶液在土体中的流速；L_g，L_s分别为土工膜与土层的厚度。

1.3   初始和边界条件

（1）初始条件

在初始阶段，土工膜及土层中的污染物浓度假设为0，即

（2）边界条件

土工膜顶部（$z = - {L_{\text{g}}}$）：

土层顶部（$z = 0$）：

土层底部（$z = {L_{\text{s}}}$）：

模型水平方向（$x = 0$与$x = L$）：

式中：g(t)为任意与时间及污染源浓度衰减相关的函数；K_g为土工膜上、下表面的分配系数。

1.4   边界转化

在实际工程中，主要关注土层中的污染物浓度的变化情况。由于土工膜较土层而言很薄，因此根据式（1），（4）~（7），z=0处的边界可转化为^[24]

2.   半解析解推导与验证

2.1   归一化处理

引入无量纲参数C=C_s/C₀，Z=z/L_s，X=x/L及T=v_st/(R_dL_s)，式（2）可转化为以下的无量纲控制方程：

式中：$\eta$为厚宽比，且$\eta = {{{L_{\text{s}}}} / L}$；$P = {{{v_{\text{s}}}{L_{\text{s}}}} / {{D_{{\text{sz}}}}}}$；$\kappa$为土层中的扩散系数比值，且$\kappa = {{{D_{{\text{sx}}}}} / {{D_{{\text{s}}z}}}}$。

另外，无量纲下的边界条件为

式中：$\lambda$为破损率，且$\lambda = D/L$；$\zeta$为体现土工膜材料阻隔污染物能力的参数，且$\zeta = {D_{{\text{sz}}}}{L_{\text{g}}}/({K_{\text{g}}}{D_{\text{g}}}{L_{\text{s}}})$。

2.2   半解析解推导

由于混合边界的特殊性，以下采用Laplace变换和Fourier变换对归一化后的控制方程进行半解析求解。首先对式（11）~（16）进行Laplace变换可得

式中：$\bar C(X, Z, s)$为Laplace变换域下的$C(X, Z, t)$，定义为$\bar C(X, Z, s) = \mathop \smallint \nolimits_0^\infty C(X, Z, T){{\rm{e}}^{ - sT}}{\rm{d}}T$；同样地，有$\bar g(s) = \mathop \smallint \nolimits_0^\infty g(T){{\rm{e}}^{ - sT}}{\rm{d}}T$。其次，根据模型在水平方向上的边界条件，对控制方程及竖向边界条件进行Fourier余弦变换，则式（17）~（19）可变换为

式中：$\tilde {\bar C}(m, Z, s) = \mathop \smallint \nolimits_0^1 \bar C(X, Z, s){\text{cos(}}{M_m}X{\text{)}}{\rm{d}}X$；${M_m}$为Fourier变换中的变量，且${M_m} = m{\rm{ \mathsf{ π}}}$，$m=\text{0}，\text{1}，\text{2}，$$\text{3}，\cdots$；$\mu _m^2 = sP + \kappa {\eta ^2}M_m^2$。

求得式（21）的通解为

式中：${A_{1m}}$，${A_{2m}}$为待定系数，且可通过代入竖向边界条件求出；$\varOmega = {P / 2}$，${\varphi _m} = \sqrt {{\varOmega ^2} + \mu _m^2}$。

将式（22），（23）代入至式（24）中，可得

因此，将式（25），（26）代入式（24），即可得到Laplace变换及Fourier变换下的浓度变化半解析解：

式中：

根据Fourier逆变换理论，可得到浓度变化在Laplace变换域下的解答：

式中：

需要说明的是，式（28）中的结果包含了未知函数$\bar f(X, s)$，无法直接对其进行计算。对此，本文通过对土层顶面与土工膜接触的部分进行离散化^[23]，再通过$\bar f(X, s)$的定义进行求解。

2.3   验证

采用数值解（numerical solution，NS）对本文的半解析结果（current solution，CS）进行验证，且破损处的浓度函数为$g(T) = 1$（即破损处浓度保持恒定）。在验证及后文的分析中，假设污染物初始浓度C₀=100 μg/L。模型中的尺寸：L_g=1.5 mm，D=0.1 m，L=0.4 m，L_s=1 m；本研究不讨论土体的吸附效应，则R_d设置为1。其他参数见表 1^{[3, 16, 25]}。

表  1  土工膜与土层的物理力学参数

Table  1.  Physical and mechanical parameters of GM and SL

参数	GM	SL
Dsx/(m2·s-1)	—	8.9×10-10
Dsz /(m2·s-1)	6×10-13	8.9×10-9
vs/ (m·s-1)	—	5.62×10-10
Kg	0.032	—
注：如文中未作特殊说明或讨论，所取参数保持不变。

下载: 导出CSV 

| 显示表格

对土层中不同空间位置处随时间变化的浓度验证结果如图 2所示。其中，图 2（a）为固定水平位置X=0.5，竖向位置为Z为0，0.5，0.75，1时的浓度曲线；图 2（b）中给出的是固定竖向位置Z=0.5，水平位置为X为0，0.5，0.75，1时的浓度曲线。从图 2可看出，本文解与数值解的结果完全一致。因此可判定本文结果的是准确且可靠的。另外基于文献[15，18，25]中对于破损处边界的处理方法推导污染物一维迁移模型，并将其结果与本文结果进行对比。结果显示，相同条件下一维模型中的浓度迁移比二维情况下的速度快。图 2（a）中，Z=0，X=0.5处一维模型的相对浓度在极短的时间内达到了0.7，可见一维模型高估了污染物在水平位置的迁移速度。图 2（b）中不同水平位置处的一维模型结果完全重合，这是因为一维模型中不同竖向位置处的污染物浓度在水平方向上是相等的。

图  2  不同位置处污染物浓度随时间变化对比

Figure  2.  Comparisons of concentration variations at different positions: (a) fixed horizontal and (b) fixed longitudinal points

下载: 全尺寸图片 幻灯片

3.   参数敏感性分析

为了探讨填埋场中污染物通过破损土工膜迁移至下伏土层的二维特性，采用浓度函数为g(T)=e^-aT（a为函数变化参数，且a > 0）^[25]，对不同参数（函数参数a，土工膜破损率$\lambda$，各向异性参数κ及模型厚宽比η）的影响进行敏感性分析；并且给出不同空间位置（X与Z）处的浓度变化曲线。其中，采用的参数值如未作特殊说明，基本参数与表 1中一致（即X=Z=0.5，λ=0.25，κ=0.1、η=2.5及T=0.1），且a=0.1。

3.1   浓度函数参数a

不同浓度函数参数a对于污染物浓度随时间变化的影响如图 3所示，其中a取值分别为0，10^-4，10^-3，10^-2，10^-1，1。从图 3中可看出，不同的浓度函数参数引起的浓度变化相差较大。具体地，当a=0时，浓度在变为1后保持恒定。当a > 0时，土中的污染物浓度均存在峰值，且过了峰值后逐渐衰减为0（当a=10^-4，10^-3时浓度也出现峰值，且最终衰减为0，考虑对比效果未取至0的时间点）；当浓度参数a逐渐增大时，浓度的峰值逐渐减小，且出现峰值的时间以及最终浓度衰减为0的时间均缩短。

图  3  不同浓度函数参数a对污染物浓度变化的影响

Figure  3.  Influences of parameter of concentration function a on variations of contaminant concentration

下载: 全尺寸图片 幻灯片

出现以上规律是因为污染源处的污染物浓度随时间进行指数衰减，且当污染物不断向土层中迁移时污染源处的浓度最终会趋于0，则土体中的污染物浓度终将由于对流作用的存在而变为0。另外，当浓度函数参数a逐渐增大时，污染源处的浓度衰减越快，则土体中污染物浓度达到的峰值也越小。

3.2   土工膜破损率$\lambda$

图 4分析了不同的土工膜破损率$\lambda$对污染物浓度变化的影响，采用$\lambda$为0，1%，5%，25%，50%，100%对浓度变化进行研究。其中$\lambda$=0表示土工膜完整，$\lambda$=100%表示土层上表面无土工膜。

图  4  不同土工膜破损率λ对污染物浓度变化的影响

Figure  4.  Influences of defect rate λ on variations of contaminant concentration

下载: 全尺寸图片 幻灯片

从图 4可知，破损率$\lambda$越大，则浓度曲线的峰值越大且峰值出现的时间越短。这是由于破损处的边界为浓度边界，破洞越大则污染源的影响范围越大，则污染物在衰减的同时进入土层中的污染物量越多，表现为土中相同位置的浓度峰值越大且峰值出现时间越短。另外，引入∆C_{p, λ}/∆$\lambda$（C_{p, λ}是破损率为$\lambda$时的浓度峰值）来分析破损的出现对浓度峰值出现的敏感性影响。(C_{p, λ=1%}-C_{p, λ=0})/(λ_1%-λ₀)≈440，(C_{p, λ=100%}-C_{p, λ=1%})/ (λ_100%-λ_1%)≈0.3，说明土工膜的破损较小（$\lambda$ < 1%）时对于浓度变化影响的敏感性更高，也可说明即使土工膜破损很小（如$\lambda$ < 1%）也可极大地提高污染物的迁移效率，这与文献中的结论一致。

3.3   各向异性参数κ

图 5描述了各向异性扩散系数κ对于污染物浓度变化的影响情况，其中采用κ为10^-2，10^-1，1，10（${D_{{{\text{s}}_{\text{z}}}}}$保持不变）进行研究。与图 4相似的是，随着κ数值的变大，污染物浓度峰值越大且出现时间越早。但是不同的是该现象是由于κ的变大意味着水平向的扩散系数${D_{{\text{sx}}}}$也变大（${D_{{\text{sz}}}}$保持恒定），因此在相同条件下污染物更容易往水平方向扩散，表现为浓度峰值的增大和峰值出现时间缩短。

图  5  各向异性参数κ对污染物浓度变化的影响

Figure  5.  Influences of anisotropic coefficient κ on variations of contaminant concentration

下载: 全尺寸图片 幻灯片

3.4   模型厚宽比η

图 6所示为污染物在具有不同尺寸的模型中迁移的曲线图，且采用模型的不同厚宽比η（其中，土层厚度L_s与破损率$\lambda$均保持恒定）进行研究，模型参数取值η为1，2，4，8。

图  6  模型厚宽比η对污染物浓度变化的影响

Figure  6.  Influences of thickness-width ratio η on variations of contaminant concentration

下载: 全尺寸图片 幻灯片

由图 6可知，污染物浓度峰值随η变大而增大，且峰值出现时间越短。这是由于η越大时，模型的宽度越大，但是破损率保持恒定，则土工膜破损的尺寸越大。结合图 4中的结果，图 6中的结果也从侧面说明了土工膜破损尺寸的变化比相对水平位置变化的影响更大。

3.5   不同相对空间位置处的浓度分布曲线

图 7展示了T=0.1时污染物浓度在水平方向的分布情况。其中，图 7为不同相对深度（Z为0，0.25，0.5，0.75，1）时污染物浓度随着水平方向分布的曲线。当Z=0时，在0 < X < 0.25（即破损处）范围内的浓度始终保持为1，0.25＜X＜1内的浓度从C_T=0.1（=e^-0.01≈ 0.99）逐渐下降。在0＜Z＜1内，浓度曲线随着深度位置增大而下移，浓度随水平位置的增加而越小。

图  7  不同深度位置污染物浓度随水平方向的分布

Figure  7.  Distribution profiles along horizontal direction at different depths

下载: 全尺寸图片 幻灯片

图 8为T=0.1时在不同水平位置X处沿着相对深度Z的浓度分布曲线，可看出浓度随着深度的变大整体变小。具体来说，水平位置越靠近破损处的土体中污染物的浓度越高。另外，当X为0，0.25时，土中顶部的污染物浓度为C_T=0.1，为破损处的浓度，其他水平位置（X为0.5，0.75，1）时的顶部浓度均比X=0与X=0.25处的浓度小。与一维模型^{[3, 16-17]}的结果相比，而本研究中污染物浓度随空间分布更符合实际情况。

图  8  不同水平位置处污染物浓度随深度的分布

Figure  8.  Distribution profiles along longitudinal direction at different horizontal positions

下载: 全尺寸图片 幻灯片

4.   结论

本文考虑了污染源浓度随时间变化的情况，建立了填埋场中污染物通过破损的土工膜并迁移到下层土体的二维计算模型。通过Laplace变换、Fourier变换及边界转换法得到了半解析解，并采用数值解验证了该半解析解的准确性。针对模型中浓度函数参数a、破损率$\lambda$、各向异性参数κ、厚宽比η进行了敏感性分析，并得到污染物浓度随深度分布曲线。

（1）不同深度处污染物浓度沿着水平方向逐渐减小，且在土层顶面浓度与破损处的浓度边界相同，比一维迁移模型更符合实际情况。

（2）污染源浓度按指数衰减时，土中污染物浓度随时间变化时出现峰值。与各向异性参数κ与厚宽比η对土中的污染物浓度影响规律相似，浓度函数参数a越大时，峰值越大，且出现的时间越短。

（3）土工膜中的破损可大幅加快污染物在土层中的迁移，且破损率$\lambda$越大，污染物的迁移越快。