HBV 可作为一个全分布或半分布模型,将流域划分为子流域。在集总模型中,假设研究区域(流域)是一个单一的单元(区域),且参数在流域内不发生空间变化。HBV 模型由四个主要模块组成:
(1) 融雪积雪模块(Snowmelt and snow accumulation module);
(2)土壤水分和有效降水模块(Soil moisture and effective precipitation module);

(3) 蒸散发模块(Evapotranspiration module);
(4) 径流响应模块(Runoff response module)。
下图说明了 HBV 模型的简化版本的一般过程。
该模型可以按日或月时间步运行;所需的输入数据包括每个时间步的降水和温度观测的时间序列,以及长期估计的月温度和潜在蒸散速率的平均值。
如上图所示:
该模型包括一个基于每个时间步的输入温度,将输入的降水处理为降雨或雪的模块。然后在土壤湿度模块中处理降雨和融雪(如果存在的话),在那里对有助于地表径流的有效降雨进行评估。降雨的剩余部分有助于土壤水分的储存,只要地下有足够的含水量,这些水分本身就可以蒸发。模型的主要输出是流域出口处的径流,由地表径流(surface runoff)、互流(interflow,来自近地表流的贡献)和基流(baseflow,来自地下水流的贡献)三部分组成。该模型有许多参数需要根据现有的观测结果进行校准。下面简要介绍每个模块的细节。
2 融雪积雪模块融雪和积雪被认为与温度成正比。
第一个模型参数为阈值温度 ;温度高于 时,雪融化,低于 时,雪堆积。
将初始 温度设置为零摄氏度(32华氏度)是一个合理的初始假设。
如果在温度低于 时发生降水事件,则降水以雪的形式积累,否则则假定输入降水为降雨(即液体)。只要温度保持在阈值温度以下,输入的降水就不会促进径流。然而,一旦温度超过阈值,融雪和降水都开始对径流做出贡献。在该模型中,用公式1估计融雪速率(snowmelt rate)作为水当量(water equivalent):
,融雪速率等于水当量 ,DD,度日因子(degree-day factor) ,T,日平均气温 ,,积雪融化起始的阈值温度,
经验参数度日因子(DD)是指积雪含水率因高于冰点一度而在一天内减少的量。度日因子范围较广,从0.7到0.9 。度日因子DD可以假设为常数或可变参数。当降雨发生在现有的雪上时,由于稍暖的雨水中提供的额外热能,融雪量会增加。DD因子是可以通过现场试验测量的模型参数。这个模型参数也可以通过对观测到的水位图(hydrograph)的校准来估计。
3 有效降水和土壤水分落在流域上的降水通常分为两部分:
第一部分有助于入渗到土壤区;第二部分有助于地表径流。第二部分通常称为有效降水,由 HBV 根据降水时的土壤含水量估计。
田间容量(Field capacity,FC) 是描述地表下最大土壤含水量的参数。一般来说,当土壤含水量或降水量越高,降水对径流产生的贡献就越大。当土壤含水量接近田间容量时,入渗减少,降雨对径流产生的贡献增加。
公式2计算了有效降水作为当前土壤含水量的函数:
,有效降水,,真实土壤湿度,,日降水深度,,模型参数(形状参数)
对于给定的土壤水分亏缺(soil moisture deficit,通过SM/FC的比值测量),称为形状系数的参数 控制了有助于径流的液态水(liquid water,P + SM)的数量。
下图绘制了土壤湿度、田间容量、形状系数和径流系数(runoff coefficient)之间的关系,径流系数定义为有效降水与总有效水深的比值()。
从图中可以看出,对于特定的土壤湿度, 越高,径流系数越低。
随着土壤湿度(SM)接近田间容量(FC),径流系数增大。田间容量FC和形状系数都被用作标定参数(calibration parameters)。
结果表明,径流系数和土壤湿度不是恒定的,它们在模拟的时间步骤中是动态变化的。需要一个土壤湿度的初始值才能开始计算。
利用式2和土壤湿度(SM)的初始值,计算出有效降水量。例如,如果径流系数估计为0.7,那么70%的降雨(rainfall)有助于径流,其余(30%)渗入地下(infiltrates into the subsurface)。然后根据入渗和蒸散量(infiltration and evapotranspiration)更新土壤湿度的初始值。对于下一个时间步骤,使用新的土壤湿度值,并使用新的降水(precipitation)重复计算。
注:降雨量是指从天空降落到地面上的雨水,未经蒸发、渗透、流失而在水面上积聚的水层深度,一般以毫米为单位,它可以直观地表示降雨的多少。气象预报把下雨、下雪都叫做降水,降水的多少叫降水量,表示降水量的单位通常用毫米。1毫米的降水量是指单位面积上水深1毫米。
4 蒸散(Evapotranspiration)为了计算流域的实际蒸散发,模型用户需要提供长期的月平均潜在蒸散发()作为输入。然后,在模拟周期内的每一天,减去根据白天平均温度与长期月平均温度的差值计算得到的潜在蒸散量得到调整后的潜在蒸散量:
,调整后的潜在蒸散量,,日平均温度,,长期月平均气温,,长期月平均潜在蒸散量,是模型参数,
当日平均温度与长期平均温度相差较大时,模型参数 可以改善模型的性能。通过使用土壤永久萎蔫点(Soil Permanent Wilting Point,PWP)来耦合土壤湿度和实际蒸散发计算。
公式4为土壤湿度与实际蒸散发的关系:
,真实潜散发,,土壤永久萎蔫点,
由式4可知,当土壤湿度高于 PWP 时,实际蒸散发与潜在蒸散发发生速率相同。PWP 是土壤水分对蒸散发的限制,即当土壤水分小于 PWP 时,实际蒸散发小于调整后的蒸散发。
换句话说,由于PWP以下土壤水分有效性不足,公式4减少了蒸散发量。
下图给出了公式4中实际蒸散发量与PWP之间的关系。从图中可以看出,当PWP接近田间容量时,实际蒸散发会更高,反之亦然。在观测的基础上,模型参数 和PWP都可以通过模型标定得到。
5 径流响应
该模块基于水库概念估算流域出口处的径流。该系统由两个概念性储层组成,一个在另一个之上,如图所示。第一个水库用于模拟近地表流(the near surface flow),第二个水库用于模拟基流(地下水贡献,the base flow,groundwater contribution)。
从时序角度来看,第一个和第二个水库分别模拟快速和缓慢的地下过程。通过使用恒定的渗流速率(percolation rate,),储层之间直接相连。
如图所示,每个水库有两个出口(和)。当上游水库水位超过阈值 时,上游水库迅速产生径流()。另外两个出口(outlets)的流量响应相对较慢。
衰退系数(recession coefficients) 、、 分别表示上、下水库的响应函数。为保证径流过程最快, 的初始值应始终大于 。第三出口()的响应应该比第二出口()慢,因此 应该小于 。三个衰退系数和渗流速率均为模型参数,经校正估计。
公式5给出了上图所示出口的响应函数:
为近地表流量(near surface flow),为互流(interflow),为基流(baseflow),为渗流(percolation),为近地表流量蓄积系数,为互流蓄积系数,为基流蓄积系数,为渗流蓄积系数,为上游水库水位,为下游水库水位,为水位阈值,为流域面积,。
总模拟径流量 可由第一水库和第二水库流量之和()得到。
而完整HBV模型还包括其他功能,如使用更详细的流量转换功能和流量路由功能。这些特征没有包括在模型的教学版本中,因为这里的重点更多地是在主要的水文过程和理解它们的相对意义和相互作用上。
6 模型校准和性能评估如前所述,在概念模型中,水文过程的数学公式用一些参数加以简化。通常,参数分为两类:
描述物理特性的物理基础参数,如田间容量(FC)、土壤永久萎蔫点(PWP)和度日因子(DD)经验性的,只有很少或没有物理基础,但用于概念性地描述过程。该模型的经验参数包括水库参数(、、、、)、形状系数()和参数 。在校准过程中,迭代改变模型参数,直到实测流量与模拟流量(observed and simulated discharge)达到满意的匹配。如果基于物理的参数是从现场测量或经验中知道的,那么使用经验参数来校准模型就足够了。否则,所有参数都将包括在模型校准中。
在模型标定和评价模型性能(即实测流量与模拟流量是否一致)时,常采用均方根误差、峰值误差、一致性指数、纳什-萨克利夫系数、相关系数和相对累积误差等不同的统计准则。
这里使用NashSutcliffe:
为Nash-Sutcliffe系数为模拟流量 ,为观测流量,为平均观测流量 ,为时间步数
Nash-Sutcliffe系数(),范围在 到1之间,其中模型越接近1,模型越准确。
另一方面,负值意味着观测值的平均值是比被评估模型更好的预测指标。
7 python实现import numpy as npimport datetimedef module_snowmelting(t, precipitations, temp, snow_pack, rainfall, F=1, temp_threshold=0): """ Snow melting module Returns the amount of rainfall and snow accumulation. ------------------------------------------------------------ Parameters: t, date - Running day precipitations, array - Expected precipitations for the running period temp, array - Expected temperatures for the running period snow_pack, array - Amount of snow accumulated F, float - Degree-day factor temp_threshold, float - Temperature threshold for snow to melt """ if temp[t] > temp_threshold: # snow melts and precipitations are rainfalls rainfall[t] = precipitations[t] + F(temp[t]-temp_threshold) snow_pack[t-1] snow_pack[t] = snow_pack[t-1] (1- F(temp[t]-temp_threshold) ) else: # there is no rainfall, precipitations contribute to the snow pack snow_pack[t] = snow_pack[t-1] + precipitations[t] rainfall[t] = 0 return rainfall, snow_packdef module_soilmoisture(t, rainfall, sm, effective_precipitation, FC=1, beta=1): """ Soil moisture module Returns the updated soil moisture and the effective precipitation that contributes to the runoff. ------------------------------------------------------------ Parameters: t, date - Running day rainfall, array - Amount of rain and melting snow that goes to the soil sm, array - Soil Moisture water amount FC, float - Field Capacity, ie maximum storage in the subsurface zone beta, positive integer - Shape Coefficient for the contribution to effective runoff """ effective_precipitation[t] = rainfall[t] (sm[t-1]/FC)beta sm[t] = sm[t-1] + rainfall[t] (1-(sm[t-1]/FC)beta) return sm, effective_precipitationdef module_evapotranspiration(t,m, sm, ea, pea, temp,temp_m, pe_m, C=1, PWP=0): """ Evapotranspiration module Returns the evapotranspiration and the soil moisture ------------------------------------------------------------ Parameters: t, date, - Running day m, integer between 1 and 12 - Running month sm, array - Soil Moisture water amount temp, array - Expected temperatures for the running period temp_m, array - Long-term monthly mean temperatures pe_m, array - Long-term mean monthly potential evapotranspiration C, float - Model parameter PWP, float - Soil Permanent Wilting Poin """ # compute the daily adjusted potential evapotranspiration pea[t] = (1+C (temp[t]-temp_m[m])) pe_m[m] # computes actual evapotranspiration if sm[t] < PWP: ea[t] = pea[t] sm[t] / PWP else: ea[t] = pea[t] sm[t]-=ea[t] return ea, smdef module_runoff(t, effective_precipitation, rl, ru, a0, a1, a2,perc, K=[1,1,1], L=200): """ Runoff module Returns the three runoff discharges. ------------------------------------------------------------ Parameters: t, date, - Running day L, float - Upper reservoir threshold K, array - Recessions coefficients for each reservoir ru, array - Upper reservoir level rl, array - Lower reservoir level a0, array - Near surface runoff a1, array - Inflow runoff a2, array - Baseflow runoff perc, float - Quantity that is percolated from the upper to the lower reservoir effective_precipitation, array - Effective contribution to the runoff """ # inflow in reservoirs ru[t] = max(0, ru[t-1] + effective_precipitation[t] - perc) rl[t] = rl[t-1] + min(ru[t-1] + effective_precipitation[t], perc) # outflows a0[t] = K[0] max(0, ru[t]-L) a1[t] = K[1] (ru[t]-a0[t]) a2[t] = K[2] rl[t] # update reservoirs lvls ru[t] -= (a0[t] + a1[t]) rl[t] -= a2[t] return rl, ru, a0, a1, a2def total_runoff(t, a0, a1, a2): """ Returns the total runoff discharge for one day. ------------------------------------------------------------ Parameters: t, date, - Running day a0, array - Near surface runoff a1, array - Inflow runoff a2, array - Baseflow runoff """ return a0[t] + a1[t] + a2[t]def HBV(days, precipitations, temp, pe_m, temp_m, init): """ Returns the total runoff waters for each day. ------------------------------------------------------------ Parameters: days, array of dates - Days of the running period precipitations, array - Expected precipitations for the running period temp, array - Expected temperatures for the running period temp_m, array - Long-term monthly mean temperatures pe_m, array - Long-term mean monthly potential evapotranspiration sm0, float - Initial soil moisutre init, list - List of the initial values in the following order: sm0, rl0, ru0, snow_pack0 """ # Parameters: F = 3 temp_threshold = 0.0 FC = 150 PWP = 150 C = 0.05 beta = 1 K = [0.2,0.1,0.05] L = 5 Perc = 0.05 T = len(days) # initialization sm = np.zeros(T) rainfall = np.zeros(T) rl = np.zeros(T) ru = np.zeros(T) total = np.zeros(T) a0 = np.zeros(T) a1 = np.zeros(T) a2 = np.zeros(T) effective_precipitation = np.zeros(T) snow_pack = np.zeros(T) ea = np.zeros(T) pea = np.zeros(T) sm0, rl0, ru0, snow_pack0 = init sm[0] = sm0 rl[0] = rl0 ru[0] = ru0 snow_pack[0] = snow_pack0 # loop over days for t in range(1, len(days)-1): m = days[t].month rainfall, snow_pack = module_snowmelting(t, precipitations, temp, snow_pack, rainfall, F, temp_threshold) sm, effective_precipitation = module_soilmoisture(t, rainfall, sm, effective_precipitation, FC, beta) ea, sm = module_evapotranspiration(t,m, sm, temp, ea, pea, temp_m, pe_m, C, PWP) rl, ru, a0, a1, a2 = module_runoff(t, effective_precipitation, rl, ru, a0, a1, a2, Perc, K, L) total[t] = total_runoff(t, a0, a1, a2) return totalsm0, rl0, ru0, snow_pack0 = 0,0,0,0init = [sm0, rl0, ru0, snow_pack0]days = [datetime.date.today() + datetime.timedelta(days=i) for i in range(10)]pe_m = [10 for i in range(12)]temp_m = [10 for i in range(12)]precipitations = [5 for i in range(len(days)+1)]temp = [5+i for i in range(len(days)+1)]a_total = HBV(days, precipitations, temp, pe_m, temp_m, init)
参考: AghaKouchak A., Habib E., 2010, Application of a Conceptual Hydrologic Model in Teaching Hydrologic Processes, International Journal of Engineering Education, 26(4), 963-973.