引言
我国渤海湾、黄河三角洲以及江苏东部的沿海地区分布着广阔的滨海盐渍土，据统计其面积为2 万多 km2，如何对这些后备的土地资源进行科学使用与管理，一直是土壤工作者关注的问题。由于滨海盐土具有含盐量高、养分低、地下水位埋藏浅、矿化度高、土质粘重等特点，需要经过改良后才能利用。改良工程中，如何使质地黏重、渗透性差的盐土快速脱盐是生产实践中最重要的关键问题。打孔注入砂性土柱(简称砂柱) 工艺技术是目前工程改良实践中一种很有效的技术，可以解决黏质低渗性盐土的难渗透问题，达到暗管排水达不到的渗透效果。砂柱对土壤的改良应用早在 1981年便有相关实验研究，该技术目前已经在中新天津生态城、江苏启东工业园等地进行了实践应用。但目前的工程技术都是采用内设竖砂柱的方式，对于其他砂柱设置方式的效果比较和内在机理，尚缺乏相关的深入研究和定量分析。
20 世纪 80 年代后，许多学者基于室内土柱实验，应用数值模拟的方法研究了土壤中水分和盐分一维运动的规律，为本研究提供了参考。但对于 2-D 土壤水盐土柱物理模拟和数值模拟的研究，由于其问题的复杂性，研究报道较少，左强曾对2-D 均质土壤排水条件下饱和非饱和水盐运动规律进行过初步数值模拟研究，采用了有限单元法迭代求解。对于江苏一带的砂质滨海盐土，张亚年曾用 HYRDUS 软件进行过暗管排水条件下水盐运移的数值模拟。
COMSOL Multiphysics 软件 (原 Finite ElementModeling Laboratory) 是一个基于偏微分方程的多物理场有限元分析软件，可以用来求解线性、非线性问题，和时间有关的稳态、瞬态问题，以及和几何形状有关的一维、二维和三维问题。与专门针对土壤水盐模拟的软件相比，COMSOL 在处理实际问题和数值计算上适用性更广，可模拟一些更复杂的工程设置方式。目前运用 COMSOL 软件于土壤中的计算和模拟，主要局限于地下水和土壤水入渗的研究，Wissmeier 将该软件应用于土壤中杀虫剂的运移模拟，对运用 COMSOL 模拟土壤溶质运移的适用性进行了验证，证明了该软件适于土壤溶质运移的模拟。
本文将土体中建构砂柱后的水盐运动简化成非均质 2-D 土壤水盐运动，应用 COMSOL 软件模拟工程实践可实现的几种砂柱构型，定量探究其加速脱盐效果，以期为滨海盐土脱盐、控盐的工程模式进行定量评价和设计优化提供指导和依据。
1、 材料与方法
1. 1 土样来源及性质
盐土柱用土取自天津大港的粘壤质滨海盐土，考虑到砂柱用砂子时，会因两个土壤导水率相差过大而引起指流或达不到实践上的脱盐效果，所以本次模拟砂柱用土选为砂壤土。两种土壤的颗粒组成(按美国制分类) 见表 1。
1. 2 试验设计
本文在粘壤质滨海盐土土柱中设置了 3 种典型形式的砂柱，与空白盐土土柱进行对比，对它们的洗盐过程进行了模拟，砂柱设置方式如图 1，图中，(a)为没有砂柱的土柱，(b) 为中竖砂柱，(c) 为单斜对角砂柱，(d) 为双斜“X”型长砂柱，左边的三维图为中竖沙柱的三维视图，蓝色部分为粘壤质盐土，黄色部分为砂壤质土柱(砂柱) 。
1. 3 COMSOL 的适用性验证与模型原理
本研究运用 COMSOL 2-D 地球科学模块的流体流动模块和溶质运移模块，模拟的水分运动采用饱和-非饱和水流的 Richards 方程，溶质运移基于 Fick定律的对流-弥散方程，将两个模式耦合在同一物理场中，来模拟淋洗条件下土壤水分和盐分(Cl－) 的二维运移过程。
为了确保 COMSOL 2-D 的地球科学模块流体流动模块和溶质运移模块耦合可以用于对土壤水盐运动的模拟，首先对 Warrick 等经典的田间咸水灌溉和淡水冲洗试验、与 Frind 所列举的 2-D 地下水污染物运动问题进行了模拟计算和比较，结果表明: 在合理的网格剖分下，应用 COMSOL 所得的结果与已有的实测值、解析解一致。由于篇幅所限，这里不再详述。
结合工程实践和便于模拟计算，模拟对象统一选为高1 m、水平断面为边长0. 3 m 正方形的长方体盐土土体，内设砂柱的水平断面为边长 0. 04 m 正方形，如图 1 所示。本模拟中所有土柱经过淹水，土壤饱和后，表层保持 2 cm 积水进行淋洗，水分在土壤中入渗的过程属于空间三维运动，由于 3-D 问题极其复杂，将问题简化为非均质、各向同性的水盐运动二维问题来模拟，此时，土壤水分运动控制方程可表达为
SS。根据文献，土柱中粘壤质盐土的 SS取 5 ×10－ 5m－ 1，砂柱砂壤土的 SS取1.3 ×10－ 5m－ 1。
溶质运移的控制方程可表达为
根据表 1 土壤容重及颗粒组成，采用 ROSETTA软件拟合了两种土壤的水力学参数，结果见表 2，其中 Ks为土壤饱和导水率。
1. 5 定解条件的确定
在整个模拟过程中，模拟对象均为高 1 m、水平断面为边长 0. 3 m 正方形的长方体盐土土体，土柱填充土壤为滨海粘壤质盐土，方砂柱的土壤为砂壤土，方砂柱的水平断面为边长 0. 04 m 正方形，土柱下端接地下水通暗管排水，地下水质量浓度为0. 1 g / L。本模拟中所有土柱先淹水饱和后再淋洗盐分，模拟初始含水率均为饱和含水率，表层保持0. 02 m积水进行淡水淋洗。
上述定解条件的数学表达式如下:
(1) 水分运动定解问题粘壤质盐土土体部分，简称 C:
水流方程初始条件为 hC(x，z，0) =0 m。 上边界条件为hC(x，0，t) =0. 02 m。 下边界条件为 hC(x，－1，t) =0 m。 内设砂柱部分，简称 S: 水流方程初始条件为 hS(x，z，0) = 0 m。
上下边界条件由砂柱位置决定，砂柱上下边界位置与土柱上下边界位置相同时，上边界条件为hS(x，0，t) =0. 02 m。下边界条件为 hS(x，－1，t) =0 m。砂柱上下边界位置低于或高于土柱上下边界位置时，砂柱上下边界设为内部连续边界。
(2) 溶质(Cl－) 运移定解问题初始条件为 c(x，z，0) =8. 389 kg/m3。
上边界为盐通量 J =0 kg/(m2·s) 。 下边界为 c(x，－1，t) =0 kg/m3。
1. 6 网格剖分COMSOL Multiphysics 的网格剖分功能可以创建自由网格、边界层网格等。利用其网格剖分工具和方法，可以生成三角形和四边形(2D) 网格，自动建立一致性的网格。当模拟的对象是由不同材料组成时，边界的一致性网格对得到精确解是非常重要的，COMSOL 的网格剖分可以使通过界面的解分量及其通量连续。4 种设置方式的网格模式详见 1. 2 节中图 1，在粗和细两种质地土壤边界处的网格密度均加密到 0. 005 m，因此保证了整个模拟的可靠性，4 种设置方式的网格局部放大图见图 2。
4 种设置情况下的网格剖分具体参数详见表 3。
2、 结果及分析
图 3、图 4 和图 5 比较了 4 种设置情况下，经淋洗 12 h、24 h 和36 h 后，土柱中的 Cl－离子浓度分布状况。表 4 列出了 4 种配置条件下不同时间的盐分淋出量及百分比。
从图 3 可以看出: 经过 12 h 淋洗后，不设砂柱、内设中竖砂柱、内设单斜对角砂柱、内设双斜“X”型砂柱的土柱盐分淋出量分别为各自盐土含盐总量的22% 、29% 、29% 、33% ，内设中竖砂柱、内设单斜对角砂柱、内设双斜“X”型砂柱的土柱淋洗效率分别为不设砂柱的 1. 32、1. 32 和 1. 50 倍。从图 4 可以看出: 经过 24 h 淋洗后，各土柱盐分淋出量分别为各自盐土含盐总量的 42%、52%、54% 和 61%，内设中竖砂柱、内设单斜对角砂柱、内设双斜“X”型砂柱的土柱盐分淋洗效率分别为不设砂柱的 1. 24、1. 29和 1. 45 倍。从图 5 可以看出: 经过 36 h 淋洗后，各土柱的盐分淋出量分别为各自盐土含盐总量的61% 、73% 、77% 和 84% ，内设中竖砂柱、内设单斜对角砂柱、内设双斜“X”型砂柱的土柱的淋洗效率分别为不设砂柱的 1. 20、1. 26 和 1. 38 倍。
因此，综合以上数据，在内设中竖砂柱、内设单斜对角砂柱、内设双斜“X”型砂柱 3 种砂柱设置方式中，双斜“X”型长砂柱的构型设置最有利于盐分的淋洗。
图 6 为各土柱随淋洗时间的变化，底层截面的盐分累计通量与各自盐土初始盐分总量的比值变化图，可以定量地判断各设置中盐分淋出量随时间的变化。图 1 所示各土柱剖面中，各自盐土的初始总量分别为 8. 389、7. 27、7. 27 和 6. 152 kg/m2，由图 6可以看出: 同时间下内设双斜“X”型砂柱的土柱盐分淋洗效率是最高的，内设单斜对角砂柱的土柱次之，内设中竖砂柱的土柱第三，三者盐分淋洗效率均高于不设砂柱的土柱。
图 7 为各个土柱随淋洗时间的变化，二维剖面盐分淋出比与水通量的比率变化图。盐土中水的流速为 2. 362 × 10－ 6m / s，纯砂土中水的流速为5. 376 × 10－ 6m / s，各个土柱的二维剖面中水的通量分别 为 6. 96 × 10－ 7、8. 14 × 10－ 7、8. 14 × 10－ 7、9. 32 ×10－ 7m2/ s。由图 7 可衡量出，60 h 之前，内设双斜“X”型砂柱的方式的水资源利用效率是最佳的。
表 4 为每种情形下，当淋出的盐分分别达到0. 50、0. 75、0. 95 时所用时间。各土柱的 50% 盐分淋出时长分别为 31、23、22、19 h，内设中竖砂柱、内设单斜对角砂柱、内设双斜“X”型砂柱比不设砂柱分别节省时间 8、9、12 h，淋出比与水通量的比率分别为 6. 652、7. 418、7. 756、7. 843 m－ 2; 各 土 柱 的75% 盐分的淋出时长分别为 45、37、35、31 h; 内设中竖砂柱、内设单斜对角砂柱、内设双斜“X”型砂柱比不设砂柱分别节省时间 8、10、14 h，淋出比与水通量的比率分别为 6. 437、6. 917、7. 313、7. 211 m－ 2; 各土柱的 95% 盐分的淋出时长分别为 59、52、49、46 h，内设中竖砂柱、内设单斜对角砂柱、内设双斜“X”型砂柱比不设砂柱分别节省 7、10、13 h，淋出比与水通量的比率分别为 6. 426、6. 234、6. 616、6. 155 m－ 2。
综合以上数据结果，表明: 双斜“X”型长砂柱和单斜砂柱对盐分的淋洗效果提高较明显，水资源利用效率也更佳。原因如下:根据土壤溶质运移规律，氯离子的运动是对流和水动力弥散(包括扩散和机械弥散) 物理过程综合作用的结果。对流的主要影响因素是水流的流速，从表2可以看出，砂柱的纵向导水率45. 54 cm / d(5. 27 × 10－ 6m / s) 明显高于盐土的导水率 20. 01 cm/d(2. 32 ×10－ 6m / s) ，从 COMSOL 模拟结果可以看出，纯盐土中水的流速低于纯砂土中水的流速，所以砂柱中氯离子的对流明显高于盐土中氯离子的对流。
通过对水势场和相应流速场的分析，在中竖砂柱设置中，同一水平线上的砂柱和盐土中的水压相同，因此斜向压力分量为零，土柱水流和砂柱水流均为纵向，砂柱和土柱之间水分交换近似为零。在单斜对角砂柱和双斜“X”型长砂柱中，砂柱中水分流动较快，同水平线上砂柱水压低于盐土的水压，因此在边界处水不断从盐土流向砂柱。
影响氯离子运移的另外一个因素是盐分浓度差和水动力弥散，水动力弥散是土壤水微观流速变化而引起的溶质的弥散，从溶质运移控制方程(2) 可以看出，溶质的弥散强度与曲折因子 τL、纵向弥散度 αT、横向弥散度 αL呈正相关关系，表 2 中砂土的τL值、αT值、αL值均高于盐土的值，所以砂柱中氯离子的水动力弥散强度明显高于盐土中氯离子的水动力弥散强度。从模拟过程中的速度场和盐分浓度场动态变化分析可知，单斜对角砂柱和双斜“X”型长砂柱的设置与中竖砂柱设置相比，会有更多的盐分从盐土中向砂柱中运动。
综合盐分的对流和水动力弥散作用，再加上砂柱的形态设置和大小的不同，就会导致图 6 所示的结果:
(1) 经过 12 h 淋洗后，不设砂柱、内设中竖砂柱、单斜对角砂柱和双斜“X”型长砂柱的盐分淋出量分别为各自盐土含盐总量的 22%、29%、29%、33% ，后 3 种砂柱设置方式的 12 h 盐分淋洗效率分别为不设砂柱的1. 32、1. 32 和1. 50 倍; 经过24 h 淋洗后，4 种设置的盐分淋出量分别为各自盐土含盐总量的 42%、52%、54%和 61%，后 3 种砂柱设置方式的 24 h 淋洗效率分别为不设砂柱的 1. 24、1. 29和 1. 45 倍; 经过 36 h 淋洗后，4 种设置的盐分淋出量分别为各自盐土含盐总量的 61%、73%、77% 和84% ，后 3 种砂柱设置方式的 36 h 淋洗效率分别为不设砂柱的 1. 20、1. 26 和 1. 38 倍。
(2) 不设砂柱、内设中竖砂柱、单斜对角砂柱和双斜“X”型长砂柱4 种设置的50%盐分淋出时长分别为31、23、22 和19 h，后3 种砂柱设置比不设砂柱分别节省时间 8、9、12 h; 4 种设置的 75% 盐分的淋出时长分别为 45、37、35、31 h，后 3 种砂柱设置比不设砂柱分别节省淋洗时间 8、10、14 h; 4 种设置的95% 盐分的淋出时长分别为 59、52、49、46 h，后 3 种砂柱设置比不设砂柱分别节省淋洗时间 7、10、13 h。
3、 结论
(1) 盐土内设砂柱可促进盐土的水盐运移，有利于提高盐分的淋洗效率，其效果受砂柱的角度和位置影响，目前中竖砂柱、单斜对角砂柱和“X”型长砂柱 3 种砂柱设置方式中，双斜“X”型长砂柱的构型最有利于盐分的淋洗，单斜对角砂柱次之，中竖砂柱的效率比不设砂柱稍好，但不及单斜对角砂柱和双斜“X”型长砂柱。
(2) 本模拟条件下，经过 24 h 淋洗后，不设砂柱、内设中竖砂柱、单斜对角砂柱和双斜“X”型长砂柱 4 种设置的盐分淋出量分别为各自盐土含盐总量的 42%、52%、54%和 61%; 4 种处理下 75% 盐分的淋出时长分别为 45、37、35、31 h。
(3) 结合盐分淋洗效率、相同盐分淋出比下的水资源利用效率两种因素来综合考虑，双斜“X”型砂柱的设置方式是最佳的设置。
参考文献： 1. 徐攸在． 盐渍土地基［M］． 北京: 中国建筑工业出版社，1993: 1 － 10． 2.王遵亲． 中国盐渍土［M］． 北京: 科学出版社，1993． 3.晓洋，陈效民，李孝良，等． 不同改良剂与石膏配施对滨海盐渍土的改良效果研究［J］． 水土保持通报，2012，32(3) : 128 －132．