0 引 言
近年来，受温室效应的影响，我国极端天气出现频率屡增不减。其中，特大暴雨出现次数占据主要位置，且伴随城市化扩张及新型城镇化建设，大规模土地被现代化建筑物与硬质道路所替代，由此改变了城市原有区域土壤的自然本底与水文特征[1-2]，造成城市区域内部应对洪涝灾害的能力显著降低，给居民生活与财产带来巨大威胁。为此，2015年国务院办公厅出台《关于推进海绵城市建设的指导意见》指出，采用渗、滞、蓄、净、用、排等措施，将70%的降雨就地消纳和利用。透水砖铺装是海绵城市建设中应用较广的工程措施之一，因其造价低、适用范围广的特点备受国内外专家学者关注。现阶段在研究透水砖铺装措施对区域降雨入渗及产流控制方面，主要针对透水砖铺装系统中结构层所用材料、设计厚度、所处顺序及透水砖铺设类型等问题进行研究[3-5]，而针对地基层土壤特点对透水砖铺装系统的入渗及产流机制研究较少。
本文以海绵城市试点城市庄河市为例，通过渗流试验，构建不同特性土壤的土水特征曲线，依据土水特征曲线模型，得到各土壤入渗物性参数。建立不同雨强条件下的透水砖铺装系统雨水入渗模型，结合铺装系统各层几何和物理参数，形成透水砖铺装系统体积含水率时空分布和产流时间计算方法，对产流控制效果最佳时基层和垫层厚度进行优化，阐明地基层土壤特性对产流影响程度，从而为海绵城市透水砖铺装工程建设提供理论依据。
1 地基层土壤的土水特征曲线测定
1.1 试验土料与试验方法
本试验土壤取自海绵城市试点城市庄河市，将采集好的土壤进行自然晾晒、烘干，捣碎后放置筛分仪中，筛分出粒径为0.5～2 mm、0.1～0.25 mm、<0.075 mm的土样若干，表1为根据GBJ 123—1988《土工试验方法标准》测定得到的各类土壤的基本土性指标。
表1 各类土壤基本土性指标
Table 1 Basic properties of all types of soil
土壤编号土壤类型土壤相对密度土壤孔隙比土壤干密度/(g·cm-3)液限/%塑限/%饱和导水系数/(cm·s-1)饱和含水率/(cm3·cm-3)1号粗砂2.5781.321.11——9.60E-0456.942号细砂2.3151.251.03——2.87E-0455.513号细砂2.3961.201.09——1.75E-0452.924号细砂2.4331.121.15——1.19E-0450.325号粉质黏土2.5661.491.030.3700.2071.55E-0459.86
选用滤纸法对各类土壤进行土水特征曲线的测定，主要试验步骤按照土壤烘干—土样装填—水量配置—水分均匀—滤纸夹置—试样包裹—水量平衡—滤纸称重8个步骤依次进行。所用试验仪器主要为杭州新华造纸厂生产的“双圈牌”No.203型慢速滤纸、环刀(φ61.8 mm×20 mm)。
1.2 试验结果与分析
本文采用van Genuchten[6](简称VG模型)进行土水特征曲线拟合。VG模型如式(1)所示：
式中：θ为体积含水率，cm3/cm3；θr为残余含水率，cm3/cm3；θs为饱和含水率，cm3/cm3；α为进气值的倒数，cm-1；ψ为土壤基质吸力，cm；n和m为经验拟合参数，其中
图1为各类土壤土水特征曲线实测点与拟合曲线汇总图。可知：各土壤土水特征曲线经拟合后均呈现倒“s”形，体积含水率随着基质吸力的增大而减小。表2为各类土壤土水特征曲线的拟合参数结果。可知：各类土壤土水特征曲线经VG模型拟合后得到的R2值均>0.88，说明按该模型拟合效果较好，能够为后续模拟土壤入渗规律及产流控制研究奠定基础。
1号土实测点； ——1号土拟合曲线；
2号土实测点； —-·2号土拟合曲线；
3号土实测点； ……3号土拟合曲线； 
4号土实测点； —-—4号土拟合曲线； 5号土实测点； 
5号土拟合曲线。
图1 各类土壤土水特征曲线实测结果及拟合曲线
Fig.1 Experimental results and fitting curves of soil-water characteristic curves for various soils
表2 各类土壤土水特征曲线拟合参数
Table 2 Fitting parameters of soil-water characteristic curves for various soils
土壤编号αnmR21号0.065381.9430.485 0.88872号0.015962.2540.556 0.96623号0.017432.0570.514 0.97244号0.017241.9140.478 0.95925号0.010591.6550.396 0.9503
2 透水砖铺装产流模型的建立
2.1 雨水产流机理
基于CJJ/T 188—2012《透水砖路面技术规程》[7]，对透水砖铺装系统进行设计。图2为透水砖铺装系统示意图。可知：该系统自上而下依次为透水砖面层、找平层、基层、垫层、地基层，其中除地基层外，其余部分为系统的结构层，各层的设计高度在相应的位置标出。
图2 透水砖铺装系统示意
Fig.2 Schematic of the porous brick paving system
表3为透水砖饱和导水系数试验结果，该结果通过恒定水头法测得。可知：透水砖饱和导水系数均>10-2 cm/s，符合CJJ/T 188—2012。假设找平层与基层分别由干硬性水泥砂浆、无砂大孔隙混凝土构成，垫层采用1号土壤。表4为系统结构层设定的物理参数。
表3 透水砖饱和导水系数测量结果
Table 3 Experimental results of saturated water conductivity of porous brick
试验编号温度/℃时间/s高度/cm流量/cm3饱和导水系数/(×10-2 cm·s-1)12815143.01101.0422825141.61921.1132835140.12811.1742845138.33461.1352855136.34451.2162865133.95101.19
表4 透水砖铺装系统结构层物理参数
Table 4 Physical parameters of porous brick pavement structure layer
名称雨水吸附率/%初始含水率/(cm3·cm-3)有效孔隙率/%饱和含水率/(cm3·cm-3)透水砖面层9.89—30.00—找平层9.91—32.00—基层9.78—35.00—垫层—0.1490—0.5694
透水砖铺装系统内产流出现方式分为2种：一种是适用系统结构层中透水砖面层、找平层、基层各部分的蓄满产流模式；另一种是适用垫层与地基层各部分的超渗产流模式[8-9]。故以土壤为入渗介质的垫层与地基层，当该层表面入渗速率>外界雨强值时，不出现产流，雨强值为控制该层入渗过程的因素；当该层表面入渗速率<外界雨强值时，出现超渗产流情况，入渗速率为该层控制入渗过程的因素。而以大孔隙多孔材料为入渗介质的透水砖面层、找平层、基层，雨水通过其中的有效孔隙得以下渗，当水流初次通过该孔隙表面时，会形成少量水膜附着于表面，此刻该位置的体积含水率可近似为雨水吸附率，若底部排水及时，则该位置的体积含水率恒等于该层的雨水吸附率；反之，雨水则将该层的有效孔隙自下而上依次填充。
为方便模拟透水砖铺装系统的产流情况，入渗时只考虑各层湿润锋垂直运移变化情况，将整个系统简化为一维入渗。图3为系统产流时间计算的具体思路与步骤。可知：根据初始条件不同，系统产流计算大体分为2种情况，一种为只考虑垫层或地基层界面产流情况，另一种为考虑垫层与基层、地基层2个界面共同的产流情况。
图3 透水砖铺装系统产流时间计算流程
Fig.3 Flow chart for calculating runoff time in porous brick pavement system
2.2 入渗模型的建立
基于透水砖铺装系统各层的产流模式，该系统入渗模型分为2种：一种适用于透水砖面层、找平层、基层各部分，另一种适用于垫层与地基层各部分。
研究透水砖面层、找平层、基层部分入渗规律时，该介质层中的有效孔隙为雨水提供入渗通道，由于该介质层材料特性，对渗入的雨水只有初始时的吸附作用，完全依靠重力作用对该层各位置进行蓄满填充。依据水量平衡原理，出现降雨时，降雨总量主要由三部分分配：第1部分以基层与垫层界面为起点，雨水自该位置以重力与土壤基质吸力双重作用进行垂直入渗，该部分作为系统的入渗量；第2部分是当基层与垫层界面处入渗能力不及外界雨强时，雨水由此向上逐渐填充除垫层外的各结构层的有效孔隙，该部分作为系统的储水量；第3部分是当雨水填充至透水砖面层时，该层上方出现汇水，该部分作为系统的产流量。由此得出式(2)：
Q=W+S+R
(2)
式中：Q、W、S、R分别为垂直方向上累计降雨量，累计入渗量，累计储水量，累计产流量，cm。根据式(2)，可得入渗时间与各层雨水蓄满高度间的关系：
Li=S/(ei-ai)
(3)
式中：Li为垂直方向上蓄满高度，cm；ei、ai分别为有效孔隙率、雨水吸附率，%；i=1,2,3分别指代透水砖面层、找平层、基层。
研究垫层与地基层部分入渗规律时，可采用的入渗模型已发展相对成熟，主要包括Green-Ampt模型与Richards方程[10-13]。前者形式简易，但其假定湿润锋为恒值，故不能较好地反映土壤的实际入渗情况；而后者可通过解析解模拟得出土壤在入渗过程中各个阶段的入渗特征。本文在Richards方程的基础上加以改进，更便于模拟透水砖铺装系统在不同雨强情况下的入渗规律。
Richards方程是非饱和土壤水分运动的基本方程，它基于达西定律与能量守恒定理而得出，通过求解得出土壤剖面体积含水率随时空的分布情况、湿润锋位置的运移及土壤入渗速率随入渗时间的变化等结果[14-15]。
根据能量守恒定理，流入与流出微分单元体的质量变化应等于其水量的变化，得到方程：
式中：θj为土壤的体积含水率，cm3/cm3；t为水流通过微分单元体的时间，s；vnj为n方向上的流速，cm/s；n=1,2,3分别表示x、y、z方向；j=4,5分别指代垫层与地基层。根据达西定律可得式(5)：
式中：K(θj)为土壤的导水系数方程，是土壤基质吸力(或体积含水率)的函数；Hj为总土水势，为重力势zj与基质势ψj之和，cm。将式(4)与式(5)联立可得方程(6)：
根据式(6)，将其简化成一维垂直方向的土壤入渗水运动方程：
式中：zj为土壤深度，cm，向下为正。式(7)经变换可得式(8)：
式中：
表示扩散度。土壤扩散度是土壤导水系数与受土壤单位体积含水率变化引起的基质吸力变化率之积，该值主要取决于土壤体积含水率。将VG模型代入Burdine和Mualem提出的统计传导率模型，得到导水系数方程[16]：
K(θj)=ksj·Θj0.5[1-(1-Θj1/mj)mj]2
(9)
式中：ksj为饱和导水系数，cm/s；Θj为有效含水饱和度，其中
由于
式(1)求导后，可得：
故扩散度又可表示为式(9)与式(10)乘积的形式：
2.3 模型的求解
上述针对透水砖铺装系统各层材料介质的差异，建立了相应的入渗模型。通过式(2)可看出：为表征系统产流前(R=0)透水砖面层、找平层、基层处入渗规律，首先需得到垫层与地基处入渗时间与入渗量之间的关系。根据式(8)，得到土壤入渗速率：
式中：Ij为土壤入渗速率，cm/s。根据图3得出透水砖面层、找平层、基层处入渗时间与蓄满高度之间的表达式如式(13)、式(14)所示：
Li=(Q-Wj)/(ei-ai)
(13)
(14)
注：1)为地基层上界面产流；2)基层、地基层上界面共同产流；3)基层上界面产流。
式中：P为外界雨强值，cm/s。假设该介质层雨水入渗前各剖面处体积含水率为雨水吸附率，则初始条件为：
θi(z,0)=ai， li<z<l′i
(15)
式中：li与l′i分别为各层上、下界面处标高，cm。假设透水砖面层上界面处输入持续恒定的雨强P，则边界条件为：
vi(0,t)=P
(16)
综合式(14)求得透水砖面层、找平层和面层中体积含水率随时空分布关系式：
假设雨水入渗前垫层与地基层土壤各剖面处的体积含水率为初始含水率，则初始条件为：
θj(z,0)=θ0j， 0<z<lj
(18)
式中：lj为土壤最大埋深深度，cm；θ0j为土壤初始含水率，cm3/cm3。假设地表处标高为0，雨水入渗时地表处的体积含水率为饱和含水率，并当地下埋深趋向∞时，各剖面处的体积含水率为恒值，恒等于土壤的初始含水率，则边界条件为：
式(6)为非线性偏微分方程，求解复杂，且无法得到严格的解析解。现将扩散度函数作为恒值用平均扩散度Dj替代，然后根据
以
替代
将式(6)变为常系数的线性方程，便于求解，如式(20)所示：
结合式(20)，联立式(18)与式(19)，得到土壤剖面体积含水率随入渗时间与空间分布变化的表达式[17]：
式中：
为补余误差函数[18]。在土壤入渗时，外界雨强值是土壤剖面体积含水率随时空分布变化的重要影响因素，故将式(21)中的θsj以θtj替代，θtj以分段函数的形式给出，以饱和导水系数为分界点，根据外界雨强值的大小，作出相应判断：
式中：P为雨强值；tf为时间步长值，s；z为空间步长值，cm。
3 模拟结果与分析
3.1 透水砖铺装系统剖面体积含水率的时空分布
本文假定地基土壤层最大埋深为35 cm，以2号土壤作为地基层为例，模拟雨强在4 cm/h、0～4 h内系统降雨入渗及产流情况。图4为透水砖铺装系统各层剖面体积含水率的时空分布情况。可知：随入渗时间增大，地基层垂直各点处的体积含水率自土壤表层向下依次趋近于土壤饱和含水率；而地基层之上的结构层垂直各点处的体积含水率则自垫面层与地基层交界面位置向上依次蓄满，当蓄满至地表处的透水砖面层，此刻系统内出现产流。结合式(12)与式(13)可得，地基层土壤特性与结构层厚度是影响产流时间的主要因素。
A—透水砖部分； B—找平层部分； C—基层部分； D—垫层部分； E—地基部分。
图4 透水砖铺装系统各层剖面体积含水率时空分布
Fig.4 Temporal-spatial distribution of volume water content for different layers in porous brick pavement system
3.2 地基层土壤特性对透水砖铺装系统内产流时间的影响
地基层土壤特性包括土壤粒径、干密度值、初始含水率值等特性。假定基层厚度为15 cm，图5为地基层中不同粒径土壤在初始含水率为0.15 cm3/cm3时系统产流时间的对比。可知：在相同雨强下，土壤粒径越大，系统产流出现越滞后；随着雨强值的增大，系统内以同种土壤作为地基层的产流出现越提前。其中，以1号土壤作为地基层时，不同雨强下产流时间涨幅程度最大，其在雨强值为4 cm/h时，产流出现最迟，相比雨强值为6，8 cm/h时，产流时间分别各延长了8.1，16.7倍；而以2号、5号土壤作为地基层时，在不同雨强下产流时间差异均不明显。
1号土； 
2号土； 
5号土。
图5 地基层土壤粒径对系统产流时间的影响
Fig.5 Effect of soil particle size on system runoff time
图6为地基层中不同干密度的细砂在初始含水率为0.15 cm3/cm3时系统产流时间的对比。可知：在相同雨强下，随着土壤干密度值的增大，系统产流出现越滞后；在不同雨强下，以同种干密度细砂作为地基层时，系统产流时间随雨强的增大而提前。其中，雨强为4 cm/h时，以各干密度值细砂作为地基层时，产流时间差异较大，3号、4号土壤相对于2号土壤，系统产流时间分别提前42.25%、62.21%。
2号土； 3号土； 4号土。
图6 地基层土壤干密度值对系统产流时间的影响
Fig.6 Effect of soil dry density on system runoff time
图7为地基层中初始含水率不同的细砂，在干密度为1.03 g/cm3时系统产流时间的对比。可知：在相同雨强下，随着土壤初始含水率增加，系统产流时间越提前；在不同雨强下，以相同初始含水率细砂作为地基层时，系统产流时间随雨强值增大而越提前。其中，图中各初始含水率下的细砂在相同雨强下产流时间差异均不明显，雨强值为4 cm/h时对产流时间影响最大，初始含水率为0.15 cm3/cm3的细砂作为地基层，相比图中其余初始含水率的细砂作为地基层，产流时间仅分别增大9.79%、21.40%。
θc=0.15 cm3/cm3； θc=0.20 cm3/cm3； θc=0.25 cm3/cm3。
图7 地基层土壤初始含水率值对系统产流时间的影响
Fig.7 Effect of initial soil moisture content on system runoff time
综上，为研究各土壤特性对系统产流时间的影响程度，地基层以2号土为例作为参照土壤。由图5—图7可知：维持其余条件恒定，地基层土壤粒径对系统产流时间的影响程度最大，干密度值次之，初始含水率值的影响程度最小。
3.3 基层与垫层厚度对透水砖铺装系统内产流时间的影响
基层厚度在透水砖铺装系统结构层中占据较大比重，设计最优的基层厚度对于系统内产流控制起到关键作用。若维持图2中系统结构层中基层与垫层的总厚度恒定，考虑基层对透水砖铺装路面的承压作用，则根据CJJ/T 188—2012，基层厚度在10～16 cm为宜。依据此范围，模拟得出以不同土壤作为地基层时，系统内产流时间与基层厚度间的关系。图8为系统以不同土壤作为地基层时基层最优厚度的确定图。可知：随着地基层土壤粒径减小或土壤干密度增加，基层与垫层厚度比对系统产流控制的影响受到削减；在各雨强下、系统产流时间在4 h内，为使系统产流出现最迟，以1—4号土作为地基层时，基层与垫层厚度比约为4∶1；而以5号土作为地基层时，基层与垫层厚度比约为3∶2。可以看出，当系统产流控制效果达到最优时，粉质黏土作为地基层相比砂质土壤作为地基层，基层厚度减小约27%、垫层厚度增大约1倍。二者基层最优厚度差异较大，是因为粉质黏土粒径小且质量轻，当垫层厚度较小时，极易随雨水穿透垫层，从而造成基层中孔隙堵塞，不利于系统的产流控制。故增大垫层厚度能减小基层孔隙被小粒径土壤填充的概率，进而增强系统产流控制能力。
—P=4 cm/h； 
—P=6 cm/h； 
—P=8 cm/h。
图8 透水砖铺装系统产流时间随基层厚度的变化
Fig.8 The variation of runoff time with thickness of base course in porous brick pavement system
4 结 论
本文通过渗流试验，构建庄河市5类土壤的土水特征曲线，进而建立不同雨强条件下的透水砖铺装系统雨水入渗模型，形成透水砖铺装系统体积含水率时空分布和产流时间计算方法，阐明透水砖铺装系统结构和土壤特性对产流的影响，优化透水砖铺装系统基层和垫层厚度，为海绵城市透水砖铺装工程建设提供理论依据。具体结论如下：
1)土壤粒径对产流控制影响程度最大。对于土壤粒径较小的地基层，随着透水砖铺装系统中垫层厚度比的增加，可以适当延迟系统的产流时间，但延迟幅度较小。
2)在雨强为4 cm/h的特大暴雨时，庄河市粗砂土壤作为地基层相比细砂土壤和粉质黏土作为地基层，可以延长产流时间至少8倍；随着雨强继续增大，延长产流时间逐渐减弱。
3)对于庄河市特大暴雨时：以粉质黏土作为地基层时，基层与垫层厚度比约为3∶2；以砂质土壤作为地基层时，基层与垫层厚度比约为4∶1时，透水砖铺装系统产流控制效果最好。
[1] 黄敬军, 武鑫, 姜素, 等. 海绵城市建设的地质影响及适宜性评价：以徐州为例[J]. 地质论评, 2018, 64(6): 1472-1479.
[2] 陈珂珂, 何瑞珍, 梁涛, 等. 基于“海绵城市”理念的城市绿地优化途径[J]. 水土保持通报, 2016, 36(3): 258-264.
[3] Heweidak M, Amin S. Effects of OASIS®phenolic foam on hydraulic behaviour of permeable pavement systems[J]. Journal of Environmental Management, 2019, 230(15): 212-220.
[4] Palla A, Gnecco I, Carbone M, et al. Influence of stratigraphy and slope on the drainage capacity of permeable pavements: laboratory results[J]. Urban Water Journal, 2015, 12 (5): 394-403.
[5] Niu Z G, Lv Z W, Zhang Y, et al. Stormwater infiltration and surface runoff pollution reduction performance of permeable pavement layers[J]. Environmental Science and Pollution Research, 2016, 23 (3): 2576-2587.
[6] Amnyattalab J, Rezaie H. Study of the effect of seepage through the body of earth dam on its stability by predicting the affecting hydraulic factors using models of Brooks-Corey and van Genuchten (Case study of Nazluchay and Shahrchay earth dams)[J]. International Journal of Environmental Science and Technology, 2018, 15 (12): 2625-2636.
[7] 中华人民共和国住房和城乡建设部. 透水砖路面技术规程：CJJ/T188—2012[S]. 北京： 中国建筑工业出版社, 2012.
[8] 赵飞, 陈建刚, 张书函, 等. 透水铺装地面降雨产流模型研究[J]. 给水排水, 2010, 36(5): 154-159.
[9] 霍文博, 李致家, 李巧玲. 半湿润流域水文模型比较与集合预报[J]. 湖泊科学, 2017, 29(6): 1491-1501.
[10] Richards L A. Capillary conduction of liquids through porous mediums[J]. Physics, 1931, 1 (5): 318-333.
[11] Stewart R D. A dynamic multidomain green-ampt infiltration model[J]. Water Resources Research, 2018, 54 (9): 6844-6859.
[12] Vafakhah M, Fakher Nikche A, Sadeghi S H. Comparative effectiveness of different infiltration models in estimation of watershed flood hydrograph[J]. Paddy and Water Environment, 2018, 16 (3): 411-424.
[13] Cockett R, Heagy L J, Haber E. Efficient 3D inversions using the Richards equation[J]. Computers & Geosciences, 2018, 116: 91-102.
[14] 王康. 非饱和土壤水流运动及溶质迁移[M]. 北京： 科学出版社, 2010.
[15] Farthing M W, Ogden F L. Numerical solution of richards’ equation: a review of advances and challenges[J]. Soil Science Society of America Journal, 2017, 81 (6): 1257.
[16] Weber T K D, Iden S C, Durner W. Unsaturated hydraulic properties of Sphagnum moss and peat reveal trimodal pore-size distributions[J]. Water Resources Research, 2017, 53 (1): 415-434.
[17] Tan J Q, Song H G, Zhang H L, et al. Numerical investigation on infiltration and runoff in unsaturated soils with unsteady rainfall intensity[J]. Water, 2018, 10 (7): 914.
[18] Menziani M, Pugnaghi S, Vincenzi S. Analytical solutions of the linearized Richards equation for discrete arbitrary initial and boundary conditions[J]. Journal of Hydrology, 2007, 332 (1): 214-225.