Water-salt phase transition of saline soils during evaporation
-
摘要: 伴随蒸发过程，盐渍土孔隙内的溶液会产生汽化及结晶现象，且两种不同的相变间存在相互依赖性。为了揭示蒸发过程中的水盐相变特性，对盐渍土在蒸发过程中的汽化和盐结晶行为进行了理论和试验研究。首先，从热力学理论出发，通过考虑相间化学势平衡和Young-Laplace方程给出了起始结晶半径的理论表达，并分析了起始结晶半径与温度和相对湿度的关系。结合Van Genuchten土水特征曲线模型，建立了孔隙溶液浓度在不同状态下盐的溶解-结晶模型，并分析了温度和初始含盐量对盐结晶行为的影响规律。最后，通过开展控制温度和初始含盐量的蒸发试验结果对理论模型进行了验证。结果表明：环境湿度降低的过程中，土体温度以及初始含盐量对孔隙溶液中的水分和盐分相变有着显著的影响，土体温度越高，孔隙溶液中水分蒸发越快，液相最终饱和度越低；随着初始含盐量的增加，液相最终饱和度下降，盐晶体体积比增加。初始含盐量还会改变盐的盐析湿度，并且盐晶体的产生会抑制水分的蒸发。Abstract: Along with the evaporation process, the solution in the pores of the saline soils will be vaporized and crystallized, and the two different phase transitions are interdependent. In order to reveal the water-salt phase change characteristics during evaporation, the vaporization and salt crystallization behaviors of the saline soils during evaporation are studied theoretically and experimentally. Firstly, starting with the thermodynamic theory, the theoretical expression for the initial crystallization radius is given by considering the interphase chemical potential balance and Young-Laplace equation, and the relationship among the initial crystallization radius, temperature and relative humidity is analyzed. Combined with the Van Genuchten soil-water characteristic curve model, the salt dissolving-crystallization model under different states of pore solution concentration is established, and the influences of temperature and initial salt content on the salt-crystallization behavior are analyzed. Finally, the theoretical model is verified by the evaporation experiments with controlled temperature and initial salt content. The results show that the soil temperature and initial salt content have significant effects on the phase transition of water and salt in the pore solution during the process of ambient humidity reduction. When the soil temperature is higher, the water evaporation becomes faster and the liquid phase saturation is lower. With the increase of the initial salt content, the final liquid phase saturation decreases and the salt crystal volume ratio increases. The initial salt content also changes the humidity of salting out, and the generation of salt crystals will inhibit the evaporation of water.
-
Keywords:
- saline soil /
- salt crystal /
- evaporation test /
- solution concentration /
- relative humidity
-
0. 引言
降雨入渗是造成滑坡最为重要的诱因之一[1]，降雨入渗会降低土体抗剪强度，增加土体重度，进而导致边坡下滑力增大，诱发边坡失稳。发展一种有效的降雨入渗模型进而在最大程度上降低降雨对边坡的危害，是降雨作用下边坡稳定性分析的关键。目前常用的边坡降雨入渗分析方法主要有[2]：①基于Richards方程的数值分析方法；②基于物理入渗模型的解析计算方法，如Green-Ampt（GA）入渗模型、Smith模型和Philip模型等。数值分析方法需要不断迭代求解降雨入渗控制方程，计算过程较为繁琐，对于复杂问题将十分耗时，并且若是网格尺寸和时间步长选择不当，将会导致计算不收敛。相较之下，以GA入渗模型为代表的解析计算方法，模型参数具有明确的物理意义，计算过程便捷，在降雨入渗下边坡稳定性研究中得到了广泛应用[3-14]。如Chen等[3]基于GA入渗模型研究了边坡降雨入渗过程，推导了边坡降雨入渗率计算表达式。张杰等[4]考虑气阻作用建立了分层假定的雨水入渗计算分析的GA模型。Zhang等[5]基于GA模型建立了可以计算时变破坏概率和失效时间的可靠度分析框架，并开发了基于电子表格的降雨型边坡破坏概率预测方法。Dou等[6]考虑土体饱和渗透系数变异性，基于GA入渗模型建立了雨水入渗-重分布分析模型。潘永亮等[7]建立了可综合考虑地下水位、初始含水率及不同降雨工况的改进GA模型。蒋水华等[8]提出了一种可考虑土体参数空间变异性的修正GA模型，并探讨了降雨入渗与多参数空间变异性作用下边坡失稳机理。
虽然上述研究大大地拓展了GA入渗模型的应用范围，但是忽略了土壤中客观存在的湿润区非饱和过渡层。大量试验发现湿润锋至土体入渗表面之间的土体并非处于完全饱和状态，在饱和区和天然区之间至少存在一个过渡层，并且过渡层内不同深度处土体体积含水量也不尽相同[9-10]。为此，张杰等[11]基于土体分层假定，推导了GA模型降雨历时与湿润锋深度之间的函数关系。胡海军等[12]在分层GA模型基础上，合理考虑湿润锋处吸力作用和水分剖面形状，分析了积水和降雨作用下非饱和重塑黄土的降雨入渗过程。Yao等[13]将饱和层与非饱和层厚度设定为入渗区深度的50%，分别对无限长边坡和马家沟滑坡进行了稳定性分析，并发现分层GA模型对研究强降雨事件下边坡稳定性问题十分有效。温馨等[14]针对黄土入渗中存在的非饱和区，提出了一种修正的GA模型来表征土体实际渗透系数随埋深的变化规律。
另一方面，受沉积、后沉积等地质作用的影响，岩土体参数存在固有的空间变异性。现有研究[8, 15-16]表明土体水力参数的空间变异性是影响降雨入渗下边坡渗流计算的关键因素。为客观评价降雨入渗下边坡失稳机理及稳定性，需要合理表征岩土体参数空间变异性的影响。然而，目前缺乏适用于考虑土体渗透系数空间变异性多层非均质边坡的降雨入渗模型。为此，提出了一种改进的Green-Ampt模型，通过4种常见的边坡工况验证该模型的有效性，进而探讨降雨入渗下多层非均质边坡含水率分布及稳定性问题。
1. 改进Green-Ampt入渗模型
1.1 Green-Ampt入渗模型
Green-Ampt（GA）入渗模型是1911年由Green等[17]提出的，最初用于解决初始干燥土体在薄层积水时的入渗问题。该模型假定雨水在入渗时存在可将湿润区和非湿润区明确分开的水平湿润锋面（图 1）。GA模型假设土体体积含水量呈矩形分布，湿润区为土体饱和层，其体积含水量为饱和体积含水量θs，非湿润区为天然层，其体积含水量为初始体积含水量θi。
![]()
图 1 经典GA入渗模型的含水率分布Figure 1. Distribution of water content of classic Green-Ampt model据此，土体入渗率f可表示为
f=−ksdhdz=−ks(hi−zf)−h0zf−0=ks(1+h0−hizf)=ks(1+h0+Sfzf)。 (1) 式中：ks为土体饱和渗透系数；h0为地面处积水水头；zf为湿润锋深度；hi为湿润锋下方的水头；Sf为湿润锋处土体概化基质吸力水头，是一个关于含水量的函数，可通过相关试验或下面的推导方程进行求解[18]：
Sf=∫1Kr(θi)ψ(Kr)dKr1−Kr(θi)=3λ+23λ+1ψb[1−Kr(θi)3λ+13λ+2]1−Kr(θi)≈3λ+23λ+1ψb。 (2) 式中：ψb为土体进气值；Kr为湿润锋处相对饱和渗透系数；λ为表征土体孔径分布特征的拟合系数。进而，可获得累积降雨入渗量I为
I = (θs−θi)zf。 (3) 可见，该模型计算简单，且模型参数具有明确的物理意义。然而，针对边坡降雨入渗问题，该模型存在如下不足：①在模拟降雨入渗时，该模型没有考虑到降雨强度R对入渗率f的影响，将会高估入渗量I，导致计算的湿润锋深度偏大。②该模型没有考虑客观存在的过渡层，大量试验研究[19-21]表明，湿润锋至土体入渗表面之间的土体并非完全处于饱和状态，即在饱和区和天然区之间应至少存在一个过渡层。如不考虑客观存在的过渡层，将会高估湿润锋以上的土重，低估边坡安全系数。③当只考虑参数（如饱和渗透系数）变异性时，传统的GA模型可以采用临界渗透系数划分降雨入渗边界条件，但是一旦考虑土体参数空间变异性，整个土层剖面的渗透系数呈空间非均匀分布，不同位置处数值大小不同，此时采用临界渗透系数划分降雨入渗边界条件不再适用[8]。
1.2 改进Green-Ampt入渗模型
为克服上述不足，十分必要对经典GA入渗模型进行改进，使其能够处理考虑土体参数空间变异性的多层非均质边坡降雨入渗问题。实际工程边坡或天然边坡通常存在一定的坡度，因而积水易形成坡面径流，其中的关键难题是如何确定不同时刻下边坡降雨入渗边界条件。鉴于采用临界渗透系数划分降雨入渗边界条件不适用考虑参数空间变异性的多层非均质边坡，为此蒋水华等[8]提出采用开始积水时间tp划分入渗边界，本文将该GA模型称为“修正GA模型”。由GA入渗模型[11, 13]可知，土体入渗率f随着累计入渗量I的增加而逐渐减小。对于倾角为α的边坡而言，当土体的累计入渗量达到某一累积入渗量Ip值时，边坡坡面开始积水[5, 7]，此时f = Rcosα，Ip的计算表达式为
Ip=(θs−θi)Sfcosα(R/ks−1)。 (4) 式中：Ip为积水时刻tp所对应的累积入渗量；α为边坡倾角；R为降雨强度。据此可得边坡开始积水的时间tp为
tp=Ip/(Rcosα)。 (5) 根据开始积水时间tp可划分入渗边界条件如下：
（1）当降雨历时t < tp时，入渗区土体处于非饱和状态，雨水将全部入渗至土体中，此时的入渗边界由流量控制，可得到入渗控制方程为
f=Rcosα=k(θ)+ks[ψr(θ)−ψr(θi)]zf。 (6) 式中：θ为土体体积含水量；ψr(θ)为湿润锋处土体相对吸力水头[8]，
ψr(θ)=∫−∞ψk(θ)ksdψ=∫−∞ψ(ψbψ)3λ + 2dψ=(ψbψ)3λ + 1ψb3λ + 1=Se3+1/λψb3λ + 1。 (7) 式中：θr为土体残余体积含水量；Se为有效饱和度，Se=(θ−θr)/(θs−θr)。
随后，依据GA模型的矩形入渗原理可推导相应的湿润锋深度为
zf=Iθ−θi = Rtcosαθ−θi。 (8) 将式（8）代入式（6）中，可得到非线性方程：
f=Rcosα=k(θ)+ks[ψr(θ)−ψr(θi)](θ−θi)Rtcosα。 (9) 求解式（9）所示的非线性方程，可得到不同降雨历时下土体体积含水量θ分布，进而结合式（8）可以计算得到对应的湿润锋深度。
（2）当降雨历时t ≥ tp时，此时边坡表面开始积水，因降雨强度超过土体允许入渗率，故土体表面入渗区处于饱和状态，此时的入渗边界由水头控制[2, 6, 8]。降雨入渗过程可分为自由入渗和积水入渗两个阶段。其入渗率f和累计入渗量I可为[11, 13]
f(t)={Rcosα(t⩽tp)ks(cosα+Sf/zf)(t>tp), (10) I=∫t0f(t)dt。 (11) 按照式（11）进行积分得到累积入渗量，进而得到入渗饱和区湿润锋深度zf为
zf=Iθs−θi。 (12) 因此，一旦确定了土体饱和渗透系数ks和相应的降雨历时t，便可针对上述两种降雨入渗边界，分别求解对应的土体湿润锋深度和含水量分布。
另一方面，大量相关试验研究[9-10, 19]表明：湿润锋至土体入渗表面之间的土体真实体积含水量并非完全处于饱和状态，即在饱和层和天然层之间应至少存在一个过渡层。根据文献[9~13]，将GA入渗模型进行分层处理，增加客观存在的过渡层，如图 2（b）所示。从图 2中可以看出，改进后的土体体积含水量分布包含饱和层、过渡层和天然层3部分。其中饱和层深度为zs，过渡层深度为zt。
![]()
图 2 改进的GA入渗模型的含水量分布Figure 2. Distribution of water content of improved Green-Ampt model文献[9，10]的相关试验数据表明，采用椭圆曲线分布来表征过渡区含水量分布具有较高的计算精度，故本文也采用椭圆曲线来描述过渡层含水量分布，由此可得修正的土体含水量分布为
θ(z)={θs (饱和层, 0≤z≤zs)θi+(θs−θi)√1−((z−zs)2zt2) (过渡层, zs≤z≤zh)θi (天然层, z>zh)。 (13) 式中：zh为整个入渗区深度，
zh=zs+zt。 (14) 随着降雨入渗的推进，饱和层、过渡层和天然层的深度都会发生变化。彭振阳等[10]通过试验发现过渡层深度与总入渗区深度zh存在一个比例关系，并且该比例与整个入渗区深度之间具有良好的线性关系：
zt=ηzh, η=azh+b。} (15) 式中：η为过渡层占整个入渗区深度的比例；a和b为拟合系数。
在过渡层修正过程中，应保持累计入渗量不变，即I=IM。I可通过式（8）和（11）计算得到，IM为增加过渡层后的累计入渗量，其值应等于饱和层入渗量Is和过渡层入渗量It之和。其中对于饱和层，入渗量Is =(θs−θi)zs；对于过渡层，由于其为椭圆曲线分布，故其入渗量等于1/4的椭圆面积It = 0.25π(θs−θi)zt，据此可得总的累计入渗量为
I=IM=Is+It=(θs−θi)zs+0.25π(θs−θi)zt。 (16) 联立式（14）~（16），便可计算得到饱和层深度zs和过渡层深度zt，进而可求得整个边坡土层的入渗区深度zh及体积含水量分布。本文将所提模型称为“改进GA模型”
2. 多层非均质边坡稳定性分析
降雨诱发滑坡通常是沿边坡浅层面发生，在剪切力作用下趋于破坏，滑动面深度远小于滑动面长度，故常采用图 3所示的无限长边坡模型分析浅层滑坡稳定性。根据莫尔-库仑破坏准则，基于极限平衡分析的无限长边坡模型安全系数计算公式为[15, 22]
FS = τfτm = c′+[(σn−ua)−σs]tanφ′Wsinαcosα。 (17) 式中：τf为土的抗剪强度；τm为沿潜在破坏面任意一点的剪应力；W为单位土体切片的重量；c′，φ′分别为土体有效黏聚力和有效内摩擦角；(σn−ua)为单位土条底部净法向力，ua为孔隙气压力。根据文献[2，6，8，15]，取ua = 0，σn=γtzcos2α，其中γt为土体有效重度，z为滑动面深度，α为无限长边坡倾角；σs为吸应力，
σs=−Seψ(θ)。 (18) 式中：ψ(θ)为基质吸力，是关于土体体积含水量θ的函数。
为了考虑降雨入渗过程中因体积含水量变化而引起的土体单位重量变化，对滑动面以上土体部分进行积分得到土体重量为
W=∫z0γtdz。 (19) 式中：土体有效重度γt[15]为
γt=γd+θ(z)γw。 (20) 式中：γd和γw分别为土体干重度和水重度。
综上，基于本文提出的改进GA模型进行降雨入渗下非均质边坡稳定性分析的步骤如下：
（1）收集并确定土体参数统计特征（包括均值、标准差、自相关函数、波动范围及概率分布等），沿垂直方向将边坡划分为H/d个等厚的均质土层（H为坡高，d为随机场单元尺寸）。
（2）基于土体参数统计特征，采用蒙特卡洛模拟（MCS）方法产生N组独立标准正态随机样本，随后基于Karhunen-Loève（KL）级数展开方法[15-16, 23]生成N组土体参数随机场实现值，并将其一一赋值给对应的土层。
（3）采用式（5）计算开始积水时间tp，并与降雨历时t进行比较，划分降雨入渗边界条件，再根据式（8），（12），（13）~（16），求解考虑客观存在过渡层的降雨入渗下边坡体积含水量分布。
（4）采用式（17）计算边坡每个土层及湿润锋处安全系数，从中可获得（H/d+1）个安全系数，并将其中最小的安全系数视作该边坡安全系数，对应的滑动面则为该边坡最危险滑动面。
3. 算例分析
下面以图 3所示的非饱和无限长边坡为例[6, 24]，验证提出的改进GA模型的有效性。边坡坡高H = 3 m，倾角α = 50°，下部为不透水基岩，降雨强度R恒定为5 mm/h。采用Brook等[25]提出的土水特征曲线模型描述土体渗透系数和基质吸力与体积含水量之间的函数关系为
k(θ)=ksSe3+2/λ, (21) ψ(θ)=ψbSe−1/λ。 (22) 式中：k(θ)为土体渗透系数。根据文献[6，24]，得到土体物理力学参数取值如表 1所示。另外根据文献[10]的试验结果，取式（15）中的系数a和b分别为-0.003，0.8712。
表 1 土体物理力学参数取值Table 1. Values of physical parameters of soil计算参数 取值 计算参数 取值 饱和渗透系数ks/(mm·h-1) 3 有效黏聚力c′/kPa 5 有效内摩擦角φ′/(°) 28 进气值ψb/kPa 2.752 土体干重度γ /(kN·m-3) 16.217 湿润锋处土体概化基质吸力Sf /mm 424.3 饱和含水量θs 0.335 残余含水量θr 0.068 初始含水量θi 0.148 土体孔隙分布特征参数λ 0.319 通过以下4种常见的边坡工况，验证采用本文提出的改进GA模型进行降雨入渗下边坡体积含水量计算及稳定性分析的有效性：①均质边坡；②两层边坡（上层渗透性大、下层渗透性小）；③两层边坡（上层渗透性小、下层渗透性大）；④考虑任意一次参数随机场典型实现的非均质边坡。同时，选用Richards方程数值解和修正GA入渗模型[8]计算结果进行对比分析。其中，Richards方程数值解采用一维流动模型来模拟边坡垂直入渗过程，描述一维水流在边坡上的入渗过程控制方程为[26]
∂θ∂t=∂∂z[k(θ)(∂h∂z−cosα)]。 (23) 式中：θ∈[θr,θs]；h为压力水头。
本文借助Hydrus-1D软件[27]对式（23）所示Richards方程进行数值求解，可得到不同降雨历时t下不同深度z处压力水头h(z, t)和土体体积含水量θ(z, t)分布。本文将Richards方程数值解视作精确解进行判断对比分析。
3.1 工况1：均质边坡
基于表 1中的土体参数均值，对图 3所示均质边坡模型进行降雨入渗分析。图 4比较了采用修正GA模型和提出的改进GA模型计算的体积含水量分布和安全系数及Richards方程数值解。
![]()
图 4 不同降雨历时下边坡体积含水量和安全系数分布的比较（工况1）Figure 4. Comparison of distribution of water content and factors of safety of slope under different rainfall durations (Scenario Ⅰ)由图 4可知，采用修正GA模型计算的体积含水量分布与提出的改进GA模型计算结果和Richards方程数值解均存在较大差别，这是由于修正GA模型假定雨水呈矩形入渗，没有考虑客观存在的过渡层。在此基础上进行边坡稳定性分析，其边坡FS沿坡体深度的变化如图 4（b）所示。由图 4可知，随着降雨的进行，雨水入渗区也不断下移，饱和层内土体基质吸力基本丧失，导致入渗区内边坡FS随着埋深的增大而逐渐减小。到达过渡层后，由于土体并非完全处于饱和状态，基质吸力的作用未完全丧失，边坡FS先下降而随后又呈增加趋势。接着，虽然土体天然层内土体含水量和基质吸力均保持初始状态，但由于土体深度不断增加，导致土条自重逐渐增加，坡体滑动力增大，此时的边坡FS表现为随着埋深的增大而逐渐降低，并在埋深z = 3.0 m的不透水层处达到最小值。
表 2比较了不同方法计算的降雨历时20，36，60 h下边坡安全系数和最危险滑动面深度。表中，入渗区是指雨水入渗至坡体内的区域，包括饱和层和过渡层（例如降雨历时60 h，图 4（a）中入渗区可视为距离坡面0~1.25 m区域）；整个边坡则包含入渗区和天然区。将其分开进行处理，可分析不同降雨历时下边坡最危险滑动面位置和确定边坡失稳是否因降雨所诱发。若边坡最危险滑动面位于入渗区，可认为边坡失稳因降雨造成。由表 2可知，边坡安全系数随着降雨历时的增加而减小。相较于修正GA模型，提出的改进GA模型与基于Richards方程数值方法计算的安全系数吻合。如降雨历时60 h，提出的改进GA模型计算的边坡安全系数为1.22，与基于Richards方程数值方法计算的1.19差别很小。表 2中采用修正GA模型计算的边坡安全系数偏小的原因是，该模型假定雨水呈矩形入渗，导致计算的体积含水量θ偏大，进而高估了土条重量W，造成边坡安全系数偏小。相较之下，提出的改进GA模型因考虑了客观存在的过渡层，计算的边坡体积含水量、安全系数以及最危险滑动面位置均与Richards方程数值解吻合，这表明提出的改进GA模型能够有效分析均质边坡降雨入渗问题。此外，降雨历时20，36 h时，3种方法均识别出边坡最危险滑动面位于z = 3.0 m不透水层处。降雨历时60 h，相较于修正GA模型计算结果（1.03 m），本文提出的改进GA模型结果（0.99 m）与精确解（0.99 m）更吻合，说明提出方法能够准确识别边坡最危险滑面位置。
表 2 边坡安全系数和最危险滑动面深度的比较（工况1）Table 2. Comparison of factors of safety and depths of critical slip surface of slope (Scenario Ⅰ)类别 降雨历时/h Richards方程数值解 修正GA模型 提出的改进GA模型 入渗区安全系数 20 2.79 2.35 2.78 36 1.72 1.53 1.74 60 1.19 1.10 1.22 整个边坡安全系数 20 1.36 1.36 1.36 36 1.34 1.34 1.34 60 1.19 1.10 1.22 最危险滑动面深度 20 3.00 m 3.00 m 3.00 m 36 3.00 m 3.00 m 3.00 m 60 0.99 m 1.03 m 0.99 m 3.2 工况2：两层边坡（上层渗透性大、下层渗透性小）
为了进一步说明提出的改进GA模型的有效性，本节对两层边坡（上层渗透性大、下层渗透性小）进行降雨入渗分析。第一层土饱和渗透系数为3.5 mm/h，土层厚度为0.5 m；第二层土饱和渗透系数为3 mm/h，土层厚度为2.5 m。两个土层的其他参数均相同，见表 1。图 5比较了工况2下采用上述2种模型和Richards方程数值方法计算的土体体积含水量分布。由图 5可知，提出的改进GA模型计算结果同样与Richards方程数值解吻合。表 3比较了采用不同方法计算的降雨历时20，36，60 h安全系数及最危险滑动面深度。可见，提出的改进GA模型与Richards方程数值方法计算的边坡安全系数依然吻合。如降雨历时20 h，Richards方程数值方法计算的入渗区安全系数为2.83，提出的改进GA模型结果为2.81，而修正GA模型结果为2.39，进一步体现了提出的改进GA模型的有效性及考虑客观存在的过渡层的必要性。
![]()
图 5 不同降雨历时下边坡体积含水量分布的比较（工况2）Figure 5. Comparison of distribution of water content of slope under different rainfall durations (Scenario Ⅱ)表 3 边坡安全系数和最危险滑动面深度的比较（工况2）Table 3. Comparison of factors of safety and depths of critical slip surface of slope (Scenario Ⅱ)类别 降雨历时/h Richards方程数值解 修正GA模型 提出的改进GA模型 入渗区安全系数 20 2.83 2.39 2.81 36 1.71 1.55 1.76 60 1.19 1.10 1.22 整个边坡安全系数 20 1.36 1.36 1.36 36 1.34 1.34 1.34 60 1.19 1.10 1.22 最危险滑动面深度 20 3.00 m 3.00 m 3.00 m 36 3.00 m 3.00 m 3.00 m 60 0.98 m 1.03 m 0.99 m 3.3 工况3：两层边坡（上层渗透性小、下层渗透性大）
本节进一步分析两层边坡（上层渗透性小、下层渗透性大）来验证改进GA模型的有效性。第一层土饱和渗透系数为3 mm/h，土层厚度为0.5 m；第2层土饱和渗透系数为3.5 mm/h，土层厚度为2.5 m。这两个土层的其他参数均相同，见表 1。图 6比较了采用上述2种模型和Richards方程数值方法计算的边坡体积含水量分布。表 4比较了基于不同方法计算的边坡安全系数。与工况2类似，提出的改进GA模型与Richards方程数值方法计算的安全系数一致，说明提出的改进GA模型可有效求解多层边坡降雨入渗问题。
![]()
图 6 不同降雨历时下边坡体积含水量分布的比较（工况3）Figure 6. Comparison of distribution of water content of slope under different rainfall durations (Scenario Ⅲ)表 4 边坡安全系数和最危险滑动面深度的比较（工况3）Table 4. Comparison of factors of safety and depths of critical slip surface of slope (Scenario Ⅲ)类别 降雨历时/h Richards方程数值解 修正GA模型 提出的改进GA模型 入渗区安全系数 20 2.79 2.35 2.78 36 1.72 1.53 1.74 60 1.19 1.10 1.22 整个边坡安全系数 20 1.36 1.36 1.36 36 1.34 1.34 1.34 60 1.19 1.10 1.22 最危险滑动面深度 20 3.00 m 3.00 m 3.00 m 36 3.00 m 3.00 m 3.00 m 60 1.01 m 1.04 m 1.00 m 3.4 工况4：考虑一次参数随机场典型实现的非均质边坡
工程实际中，土体饱和渗透系数ks呈空间非均质分布。本文将边坡沿垂直方向上划分为60个均质土层，厚度为d = 0.05 m[6]。基于KL级数展开方法来离散ks随机场，利用高斯型自相关函数模拟土体ks的空间自相关性，垂直自相关距离lv取0.5 m。当KL级数展开截断项数取6时，计算的期望能比率因子ε为95.67%，可以满足随机场离散精度要求（大于95%）。图 7给出了土体饱和渗透系数随机场的一次典型实现，图中颜色较深部分表示ks值较大区域，颜色较浅部分表示ks值较小区域。
![]()
图 7 饱和渗透系数随机场一次典型实现分布图Figure 7. Typical realization of random filed of hydraulic conductivity图 8比较了不同方法计算的降雨历时8，36，60 h边坡体积含水量分布。由图 8可知，在3个不同降雨历时下，提出的改进GA模型计算的体积含水量分布和Richards方程数值解同样吻合，入渗区和整个边坡的安全系数和最危险滑动面深度几乎相同（见表 5），同样说明了基于提出的改进GA模型进行多层非均质边坡降雨入渗分析的有效性。由表 5可知，降雨历时8，36 h，由这2种模型和Richards方程数值方法计算的整个边坡安全系数皆小于入渗区安全系数，此时边坡最危险滑动面均位于z = 3 m的不透水层处。然而，在降雨历时60 h时，提出的改进GA模型和Richards方程数值方法均表明边坡最危险滑动面位于z = 3 m不透水层处。而基于修正GA模型计算的整个边坡安全系数等于入渗区安全系数，即修正GA模型认为此时最危险滑动面处于入渗区（距离坡面约0.85 m），显然错误评估了边坡失稳模式。这是由于该典型实现下土层表面渗透性低，湿润锋推进缓慢，虽然湿润锋处安全系数在不断减小，但是始终大于不透水层处安全系数，故而边坡仍然沿着不透水层处发生失稳破坏。综上，采用修正GA模型会造成对边坡失稳机理的错误认识，相比之下采用提出的改进GA模型能够准确识别最危险滑动面位置，获得更合理的边坡安全系数。另外需要说明的是，提出的改进Green-Ampt模型和Richards方程数值解得到的过渡层分布存在误差，并且该误差随着降雨历时的增加而增大，见图 4（a），5，6，这是因为本文没有考虑式（15）中系数a和b随着降雨历时的增加有可能发生变化。尽管如此，两者计算的湿润锋深度随着降雨历时的增加均保持一致。
![]()
图 8 工况4不同降雨历时下边坡体积含水量分布比较Figure 8. Comparison of distribution of water content of slope under different rainfall durations (Scenario Ⅳ)表 5 边坡安全系数和最危险滑动面深度的比较（工况4）Table 5. Comparison of factors of safety and depths of critical slip surface of slope (Scenario Ⅳ)类别 降雨历时/ h Richards方程数值解 修正GA模型 提出的改进GA模型 入渗区安全系数 8 6.58 5.62 6.43 36 1.95 1.72 1.98 60 1.40 1.24 1.39 整个边坡安全系数 8 1.37 1.37 1.37 36 1.35 1.35 1.35 60 1.33 1.24 1.33 最危险滑动面深度 20 3.00 m 3.00 m 3.00 m 36 3.00 m 3.00 m 3.00 m 60 h 3.00 m 0.85 m 3.00 m 为了说明这4种工况下边坡体积含水量分布之间的差别，图 9比较了这4种工况降雨历时36，60 h的边坡含水量分布曲线。由图 9（a）可知，工况1，2，3的边坡含水量分布曲线重合，而工况4入渗区较小。这是由于工况4边坡表面渗透性小，雨水入渗率较低，进而导致入渗区小。工况1，3由于上层饱和渗透系数相同，故其在降雨历时36 h下含水量分布曲线重合。工况2与工况1，3的边坡含水量分布曲线重合的原因是，在降雨初期，降雨历时未超过积水时刻tp，雨水入渗率等于Rcosα，此时的入渗量仅与降雨历时相关，与土体饱和渗透系数大小无关。因此，工况1，2，3在降雨历时36 h入渗量相差不大，对应的边坡含水量分布曲线基本重合。相比之下，降雨历时60 h，工况4由于其自身的渗透性较低，故入渗区小；工况3下由于下层土体渗透性大，导致入渗量增大，进而入渗深度更大；工况1，2由于下层土体的饱和渗透系数一致，故这两种工况不仅入渗率和入渗量相等，而且边坡体积含水量分布也相同，见图 9（b）。
![]()
图 9 不同工况下边坡体积含水量分布的比较Figure 9. Comparison of distribution of water content of slope for different scenarios4. 结论
考虑降雨入渗作用下边坡土体中客观存在的过渡层，提出了一种适用于多层非均质边坡降雨入渗分析的改进GA模型，为合理确定不同降雨历时下多层非均质边坡的入渗区深度和土体体积含水量分布奠定了基础。同时，以无限长边坡模型为例，进行降雨入渗下边坡稳定性分析，并与采用Richards方程数值解和修正GA模型计算的边坡安全系数和最危险滑动面深度进行了对比分析，验证了提出的改进GA模型的有效性。
（1）修正GA模型局限于雨水呈矩形入渗，没有考虑雨水入渗过程中土体客观存在的过渡层，导致降雨入渗过程中单位土体重度被高估，计算的边坡安全系数偏小，造成对边坡失稳机理的错误认识。
（2）提出的改进GA模型因考虑了边坡土体入渗区内客观存在的过渡层，故计算的土体含水量分布、安全系数和滑动面分布不仅贴近客观实际，而且与Richards方程数值解吻合，能够较好地解决考虑土体饱和渗透系数空间变异性的多层非均质边坡降雨入渗问题。
（3）相比于修正GA模型，提出的改进GA模型能够准确识别非均质边坡最危险滑动面位置，为降雨型滑坡加固设计及风险防控提供了一条有效途径。如工况4中，降雨历时60 h时，修正GA模型认为最危险滑动面位于入渗区，而改进GA模型和Richards方程计算结果均证明边坡最危险滑动面位于不透水层处。
（4）工程实际中存在多种环境因素对边坡降雨入渗过程具有重要的影响，包括边坡初始状态（如含水量，土壤类型等）的不均匀分布、边坡倾角大小、覆土厚度和降雨雨型等，融入多种环境因素影响的边坡降雨入渗模型还有待进一步研究。
-
![]()
图 1 不同温度下起始结晶半径随相对湿度的变化曲线
Figure 1. Curves of initial crystallization radius with relative humidity at different temperatures
![]()
图 2 不同温度下液相饱和度和盐晶体体积比与相对湿度的变化曲线
Figure 2. Variation curves of liquid phase saturation, salt crystal volume ratio and relative humidity at different temperatures
![]()
图 3 不同摩尔分数下液相饱和度和盐晶体体积比随相对湿度的变化情况
Figure 3. Variation of liquid saturation and salt crystal volume ratio with relative humidity under different mole fractions
![]()
图 4 不同温度条件下盐析湿度与初始摩尔分数的变化情况
Figure 4. Variation of salting-out humidity and initial mole fraction under different temperatures
![]()
图 5 在不同温度条件下试样溶液饱和度随着蒸发时间的变化
Figure 5. Variation of saturation of sample solution with evaporation time at different temperatures
![]()
图 6 在不同温度条件下试样液相饱和度和盐晶体体积比随相对湿度的变化情况
Figure 6. Variation of saturation of sample liquid and salt crystal volume ratio with relative humidity under different temperatures
![]()
图 7 不同初始含盐量下试样液相饱和度和盐晶体体积比随相对湿度的变化曲线
Figure 7. Variation curves of saturation of sample liquid and salt crystal volume ratio with relative humidity under different initial salt contents
![]()
图 8 试样液相饱和度计算结果和试验值的对比情况
Figure 8. Comparison between calculated results and test values of liquid phase saturation of sample
表 1 土的物理性质指标
Table 1 Physical properties of soil
土粒相对密度 液限
wL/%塑限
wP/%塑性指数IP 天然含水率w/% 不均匀系数
CU=d60/d102.71 26.8 17.6 9 9.3 6.87 
下载: 导出CSV
表 2 试验方案
Table 2 Test schemes
试验 溶质初始摩尔分数 温度/℃ ΔhR/% 试验Ⅰ 0 10, 20, 30 -5 试验Ⅱ 对应温度饱和初始摩尔分数 30, 20, 10 -5 试验Ⅲ 0.04, 0.03, 0.02, 0.01, 0 30 -5
下载: 导出CSV
表 3 在饱和溶液条件下试样液相饱和度计算结果和试验值的对比情况
Table 3 Comparison between calculated results and test values of saturation of liquid phase of samples in saturated solution
温度 相对湿度/% 90 80 70 60 50 40 30 20 30℃ 计算值 0.9623 0.8462 0.7126 0.5685 0.4726 0.3809 0.2843 0.2097 试验值 0.9725 0.8492 0.7042 0.5742 0.4642 0.3717 0.2892 0.2133 误差/% 1.06 0.35 1.18 1.01 1.78 2.41 1.72 1.71 20℃ 计算值 0.9479 0.8356 0.7175 0.5846 0.4943 0.4128 0.3465 0.2756 试验值 0.9541 0.8383 0.7133 0.5967 0.4992 0.415 0.3416 0.2717 误差/% 0.65 0.32 0.58 2.07 0.99 0.53 1.41 1.41 10℃ 计算值 0.9498 0.8316 0.7282 0.6175 0.5243 0.4462 0.3758 0.3062 试验值 0.9475 0.8383 0.7208 0.6133 0.5217 0.4417 0.3708 0.3033 误差/% 0.24 0.81 1.02 0.68 0.50 1.01 1.33 0.95
下载: 导出CSV
表 4 在欠饱和溶液条件下试样液相饱和度计算结果和试验值的对比情况
Table 4 Comparison between calculated results and test values of saturation of liquid phase of samples in undersaturated solution
摩尔分数 相对湿度/% 90 80 70 60 50 40 30 20 0.04 计算值 0.9142 0.7621 0.6376 0.5209 0.4359 0.3572 0.2855 0.2102 试验值 0.9225 0.7533 0.6317 0.5242 0.4317 0.3525 0.2800 0.2125 误差/% 0.91 1.15 0.92 0.63 0.96 1.32 1.93 1.09 0.03 计算值 0.9104 0.7023 0.5786 0.4814 0.4172 0.3413 0.2851 0.2251 试验值 0.9050 0.7075 0.5725 0.4858 0.4108 0.3433 0.2821 0.2208 误差/% 0.59 0.74 1.05 0.91 1.53 0.59 1.05 1.91 0.02 计算值 0.8811 0.6854 0.5401 0.4321 0.3772 0.3272 0.2764 0.2275 试验值 0.8858 0.6933 0.5450 0.4367 0.3725 0.3233 0.2733 0.2258 误差/% 0.53 1.15 0.91 1.06 1.25 1.19 1.12 0.75 0.01 计算值 0.8761 0.6708 0.5306 0.4387 0.3508 0.2907 0.2495 0.2041 试验值 0.8683 0.6783 0.5358 0.4350 0.3583 0.2967 0.2458 0.2008 误差/% 0.89 1.09 0.98 0.84 2.14 2.06 1.48 1.62
下载: 导出CSV
-
[1] HILLEL D. Introduction to Soil Physics[M]. Pittsburgh: Academic Press, 1982.
[2] 唐朝生, 施斌, 顾凯. 土中水分的蒸发过程试验研究[J]. 工程地质学报, 2011, 19(6): 875-881. doi: 10.3969/j.issn.1004-9665.2011.06.012 TANG Chaosheng, SHI Bin, GU Kai. Experimental investigation on evaporation process of water in soil during drying[J]. Journal of Engineering Geology, 2011, 19(6): 875-881. (in Chinese) doi: 10.3969/j.issn.1004-9665.2011.06.012
[3] GRAY D M. Handbook on the Principles of Hydrology[M]. NewYork: Water Information Center, inc, 1973.
[4] DENISOV Y M, SERGEEV A I, BEZBORODOV G A, et al. Moisture evaporation from bare soils[J]. Irrigation and Drainage Systems, 2002, 16(3): 175-182. doi: 10.1023/A:1021247218535
[5] NASSAR I N, HORTON R. Water transport in unsaturated nonisothermal salty soil: Ⅱ theoretical development[J]. Soil Science Society of America Journal, 1989, 53(5): 1330-1337. doi: 10.2136/sssaj1989.03615995005300050005x
[6] SHIMOJIMAA E, YOSHIOKA R, TAMAGAWA I. Salinization owing to evaporation from bare-soil surfaces and its influences on the evaporation[J]. Journal of Hydrology, 1996, 178(1/2/3/4): 109-136.
[7] SHOKRI N, OR D. What determines drying rates at the onset of diffusion controlled stage-2 evaporation from porous media?[J]. Water Resources Research, 2011, 47(9): W09513.
[8] NOROUZI RAD M, SHOKRI N. Nonlinear effects of salt concentrations on evaporation from porous media[J]. Geophysical Research Letters, 2012, 39(4): 4403.
[9] GUPTA S, HUININK H P, PRAT M, et al. Paradoxical drying of a fired-clay brick due to salt crystallization[J]. Chemical Engineering Science, 2014, 109: 204-211. doi: 10.1016/j.ces.2014.01.023
[10] SHOKRI-KUEHNI S M S, NOROUZI RAD M, WEBB C, et al. Impact of type of salt and ambient conditions on saline water evaporation from porous media[J]. Advances in Water Resources, 2017, 105: 154-161. doi: 10.1016/j.advwatres.2017.05.004
[11] NOROUZI RAD M, SHOKRI N, SAHIMI M. Pore-scale dynamics of salt precipitation in drying porous media[J]. Physical Review E, 2013, 88(3): 032404. doi: 10.1103/PhysRevE.88.032404
[12] SHOKRI N. Pore-scale dynamics of salt transport and distribution in drying porous media[J]. Physics of Fluids, 2014, 26(1): 191-209.
[13] JAMBHEKAR V A, HELMIG R, SCHRÖDER N, et al. Free-flow-porous-media coupling for evaporation-driven transport and precipitation of salt in soil[J]. Transport in Porous Media, 2015, 110(2): 251-280. doi: 10.1007/s11242-015-0516-7
[14] HUININK H P, PEL L, MICHELS M A J. How ions distribute in a drying porous medium: a simple model[J]. Physics of Fluids, 2002, 14(4): 1389-1395. doi: 10.1063/1.1451081
[15] SGHAIER N, GEOFFROY S, PRAT M, et al. Evaporation-driven growth of large crystallized salt structures in a porous medium[J]. Physical Review E- Statistical, Nonlinear, and Soft Matter Physics, 2014, 90(4): 042402. doi: 10.1103/PhysRevE.90.042402
[16] ESPINOSA-MARZAL R M, SCHERER G W. Mechanisms of damage by salt[J]. Geological Society, London, Special Publications, 2010, 331(1): 61-77. doi: 10.1144/SP331.5
[17] ZHANG C M, LI L, LOCKINGTON D. Numerical study of evaporation-induced salt accumulation and precipitation in bare saline soils: mechanism and feedback[J]. Water Resources Research, 2014, 50(10): 8084-8106. doi: 10.1002/2013WR015127
[18] COUSSY O. Deformation and stress from in-pore drying-induced crystallization of salt[J]. Journal of the Mechanics and Physics of Solids, 2006, 54(8): 1517-1547. doi: 10.1016/j.jmps.2006.03.002
[19] 万旭升, 赖远明, 张明义, 等. 土中未冻含水量与温度关系研究[J]. 铁道学报, 2018, 40(1): 123-129. doi: 10.3969/j.issn.1001-8360.2018.01.019 WAN Xusheng, LAI Yuanming, ZHANG Mingyi, et al. Research on relationship between unfrozen water content in soil and temperature[J]. Journal of the China Railway Society, 2018, 40(1): 123-129. (in Chinese) doi: 10.3969/j.issn.1001-8360.2018.01.019
[20] QIU E X, WAN X S, QU M F, et al. Estimating unfrozen water content in frozen soils based on soil particle distribution[J]. Journal of Cold Regions Engineering, 2020, 34(2): 04020002. doi: 10.1061/(ASCE)CR.1943-5495.0000208
[21] WILSON G W, FREDLUND D G, BARBOUR S L. The effect of soil suction on evaporative fluxes from soil surfaces[J]. Canadian Geotechnical Journal, 1997, 34(1): 145-155. doi: 10.1139/t96-078
-
期刊类型引用(5)
1. 郭旭辉，朱鸿鹄，吴冰，高宇新，胡乐乐，曹鼎峰. 基于人工神经网络的黄土含水率光纤被动感测技术研究. 岩土力学. 2025(02): 653-664 . 
百度学术
2. 高磊，袁泽，王勤，高明军. 基于OFDR技术的土体含水率模型研究. 工程地质学报. 2025(01): 96-105 . 百度学术
3. 严辉，林沛元. 深圳市岩溶地层标准贯入击数神经网络模型. 地质科技通报. 2025(02): 305-321 . 百度学术
4. 李杰，朱鸿鹄，吴冰，刘喜凤，王家琛，曹鼎峰，施斌. 下蜀土降雨入渗光纤感测及渗透系数估算研究. 工程地质学报. 2024(02): 601-611 . 百度学术
5. 朱鸿鹄. 工程地质界面：从多元表征到演化机理. 地质科技通报. 2023(01): 1-19 . 百度学术
其他类型引用(0)
-
其他相关附件
计量
- 文章访问数: 258
- HTML全文浏览量: 85
- PDF下载量: 70
- 被引次数: 5
