SWAT模型的主要子模型介绍

SWAT主要通过水文过程模型、土壤侵蚀模型和污染负荷模拟模型等三类子模块进行水、泥沙和化学物质的模拟[99]。下面将结合2005版的SWAT理论[88]和相关文献简要介绍相应的子模型。

8.7.3.1　水文过程模型

SWAT采用的水量平衡方程如下

式中：SWt为时段末土壤含水量，mm；SW0为时段初土壤含水量，mm；t为计算时段，d；Rday，i为第i天的降雨量，mm；Qsurf，i为第i天的地表径流量，mm；Ea，i为第i天的蒸发量，mm；wseep，i为第i天的渗透量，mm；Qgw，i为第i天的地下径流量，mm。

1.地表径流

在产流计算中，SWAT利用较多是SCS曲线法，它是根据日土壤湿度状况来预测日降水所产生的地表径流。其降雨—径流关系表达式如下

式中：Qsurf为径流累计量，mm；Rday为日降雨量，mm；Ia为初损，mm，即地表径流产生之前的降雨损失；S为降雨前土壤潜在最大滞留量，mm。

S是一个描述土壤持水能力的参数，因土壤类型、土地利用情况、管理水平、地形坡度和不同时间等因素的不同而不同，该参数被定义为

其中CN是一个无量纲参数，称为曲线数，是反映降水前流域特征的一个综合参数，与流域前期土壤湿度（Antecedent Moisture Condition，AMC）、坡度、土地利用状况、土壤类型等因素有关。CN值越大，越容易产生径流，反之则产生径流越困难。CN值在0～100之间，当CN＝100，则S＝0，表示地表存在一层不透水层，如城市的水泥地面；当CN＝0，则S→∞，表示地表存在完全透水层，这种情况下不产生径流。

通过大量数据资料的分析，Ia与S有如下统计关系

于是，可得径流量计算公式为

如果Rday＜0.2S，则Qsurf＝0。

CN值可根据不同的土壤类型、土地利用和植被覆盖的组合根据所提供的表格进行查询[88]。但所给出的表中CN值是在坡度为5%和前期土壤湿度（AMC）为正常湿度状态的前提下得出的，于是，SWAT引入了CN值修正公式。

SCS定义了三类前期土壤湿度（AMC）状态：Ⅰ类为干旱（凋萎点），对应的CN值最小；Ⅱ类为正常湿度；Ⅲ类为湿润（田间持水量状态）。Ⅰ类和Ⅲ类状态下的CN值计算公式如下

式中：CN1、CN2和CN3分别为前期土壤湿度（AMC）在Ⅰ类、Ⅱ类和Ⅲ类状态下坡度为5%时的CN值。

SCS模型坡度调整公式为

式中：CN2s为前期土壤湿度（AMC）在Ⅱ类状态下对应区域平均坡度的CN值；CN2和CN3分别为前期土壤湿度（AMC）在Ⅱ类和Ⅲ类状态下坡度为5%时的CN值；slp为计算区域平均坡度实际值，%或m/m。

在汇流计算中，SWAT针对水文响应单元（HRU）分别计算河道和坡面汇流时间。

河道汇流时间计算公式为

式中：ct为河道汇流时间，h；L为河道长度，km；n为河道曼宁系数；A为水文响应单元（HRU）面积，km2；cs为河道坡度，%或m/m。

坡面汇流时间计算公式为

式中：ot为坡面汇流时间，h；sl为计算区域平均坡长，m；n为坡面曼宁系数；s为坡面坡度，%或m/m。

2.蒸散发

蒸散发包括冠层截留水蒸发、蒸腾和升华及土壤水的蒸发。蒸散发是水分转移出流域的主要途径，在许多流域与大陆（南极洲除外）蒸发量都大于径流量。准确地计算蒸散发量是估算水资源量的关键，也是研究气候和土地覆盖变化对河川径流影响的关键[99]。

（1）潜在蒸散发。SWAT提供了Penmam-Monteith、Priestley-Taylor和Hargreaves-Samani等三种计算潜在蒸散发能力的方法。另一方面，在SWAT中也可以直接使用实测资料或已经计算好的逐日潜在蒸散发资料。

Penmam-Monteith计算公式请参阅本文第三章相关章节。

Priestley-Taylor公式由Priestley和Taylor1972年提出，计算公式为[100]

式中：ET0为蒸散发能力，mm；αpet为常数，取αpet＝1.28；Δ为饱和水汽压斜率，k Pa/℃；λ为水的气化潜热，MJ/kg；γ为湿度计常数；Hnet为净辐射，MJ/（m2·d）；G为土壤热通量，MJ/（m2·d），由于土壤热能的储存和释放只有十几小时（早晨或上午吸收，下午或晚上就全部释放），对于以10～30d为计算尺度的模型模拟来说，土壤的热通量很小，通常可以忽略，因此在SWAT中，G＝0。

水的气化潜热λ计算公式为

式中：T为气温，℃。

饱和水汽压斜率Δ计算公式为

式中：ea为饱和气压，k Pa；T为气温，℃。

净辐射Hnet计算公式为

式中：Rs为太阳辐射量，MJ/m2；β为反照率。

β与土壤、作物和雪的覆盖有关：

1）当雪的厚度大于5mm时，β＝0.6。

2）当雪的厚度小于5mm且无作物生长时，β值则为土壤的反照率。

3）当有作物生长时，β为

其中

式中：COVs为土壤覆盖指数；CV为地面生物量与作物残余物之和，t/hm2；βs为土壤反照率。

Hargreaves-Samani计算公式为[101]

（2）实际蒸散发。在总的潜在蒸散发量计算结束后，就可以计算实际蒸散发量。SWAT首先计算植被冠层截留蒸发，然后计算最大蒸腾量、最大升华量和最大土壤水分蒸发量，最后计算实际的升华量和土壤蒸发量。

1）冠层截留蒸发。SWAT在计算实际蒸法时假定最大可能地蒸发冠层截留的水分，如果潜在蒸散发E0小于冠层截留的自由水量，则

式中：Ea为实际日蒸散发量，mm；Ecan为日冠层中自由水蒸发量，mm；E0为日潜在蒸散发量，mm；RINT（f）为日冠层中最终自由水量，mm；RINT（i）为日冠层中原有自由水量，mm。

如果潜在蒸散发E0大于冠层截留的自由水量，则

一旦冠层中截留的自由水分被全部蒸发后，继续蒸发所需要的水分将来至于植被、雪或土壤中。

2）植物蒸腾。如果假设植被生长在一个理想状态下，则植物蒸腾计算式可表示为

3）升华。如果水文单元中有积雪，将会发生升华，只有无积雪时，土壤中才会发生蒸发。因此，当水文单元中存在有积雪时，SWAT首先从积雪中升华水分以满足蒸发需求。如果积雪中的水含量大于最大的升华或土壤蒸发量时，则

如果积雪中的水含量小于最大的升华或土壤蒸发量时，则

4）土壤水分蒸发。当存在土壤水分蒸发时，SWAT首先要在不同土壤层中分配蒸发量，因此，土壤深度层次的划分将决定土壤允许的最大蒸发量，其计算公式为

式中：Esoil，z为在土壤深度z（mm）处的蒸发需求量，mm。

式（8-23）中系数确定的依据是使50%的蒸发需求量来自0～10mm土壤深度范围，而95%的蒸发需求量来自0～100mm土壤深度范围。

土壤水分蒸发量是由土壤上层蒸发量与决定的，即

式中：Esoil，ly为ly层的土壤蒸发量，mm；Esoil，zu为土壤上层蒸发量，mm；Esoil，zl为土壤下层蒸发量，mm。

由于前面已假设50%的蒸发量由0～10mm内的土壤层提供，如果是100mm的蒸发量，那么其中50mm的蒸发量都0～10mm内的土壤层来提供，试验表明这是无法满足的[88]。另一方面，当某一层土壤的蒸发需求未能满足时，SWAT不允许有另一层土壤来补偿这一需求，这必将导致蒸散发量的减少。为了满足蒸发需水量，SWAT给出了一个反映土壤毛细作用和裂隙等对修正系数，使土壤水分蒸发量计算公式调整为

式中：esco为土壤蒸发补偿系数。

esco反映了土壤毛细作用和裂隙等对不同土壤层蒸发需求量的影响，不同的esco对应相应的土壤层划分深度，具体参见文献[88]。随着esco的减小，SWAT允许从更深的土壤获得更多水分以供给蒸发。当土壤层含水量低于田间持水量时，蒸发需水量也相应减少，相应的计算公式为

另外，为了限制干旱状态下的蒸发量，SWAT给出了一个最大允许蒸发量，其值为植被可利用水量的80%，植被可利用水量即是土壤层中总含水量与凋萎点时土壤层含水量之差。最大允许蒸发量计算公式如下

5）融雪。在SWAT中，融雪过程通过积雪温度与大气最高温度的平均值和融雪开始时的最低温度之差的线性函数来描述，即

式中：SNOmlt为日积雪融化水量，mm；bmlt日积雪融化系数，mm/（d·℃）；snocov为水文单元中积雪覆盖面积比；Tsnow为日积雪温度，℃；Tmx为日大气最高温度，℃；Tmlt为融雪开始时的最低温度，℃。

积雪融化系数bmlt随季节变化而变化，最大值出现在夏天，最小值出现在冬天，其计算公式为

式中：bmlt6为6月21日的积雪融化系数；bmlt12为12月21日的积雪融化系数；dn为公历的第n天。

3.土壤水

土壤水可以因植物吸收或蒸腾而损耗，也可以渗漏到土壤底层最终补给地下水，还可以壤中流的形式补给地表径流。SWAT壤中流计算公式如下。

壤中流（Qlat）计算基本公式为

式中：Qlat为ly层土壤壤中流水量，mm/d；SWly，excess为ly层土壤可流出水量，mm；SWly为ly层土壤含水量，mm；FCly为ly层土壤田间持水量，mm；Ksat为土壤饱和导水率，mm/h；slp为流域平均坡度；φd为土壤层总空隙率φsoil与土壤水分含量达到田间持水量的空隙率φfc之差；Lhill为山坡坡长，m。

4.地下水

采用下列方程式计算地下水量

式中：Qgw，i为第i日进入河道的地下水量，mm；Qgw，i-1为第（i-1）天进入河道的地下水量，mm；Δt为时间步长，d；wrchrg，i为第i日蓄水层的补给量，mm；αgw为基流的消退系数；N为从退水开始所持续的时间，d；Qgw，N为第N天的地下水，mm；Qgw，0为开始消退时的地下水，mm；δgw为补给滞后时间，d；wseep为第i日通过土壤剖面底部进入地下含水层的水量，mm。

8.7.3.2　土壤侵蚀模型

SWAT对于土壤侵蚀过程是通过MUSLE（Modified Universal Soil Loss Equation）[102]来进行模拟计算的，在MUSLE中，用径流因子代替了USLE[97]中的降水动能函数，因此无需泥沙输移系数就可以预测泥沙产量，且可以用于单此暴雨的产沙预测分析。MUSLE水土流失计算方程为

式中：sed为土壤侵蚀量，t；Qsurf为地表径流量，mm/hm；qpeak洪峰流量，m3/s；areahru为水文单元面积，hm；KUSLE为土壤侵蚀因子；CUSLE植被覆盖和管理因子；PUSLE保持措施因子；LSUSLE为地形因子；CFRG为土壤粗糙度因子。

1.土壤侵蚀因子KUSLE

KUSLE反映土壤抵抗侵蚀能力的高低，它与土壤物理性质和组成结构密切相关。当土壤细沙（粒径0.002～0.05mm）和极细沙（粒径0.05～0.10mm）含量低于70%时，Wischmeier等[103]在1971年提出了计算土壤侵蚀因子公式

式中：M为土壤粒径参数；OM为有机物百分含量；csoilstr为土壤结构代码；cperm为土壤渗透特性等级；msilt为土壤细沙（粒径0.002～0.05mm）百分含量；mvfs为土壤极细沙（粒径0.05～0.10mm）百分含量；mc为土壤黏粒（粒径小于0.002mm）百分含量；org C为土层有机碳容量，%。

根据土壤自然结构和形状，在SWAT中给出了4种土壤结构代码csoilstr，具体可根据表8-7对照查询相应的csoilstr。对于极细颗粒的粒状土壤，csoilstr＝1；对于细颗粒的粒状土壤，csoilstr＝2；对于中等或粗颗粒的粒状土壤，csoilstr＝3；对于块状、板块状、棱形或大块状土壤，csoilstr＝4。

表8-7　土壤结构代码　单位：mm

cperm是根据土壤剖面最低饱和水力传导率给定的，具体见表8-8。

表8-8　土壤渗透特性等级

1995年Williams[102]又提出了计算KUSLE的一个替换公式，即

式中：fcsand为高粗沙含量土壤侵蚀因子低值或少沙土壤侵蚀因子高值；fcl-si为高细沙或粘粒含量土壤的侵蚀因子低值；forgc为高有机碳含量土壤侵蚀修正因子；fhisand为高含沙土壤侵蚀修正因子；ms为土壤沙粒（粒径0.05～2.00mm）百分含量，%；msilt为土壤细沙（粒径0.002～0.05mm）百分含量，%；mc为土壤黏粒（粒径小于0.002mm）百分含量，%；org C为土层有机碳容量，%。

2.植被覆盖和管理因子CUSLE

植被覆盖和管理因子CUSLE表示植物覆盖和作物栽培措施对防止土壤侵蚀的综合效应，其含义是在地形、土壤类型和降水条件相同的情况下，种植作物或林草的土地与连续休闲土地土壤流失量的比值，其最大值为1。SWAT采用的植被覆盖和管理因子CUSLE计算公式为

式中：CUSLE，mn为植被覆盖和管理因子最小值；CUSLE，aa为根据历史资料计算得到的植被覆盖和管理因子年平均值；rsdsurf为地表植物残余量，kg/hm2。

表8-9　等高耕耘和种植PUSLE值[97]

3.保持措施因子PUSLE

保持措施因子PUSLE是指有保护措施的地表土壤流失与不采取任何措施的土壤流失的比值，这里的保护措施包括等高耕耘、等高带状种植和梯田种植。等高耕耘和种植几乎能够完全抵御中低强度暴雨侵蚀，而对于超强暴雨侵蚀几乎是无能为力。等高耕耘和种植抵御土壤侵蚀的有效坡度范围为3%～8%。不同坡度等高耕耘和种植和等高带状种植的PUSLE值分别见表8-9和表8-10。

表8-10　等高带状种植的PUSLE值[97]

表8-10中，A是指4年轮作的行（排）农作物、小粒谷类作物和草地（2年）；B是指4年轮作的行（排）农作物（2年）、冬季谷类作物（2年）和草地（1年）；C是指其他的带状行（排）农作物和冬季谷类作物。

4.地形因子LSUSLE

地形因子LSUSLE计算公式为

式中：Lhill为坡长，m；αhill为坡度，（°）或rad；m坡长指数；slp为HRU的坡度，%

或m/m。

5.土壤粗糙度因子CFRG

土壤粗糙度因子CFRG计算公式为

式中：rock为第一层土壤岩石百分含量，%。

8.7.3.3　污染负荷模型

1.溶解态氮污染负荷模型(www.xing528.com)

硝态氮可随地表径流、侧向流和渗流在水体中迁移，因此，土壤流动水体中硝态氮浓度乘以土壤流动水体总量便是土壤中硝态氮流失总量。

硝态氮浓度计算公式为

式中：concNO3，mobile为ly层流动水体中硝态氮浓度，kg/mm；NO3ly为土壤ly层中硝态氮量，kg/hm；wmobile为土壤ly层中流动水体总量，mm；θe为孔隙度；SATly土壤ly层饱和水量，mm。

对wmobile分两种情况计算：

（1）对于地表10mm土层，有

（2）对于10mm以下土层，有

式中：Qsurf为地表径流，mm；Qlat，ly为土壤ly层侧向流，mm；wperc，ly为渗透到下层土壤水流（mm）。

随地表径流流失的硝态氮NO3surf可通过下式计算

式中：NO3surf为随地表径流流失的硝态氮，kg/hm；βNO3为硝态氮渗透系数；concNO3，mobile这里是指地表10mm土层的流动水体中硝态氮浓度，kg/mm。

随侧向流流失的硝态氮NO3lat，ly（kg/ha）可以分两种情况计算：

（1）对于地表10mm土层，有

（2）对于10mm以下土层，有

随渗流流失的硝态氮NO3perc，ly（kg/hm）可通过下式计算

2.吸附态氮污染负荷模型

吸附在土壤颗粒上的有机氮可随地表径流迁移，SWAT中采用1976年McElroy等[104]发展并由Williams和Hann[105]修正的负荷函数计算随水土流失的有机氮，其公式为

式中：orgNsurf为由地表径流传输到主干河道的有机氮量，kg/hm2；concorgN为地表10mm土壤中有机氮浓度，g/t；sed为产沙量，t；areahru水文响应单元面积，hm2；εN：sed为氮富集率；orgNfrsh，surf为地表10mm土壤中新鲜有机物中氮含量，kg/hm2；orgNsta，surf为土壤中稳定有机物中氮含量，kg/hm2；orgNact，surf为土壤中活性有机物中氮含量，kg/hm2；ρb第一层土壤容重，t/m3；depthsurf为土壤表层深度，即10mm；concsed，surq为地表径流中的泥沙浓度，t/m3；sed为产沙量，t；areahru水文响应单元面积，hm2，Qsurf为地表径流，mm。

3.溶解态磷污染负荷模型

溶解态磷在土壤中的迁移主要是通过扩散进行。地表10mm土壤溶解态磷输移计算公式为

式中：Psurf为通过地表径流流失的溶解态磷数量，kg/hm；Psolution，surf为地表10mm土壤溶解态磷数量，kg/hm2；Qsurf为地表径流，mm；ρb为地表10mm土壤容重，t/m3；depthsurf为土壤表层深度，即10mm；kd，surf为土壤磷分配系数，m3/t，即表层土壤（10mm）中溶解态磷的浓度与地表径流中溶解态磷的浓度的比值。

4.吸附态磷污染负荷模型

有机磷和矿物质磷是通过吸附在土壤颗粒上并通过径流迁移，SWAT中同样采用1976年McElroy等[104]发展并由Williams和Hann[105]修正的负荷函数计算随水土流失的有机磷，其公式为

式中：sed Psurf为通过地表径流流失的有机磷数量，kg/hm2；concsedP为地表10mm土壤中吸附在泥沙上有机磷的浓度，g/t；sed为产沙量，t；areahru水文响应单元面积，hm2；εP：sed为磷富集率；min Pact，surf为地表10mm土壤中活性矿物中磷含量，kg/hm2；min Psta，surf为地表10mm土壤中稳定矿物中磷含量，kg/hm2；org Phum，surf为地表10mm土壤中腐殖有机体中磷含量，kg/hm2；org Pfrsh，surf为地表10mm土壤中新鲜有机物中磷含量，kg/hm2；ρb第一层土壤容重，Mg/m3；depthsurf为土壤表层深度，即10mm；concsed，surq为地表径流中的泥沙浓度，Mg/m3；sed为产沙量，t；areahru水文响应单元面积，hm2；Qsurf为地表径流，mm。

5.河道中各种氮的变化

在有氧的水环境中，氮的存在形式可以从有机氮转化到氨，然后到亚硝酸盐、硝酸盐。

河道中有机氮的数量可随海藻中氮转化为有机氮而增加，而河道中有机氮浓度也会因有机氮转化为氨氮或有机氮吸附到沉淀物上而减少。一天内有机氮的变化为

式中：ΔorgNstr为有机氮变化量，mg/L；α1为海藻生物量中的氮含量，mg/mg；ρa为海藻死亡率，d-1；algae为一天开始时藻类生物含量，mg/L；βN，3为有机氮转化为氨的速度常数，d-1；orgNstr为一天开始时有机氮含量，mg/L；σ4为有机氮沉淀系数，d-1；TT为河段水流运动时间，d；βN，3，20为20℃时有机氮转化为氨的速度常数，d-1；Twater为日（或小时）平均水温，℃；σ4，20为20℃时有机氮沉淀系数，d-1；Twater为日（或小时）平均水温，℃。

河道中氨氮的数量可随有机氮氧化或来之河道沉积物扩散氨氮而增加，而河道中氨氮浓度也会因氨氮转化为亚硝酸盐或海藻的吸收而减少。一天内氨氮的变化为

式中：ΔN H4str为氨氮变化量，mg/L；βN，3为有机氮转化为氨氮的速度常数，d-1；orgNstr为一天开始时有机氮含量，mg/L；βN，1为氨氮生物氧化速度常数，d-1；NH4str为一天开始时氨氮含量，mg/L；σ3为沉淀物氨氮释放率，mg/（m2·d）；depth为河道水深，m；frNH4为藻类氨氮吸收系数；α1为海藻生物量中的氮含量，mg/mg；μa为藻类生长速度，d-1；algae为一天开始时藻类生物含量，mg/L；TT为河段水流运动时间，d；βN，1，20为20℃时氨氮氧化速度常数，d-1；Oxstr为河道溶解氧浓度，mg/L；Twater为日（或小时）平均水温，℃；σ3，20为20℃时为沉淀物氨氮释放率，mg/（m2·d）；fNH4为氨氮偏好因子；NH4str为河道氨氮浓度，mg/L；NO3str为河道硝酸盐浓度，mg/L。

河道中亚硝酸盐的数量可随氨氮转化为亚硝酸盐而增加，而随亚硝酸盐向硝酸盐的转化而减少。一天内亚硝酸盐的变化为

式中：ΔNO2str为亚硝酸盐变化量，mg/L；βN，1为氨氮生物氧化速度常数，d-1；NH4str为一天开始时氨氮含量，mg/L；βN?，2为由亚硝酸盐氧化为硝酸盐的速度，d-1；NO2str为一天开始时亚硝酸盐含量，mg/L；TT为河段水流运动时间，d；βN，2，20为20℃时由亚硝酸盐氧化为硝酸盐的速度常数，d-1；Oxstr为河道溶解氧浓度，mg/L；Twater为日（或小时）平均水温，℃。

河道中硝酸盐的数量可随亚硝酸盐的氧化而增加，而河道中硝酸盐的浓度会因海藻的吸收而减小。一天内硝酸盐的变化为

式中：ΔNO3str为硝酸盐变化量，mg/L；βN，2为由亚硝酸盐氧化为硝酸盐的速度，d-1；NO2str为一天开始时亚硝酸盐含量，mg/L；frNH4为藻类氨氮吸收系数；α1为海藻生物量中的氮含量，mg/mg；μa为藻类生长速度，d-1；algae为一天开始时藻类生物含量，mg/L；TT为河段水流运动时间，d。

6.河道中各种磷的变化

死亡海藻将海藻磷转化为有机磷，有机磷矿化后又成为能够被海藻吸收的可溶解态磷，这就是磷循环。有机磷可通过吸附在沉淀物上在河道进行迁移，一天内有机磷的变化为

式中：Δorg Pstr为有机磷变化量，mg/L；α2为海藻生物量中的磷含量，mg/mg；ρa 为海藻死亡率，d-1；algae为一天开始时藻类生物含量，mg/L；βP，4为有机磷矿化的速度，d-1；org Pstr为一天开始时有机磷含量，mg/L；σ5为有机磷沉淀系数，d-1；TT为河段水流运动时间，d；βP，4，20为20℃时有机磷矿化速度常数，d-1；Twater为日（或小时）平均水温，℃；σ5，20为20℃时有机磷沉淀系数，d-1。

河道中无机磷/溶解态磷可随有机磷的矿化和河道沉淀物中无机磷的扩散的增加而增加，而河道中溶解态磷的浓度会因海藻的吸收而降低。一天内无机磷/溶解态磷的变化为

式中：Δsol Pstr为无机磷/溶解态磷变化量，mg/L；βP，4为有机磷矿化速度常数，d-1；org Pstr为一天开始时有机磷含量，mg/L；σ2为沉淀物溶解态磷释放率，mg/（m2·d）；depth为河道水深，m；α2为海藻生物量中的磷含量，mg/mg；μa为藻类生长速度，d-1；algae为一天开始时藻类生物含量，mg/L；TT为河段水流运动时间，d；σ2，20为20℃时沉淀物溶解态磷释放率，mg/（m2·d）；Twater为日（或小时）平均水温，℃。