基于COMSOL软件的绿洲盐渍化土壤中多离子耦合运移模型构建

焦会青，盛 钰，赵成义，李保国 ※

（1．中国农业大学资源与环境学院，农业部华北耕地保育重点实验室，国土资源部农用地质量与监控重点实验室，北京100193；2．中国科学院新疆生态与地理研究所，乌鲁木齐 830011）

0 引 言

水资源短缺和土壤盐渍化是制约干旱区农业发展的重要因素[1-2]，盐碱土对作物生长的危害，不仅来自于耕层土壤较高的盐分浓度，也来自于盐分离子的组成差异[3]。土壤中离子之间相互影响、相互作用，由于各类盐分离子的化学特性存在差异，使其在土壤中的迁移能力以及淋洗难易程度不同，对作物的影响机制和程度也有很大差异。因此，探求土壤中各盐分离子的迁移规律可以为盐渍化土壤的综合治理以及高效利用提供科学依据[4]。

基于过程的数值模型是研究土壤水和溶质运移的有力工具[5-7]，Robbins等[8]将化学沉淀溶解和离子交换平衡模型子程序与一维非饱和水流耦合，并用蒸渗仪[9]和田间试验[10]数据对模型进行了验证。Chen等[11-12]建立了一维稳态饱和流条件下Ca2+、Na+、SO42-和Cl-的耦合运移模型，以动力学模型描述硫酸钙的沉淀溶解，以Gapon方程描述阳离子的交换过程。Šimůnek和Suarez[13]、Suarez和Šimůnek[14]建立了UNSATCHEM，用于模拟非饱和多孔介质中离子的平衡或者动力学反应过程，模型考虑温度及CO2的影响，并且已经成功用于模拟CaSO4对钠质土的改良[15]、蒸渗仪中滴灌条件下盐分离子的运移[16]等过程，近些年，其中主要的离子反应模块已经嵌入HYDRUS，但是由于缺乏相应的试验数据，而且涉及参数较多，该模块应用相对较少[17-18]。Lekakis和Antonopoulos[19]在一维土壤水氮运移模型WANISIM中加入离子（Ca2+、Na+、Mg2+）运移模块，使模型可用于盐分管理工作，但是模型中仅考虑了阳离子的交换和吸附。

COMSOL是一款高度集成的数值仿真软件，基于有限元理论，以高效的计算性能和杰出的多场耦合分析能力实现任意多物理场的数值仿真。在地球科学领域，COMSOL可以模拟地下水流动、土壤中溶质运移和气流形成规律、地热、地面沉降等问题[20-24]，Rosenbom等[25]采用COMSOL评估了小尺度的空间变异对土壤中杀虫剂降解和淋洗过程的影响；Wang等[26]采用COMSOL模拟了一维土柱中水盐运移以及盐分结晶过程，以此研究盐分结晶对路基的影响程度；Mao等[27]基于COMSOL构建了土壤与地上部的耦合模型，模拟了杀虫剂在土壤中迁移以及挥发至大气的过程。COMSOL中提供了非饱和多孔介质溶质运移基本模块，没有考虑土壤中各组分之间的化学反应，因此，将COMSOL应用于非饱和多孔介质盐分离子耦合运移的研究较少。但是，COMSOL支持自定义偏微分方程（PDE）的求解，用户可以在该模式下自由定义所需方程并实现方程之间的耦合，而且在模型的构建中，材料属性、源汇项以及各种边界条件都可以使用任意类型的函数，比如分段、插值、解析等，这些函数既可以时间相关，也可以空间相关，其自变量既可以是独立变量，也可以是求解变量本身[28]，借助于这些优势可以在COMSOL中进行物理问题的二次开发，而不需要另外编程求解，对于实现盐分离子之间特定的化学反应有很大优势。

本研究采用COMSOL多孔介质和地下水流模块模拟土壤水分运动，以溶质运移及盐分离子反应理论为基础，基于 COMSOL自定义偏微分方程模块构建盐渍化土壤SO42-、Ca2+、Na+、Cl-、Mg2+耦合运移模型，并以绿洲棉田膜下滴灌条件下主要盐分离子的迁移为例对模型进行了验证，为盐分离子运移模拟研究提供了一种简便的模型构建方法。

1 非饱和多孔介质盐分多离子耦合运移数学模型

1.1 控制方程

土壤中水分运动以 Richards方程描述，其中，根系吸水基于潜在蒸腾量、根系分布以及土壤基质势计算，具体方法见文献[29]。

本研究主要探讨 SO42-、Ca2+、Na+、Cl-、Mg2+的互相作用和迁移过程，离子的运移过程基于对流弥散方程构建，SO42-、Ca2+、Na+、Cl-、Mg2+的运移方程如下，

式中：θ为土壤体积含水量，L3/L3；c1-c5分别为土壤溶液中 SO42-、Ca2+、Na+、Cl-、Mg2+的浓度，N/L3；s2、s3、s5分别为Ca2+、Na+、Mg2+的吸附量，N/M；t为时间，T；xi为空间坐标，L；ui为土壤孔隙流速，L/T；ρb为土壤容重，M/L3；X为CaSO4的沉淀量，N/L3；Dij为水动力弥散系数，L2/T，计算方法参考文献[30]。

1.2 离子反应

对于盐渍化土壤，土壤胶体上供吸附阳离子的剩余电荷一般较少，因此吸附作用通常较小，主要发生的是阳离子的交换过程。如果土壤溶液中离子组成或浓度发生改变，势必会影响各离子在溶液和胶体之间的平衡状态。假设土壤中的阳离子交换量不变且等于交换性Ca2+、Mg2+、Na+之和，用 Gapon方程描述Ca2+、Na+、Mg2+之间的交换过程[29]，

式中：KNa/Ca为Na+和Ca2+的交换常数，(N/L3)-1/2；KMg/Ca为 Mg2+和 Ca2+的交换常数；γ2、γ3、γ5分别为 Ca2+、Na+、Mg2+的活度系数。

硫酸钙的沉淀溶解满足溶度积原理，反应以二级动力学过程描述[31]，

式中：Ksp为溶度积常数，(N/L3)2；kd为溶解常数，L3/(N·T)；kp为沉淀常数，L3/(N·T)；α(X)为单位跃迁函数，即 X ＞ 0时为1，X ≤ 0时为0。

1.3 活度系数

土壤溶液中可溶盐浓度增加或减少会改变各离子的活度系数，进而影响阳离子交换和溶解沉淀过程，活度系数的估算可采用Davies方程[31]，但是当离子强度I ＞ 0.1 mol/L时，活度系数与离子强度不再是通用的函数关系，而是取决于特定离子的交互作用[32]，因此，这里也可以采用适用于特定条件下的函数方程。

2 基于COMSOL盐分多离子耦合运移模型实现方法

以COMSOL多孔介质和地下水流模块模拟土壤水分运动，而对于COMSOL没有定义的物理问题，可以通过自定义偏微分方程组（PDEs）实现二次开发，COMSOL中提供的系数型PDE如下，

式中：u 为求解变量；ea、da、c、α、γ、β、a、f均为自定义系数。

所有自定义系数可根据实际物理问题定义为不同类型的函数，这些函数可以是不同阶次的导数，既可以时间相关，也可以空间相关，其自变量既可以是独立变量，也可以是求解变量本身。对于各离子的控制方程，首先根据公式（9）的形式进行整理，然后在相应的编辑框中定义各系数，对于没有对应项的部分则可通过定义系数f实现，而不再需要编程求解。就离子运移方程而言，可设定为对流弥散项，f则可作为源汇项，表征离子间的化学反应。各离子运移方程通过土壤含水量和达西流速与 Richards方程耦合，又通过离子之间的化学反应实现各离子方程之间的耦合。

相应的，该模块也提供了多种开放的边界条件，比如Dirichlet边界和通量边界条件，通量边界定义如下可构成方程的时变项，

式中，g和q均为自定义系数。

3 模型检验：绿洲棉田膜下滴灌土壤盐分离子运移

以新疆绿洲膜下滴灌条件下的离子迁移过程为例，通过计算模拟值与实测值之间的平均绝对误差（MAE）、均方根误差（RMSE）、平均相对误差（MRE）以及决定系数（R2）对模型进行检验。

3.1 试验及模型设置

试验在新疆阿克苏水平衡试验站开展，共设置 3个处理，单次灌水量分别为46.8 mm（T1）、39.6 mm（T2）、28.8 mm (T3)，在无降雨的情况下每隔6 d灌溉一次。棉花的种植模式如图1所示，窄行10 cm，宽行65 cm，膜间距60 cm，一膜覆盖四行棉花，滴灌管铺设在宽行中间，滴头间距30 cm。2015年和2016年膜下滴灌期间，在宽行和膜间的中间位置采集土壤样品，采样深度分别为0～20、20～40、40～60、60～80、80～100 cm，采用烘干法测定土壤含水量，每个灌水周期测定两次。2015年8月1日至2015年8月25日在宽行和膜间深20、40、60、80 cm处埋设土壤溶液提取器，定点采集土壤溶液，共采集 6次。土壤溶液提取器主要由取样器、负压泵、采样瓶组成，所用取样器陶土头的直径为2 cm、长度为5 cm，当外加吸力大于土壤基质吸力时，使土壤溶液通过陶土头进入采样瓶，将收集的土壤溶液用于盐分离子的测定。其中Cl-采用AgNO3滴定法测定，SO42-采用EDTA间接滴定法测定，Ca2+和Mg2+采用EDTA络合滴定法测定，Na+和K+采用火焰光度计测定。灌溉水矿化度0.485 g/L，SO42-、Ca2+、Na+、Cl-、Mg2+浓度分别为 1.20、0.88、4.09、2.96、0.46 mmol/L。

模拟区域宽75 cm，高100 cm，模拟区域及水分运动的边界条件见图1，对于盐分运动，考虑灌溉带入的盐分，但是由于降水量非常小，所以忽略降水带入的盐分，所有边界均采用第三类边界条件。将灌溉前的实测含水量和各离子含量进行二维插值，作为模拟初始值。

图1 棉花种植模式及水分运动边界条件示意图Fig.1 s of planting mode of cotton and boundary conditions for soil water flow

土壤颗粒组成见表1，并根据Rosetta得到土壤水力学参数（van Genuchten方程[33]）的初始值，采用 2015年的土壤含水量数据对模型进行校准，并采用2016年的试验数据验证模型，率定的参数值见表1。计算弥散系数时所采用的纵向弥散度取12 cm，横向弥散度为纵向弥散度的1/10，溶质在纯水中的扩散系数取1.296 cm2/d，根系吸水的水分胁迫函数（Feddes模型）中h1、h2、h3、h4分别取-10、-25、-500、-14 000 cm[34]。

阳离子交换量CEC、Gapon常数KNa/Ca和KMg/Ca参考文献[17，35]分别设定为 6.996 × 10-2mol/kg、0.50(mol/L)-1/2、0.52；硫酸钙溶度积常数Ksp取2.399 × 10-5(mol/L)2[31]，硫酸钙溶解速率与土壤水流速有关[36], kd取3.13 × 10-3L/(mol·s)，沉淀速率常数 kp取 2.22 × 10-3L/(mol·s)[35]。由于试验地抽提的土壤溶液样品中 90%的离子强度都大于0.1 mol/L，且土壤水溶性SO42−占全盐比例约为 67%，Ca2+占全盐比例约为 21%，假设土壤溶液中Ca2+和SO42−处于平衡状态，根据CaSO4的溶度积常数和土壤溶液中 SO42−和 Ca2+的浓度反推出 SO42−和 Ca2+的活度系数，以离子强度为横坐标，活度系数为纵坐标，绘制各个散点，并拟合得到指数方程，与Davies方程的对比结果见图2。

表1 土壤颗粒组成及水力学参数Table 1 Particle size fractions and corresponding soil hydraulic parameter values

图2 活度系数与离子强度的关系Fig.2 Relationship between activity coefficient and ionic strength

3.2 模拟结果验证

表2给出了2015年和2016年各处理土壤含水量模拟值与实测值的统计结果，MAE值介于 0.023～0.033 cm3/cm3，RMSE 值介于0.030～0.040 cm3/cm3，与Phogat等[37]和Ramos等[18]得到的结果相似。

表2 土壤含水量模拟值与实测值的统计对比Table 2 Statistical comparison of simulated and measured soil water contents

2015年各离子模拟值与实测值统计对比结果见表3，Ca2+和SO42-的模拟效果总体上优于Na+、Cl-和Mg2+，处理一和处理二优于处理三，这是因为处理三灌溉量较小，只有少量的灌溉水能够到达膜间位置，导致膜间始终保持在相对干燥的状态，在整个观测期内膜间深度 20 cm处仅收集到一次土壤溶液，而且为了收集到能够满足测量的水量，抽提时间相对较长，因此在膜间没有合理的初始值导致了该处理模拟值与实测值较大的差异。通常情况下，RMSE ≥ MAE，而RMSE大于MAE的程度往往可以反映实测值与模拟值之间差异的变化程度，或者说在一定意义上可以反映最大差异的程度[18]，正因为处理三膜间位置较大的模拟误差，导致处理三RMSE明显大于MAE。该模型模拟值与实测值之间的差异，一方面可以归因于土壤溶液离子浓度的测定方法，土壤溶液抽提器采集的是陶土头周围一定时间和空间范围内的土壤溶液，该范围与陶土头周围土壤的水力学性质和基质势有关[38]，但是模型给出的是陶土头相应位置点的数据；另一方面可归因于模型参数的设置，阳离子交换和硫酸钙沉淀溶解反应的相关参数均根据文献给出，而这些参数与特定的试验环境有关，也必然会导致模拟值与实测值的差异。独立扰动各反应参数，当阳离子交换量 CEC、Na+和 Ca2+的交换常数 KNa/Ca、Mg2+和 Ca2+的交换常数KMg/Ca的变化幅度为±20%时，与本研究模拟结果相比，所有离子模拟值与实测值之间的 RMSE值分别波动-2.60%～3.99%、0.25%～3.45%、-1.45%～2.74%；当硫酸钙沉淀常数 kp、硫酸钙溶解常数 kd的变化幅度为±20%时，所有离子模拟值与实测值之间的RMSE值分别波动-0.29%～1.51%、-9.36%～11.10%，其中，kd对SO42-和Ca2+的模拟精度影响较大，当 kd的变化幅度为±10%时，其RMSE值波动范围介于-4.69%～7.36%，这是因为灌溉使膜下土壤溶液中SO42-和Ca2+浓度明显降低，打破了硫酸钙的沉淀溶解平衡，主要发生着硫酸钙的溶解过程，kd影响硫酸钙的溶解速度，因此会对SO42-和Ca2+的模拟产生较大的影响。另外，气温对上部土层水溶性离子含量的变化有明显的影响[39]，本研究没有考虑温度因素也可能是产生模拟误差的原因之一。

3.3 膜下滴灌条件下各离子在土壤剖面上的运移特征

表3 土壤溶液离子浓度模拟值与实测值的统计对比Table 3 Statistical comparison of simulated and observed ion concentrations in soil solution

以处理二为例，图3给出了一个灌水周期（7 d）内土壤含水量的变化过程，从第二天开始灌溉，共 24 h。滴灌结束后，滴灌管下方土壤含水量接近饱和，但是由于膜下滴灌条件下，仅有少量的灌溉水可以到达膜间位置，使膜间地表在灌溉后也保持在相对干燥的状态，随着时间的推移，滴灌头周围的土壤水分逐渐向深层和膜间运动，使膜下土壤含水量逐渐降低。由于土壤颗粒组成的差异以及地下水的补给等原因，32 cm以下土层的含水量总体上高于上层。

图3 单个灌水周期内土壤含水量（cm3·cm-3）在剖面上的变化过程Fig.3 Change of soil water contents in profile during an irrigation cycle

图 4给出了各离子在一个灌水周期内的变化过程，灌水结束后所有盐分离子在膜下位置都明显降低，主要是 0～40 cm土层，随着时间的推移，该位置的离子浓度均有不同程度的提升。灌溉导致土壤溶液中Ca2+和SO42-浓度降低，打破硫酸钙的沉淀溶解平衡，促使硫酸钙溶解，直至土壤溶液中Ca2+和SO42-活度的乘积等于溶度积常数；Ca2+浓度的逐渐增加提高了其交换能力，因而进一步促进胶体上吸附的 Na+和 Mg2+释放至土壤溶液；除此之外，根系的吸水过程也会导致离子浓度的升高，这也是 Cl-浓度升高的主要原因。Cl-、Na+在灌溉后浓度升高的幅度小于Ca2+和SO42-，淋洗效果较好，这与文献[40,41]研究结果相似，由此看来，在该试验条件下，滴灌不仅能降低根区土壤盐分浓度，也能明显降低对作物危害较大的Cl-和Na+的组成比例，为棉花的正常生长提供适宜的土壤环境。

膜间地表由于蒸发作用的影响，使各离子均在地表逐渐积累，从图4可以看到，Cl-在膜间地表积累量最大，这主要是由于Cl-不参与化学反应，受土壤固相影响相对较小，移动性较强，而Ca2+和SO42-的活度乘积在大于溶度积常数时就会产生沉淀，使其维持在相对平衡的状态。

3.4 活度系数对离子运移的影响

本研究采用指数方程（由试验数据拟合得到）计算活度系数，为了分析活度系数对离子运移的影响，与基于Davies活度系数方程的模拟结果进行了对比，图5给出了一个灌水周期内2种方法计算得到的活度系数分布。这两种方法都是以离子强度为自变量，且随着离子强度的增加而减小，因此整个剖面上活度系数的变化规律与离子浓度相反，即灌溉后膜下活度系数增加，随着时间的推移又逐渐减小，而在膜间近地表随着蒸发积盐的进行有小幅度的减小。同时，也可以看到在整个土壤剖面上，指数方程计算得到的活度系数小于Davies方程，例如在灌溉后，在整个剖面上，指数方程计算得到的活度系数介于0.07～0.37，而Davies方程计算得到的活度系数介于 0.29～0.57。

图4 单个灌水周期内土壤溶液离子浓度在剖面上的变化过程Fig.4 Changes of various ion concentrations in soil solution in profile during an irrigation cycle

图6为基于Davies活度系数估算方法得到的各离子的模拟结果，与基于指数方程得到的模拟结果（图4）相比，Ca2+和 SO42-的差异最明显。在第一天尚未灌溉时，整个剖面上土壤溶液中Ca2+和SO42-浓度小于图4中模拟结果，这是因为Davies方程计算的活度系数大于指数方程，使土壤溶液中Ca2+和SO42-活度乘积大于溶度积常数，促进了 CaSO4的沉淀过程。在第二天灌溉后，膜下浅层土壤中Ca2+和SO42-都降低，此时CaSO4开始溶解，其溶解速度与溶度积常数和这两种离子活度乘积的差异成正比，因此活度系数较高时 CaSO4的溶解速度较低，导致在灌溉以后膜下Ca2+和SO42-的浓度依然低于图4中的模拟结果。

图5 单个灌水周期内活度系数在剖面上的变化过程Fig.5 Changes of activity coefficient in the profile during an irrigation cycle

两种活度系数计算方法在膜间地表产生的差异则更为明显，基于指数方程得到的Ca2+和SO42-浓度在逐渐升高，而基于Davies方程得到的Ca2+和SO42-浓度有小幅度的减小，这是因为Davies方程计算得到的活度系数较大，模拟开始时刻该处Ca2+和SO42-的活度乘积大于溶度积常数，所以在整个灌水周期内都处于沉淀过程，使 Ca2+和SO42-浓度逐渐降低；而指数方程计算的活度系数相对较小，两种离子活度的乘积尚未达到溶度积常数，因此使离子在此聚积而尚未沉淀。

Na+和Mg2+的变化趋势与图4结果类似，但是总体上略小于图4，这是因为Na+、Mg2+和Ca2+处于动态的交换过程中，当土壤溶液中 Ca2+浓度较小时，其交换能力较小，土壤溶液中的Na+和Mg2+则更易被土壤胶体吸附。基于Davies活度系数估算方法得到的Ca2+、SO42-、Na+、Mg2+的模拟结果与实测值间的 RMSE值分别为 7.87、16.81、6.55、4.93 mmol/L，明显大于基于指数方程得到的结果。因此，活度系数的准确估算对模拟结果的准确性有重要影响，尤其是参与沉淀溶解反应的Ca2+和SO42-，当离子强度I ＞ 0.1mol/L时，活度系数与离子强度不再是通用的函数关系[32]，采用通用的函数关系可能会带来较大的模拟误差。

图6 基于Davies方程单个灌水周期内土壤溶液离子浓度在剖面上的变化Fig. 6 Changes of ion concentrations of soil solution in profile during an irrigation cycle based on Davies Equation

4 结论

1）COMSOL在材料属性、源汇项以及边界条件的设定上非常灵活，而且支持自定义偏微分方程的求解，基于这些优势，在多孔介质和地下水流模块模拟非饱和土壤水流的基础上，自定义偏微分方程组构建盐渍化土壤SO42-、Ca2+、Na+、Cl-、Mg2+耦合运移模型，考虑阳离子交换过程以及硫酸钙的沉淀溶解反应，并以绿洲棉田膜下滴灌田间试验对模型进行了验证，结果表明该模型能够较好地反映土壤中盐分离子的动态变化规律。但是，该模型未考虑碳酸钙的沉淀溶解过程，因此，对于碳酸盐质量分数较高的土壤，该模型不适用，可以进一步在COMSOL自定义偏微分方程模块下定义该过程。

2）膜下滴灌条件下，膜下土壤中的盐分离子有不同程度的淋洗，其中，Cl-、Na+在灌溉后浓度升高速度小于Ca2+和 SO42-，淋洗效果较好。膜间地表由于蒸发作用的影响，使各离子在地表逐渐积累，Cl-移动性较强，在膜间地表积累程度最大。

3）活度系数估算方法对模拟结果的准确性有重要影响，尤其是盐分含量较高时，采用通用的函数关系可能会带来较大的模拟误差。