本发明涉及地理学、气象学与物理学的交叉领域,特别涉及一种站点尺度的草原地区积雪物理过程模拟方法。
背景技术:
1、雪灾,又称白灾,是由冬春季集中降雪形成大面积持续的积雪覆盖,对交通、电力、通信和农牧业造成严重影响的一种主要自然灾害。
2、因此有必要对积雪物理过程进行模拟,以便预测积雪雪深及其随时间的变化过程(积雪持续时间),由此做出进一步的应对。
3、现有的积雪物理模型,如snowmodel,crocus,snowpack等雪盖模式都是基于欧洲、北美、西伯利亚等地区的观测数据建立起来的,不适用于中国的实际情况。首先,这些模型认为积雪深度超过3cm时,积雪表层反照率为常数,在0.7~0.85之间,但在中国草原区,由于草原立枯物等的存在,3cm并不能完全覆盖地表,反照率估计值偏高。其次,这些模型的新雪密度估算方法是基于欧洲、北美、西伯利亚等地区的观测数据建立起来的,应用在中国草原区,尤其是用于如青藏高原等草原区时,会出现较大的模拟偏差。这些偏差会导致无法精确估算中国草原区的积雪雪深以及积雪持续时间。
技术实现思路
1、本发明的目的在于提出一种站点尺度的草原地区积雪物理过程模拟方法,该方法在snowmodel基础上,明确了积雪层数的确定方法,并根据中国草原地区气候特点,考虑了积雪厚度对积雪反照率的影响,从而能够更为准确地估算积雪雪深以及积雪持续时间。
2、根据本发明的一种站点尺度的草原地区积雪物理过程模拟方法包括如下步骤:
3、步骤一:判断当前时间段内降水相态,
4、根据气象站点的小时时间段(例如是指1-12个小时的时间段,每个时间段例如为3小时或6小时)降水量、气温等观测数据判断是否发生降水,若发生降水,判断降水相态,即固态降水(降雪)、液态降水(降雨)或是固态液态混合降水(雨夹雪),
5、步骤二:确定当前时间段的初始地表积雪覆盖层数,
6、根据步骤一的结果以及前一时刻地表积雪覆盖层数,确定当前初始地表积雪覆盖层数,从表层积雪开始,表层积雪记为第1层,往下为第2层,以此类推,直到最底层,
7、步骤三:模拟第1层积雪雪深以及各物理量在当前时间段内的变化,
8、若第1层为新雪,则根据站点的气温、相对湿度、气压、地表温度、风速等观测数据,模拟新雪密度、新雪反照率、新雪粒径和雪深,若第1层积雪不是新雪,则模拟第1层积雪的密度、反照率、密度、雪深、积雪粒径、含水量等物理量在当前时间段内的变化,
9、步骤四:模拟第2层到最底层各层积雪雪深以及各物理量在当前时间段内的变化,
10、若存在2层或以上积雪层,根据上一层积雪的状态,模拟第2层到最底层积雪的密度、雪深、积雪密度、积雪粒径、含水量等物理量在当前时间段内的变化,
11、步骤五:根据模拟结果,确定积雪总雪深,并重新确定地表积雪覆盖层数,
12、总雪深为第1层至最底层各层积雪雪深总和,若第1层积雪雪深模拟结果为0,则地表积雪覆盖层数减1,原第2层积雪变成第1层,以此类推,否则保持不变,
13、步骤六:循环进入下一时间段积雪状态模拟,
14、根据下一时间段气象观测数据,参照步骤一到五模拟积雪物理参数变化。
15、优选地,在步骤一中,以下述方法判断当前时间段内降水相态方式为:
16、
17、prwater:降水概率,0表示没有降水,0~0.3为降雨,0.3~0.7之间为雨夹雪,大于等于0.7为降雪;
18、pw:站点观测的降水量(m);
19、ta:站点观测的气温(k);
20、tf:纯水冰点温度,273.15k;
21、rh:站点观测的相对湿度(%);
22、α,β,γ分别为拟合参数,例如采用-8.605,0.737,0.021。
23、优选地,在步骤二中,以下述方法确定当前时间段初始地表积雪覆盖层数:
24、
25、nsccurrent:当前时间段内地表积雪层数;
26、nscpre:上一时间段地表积雪层数;
27、prwater:含义同上。
28、优选地,在步骤三中,模拟第1层积雪在当前时间段内各物理量状态,以下述方法确定积雪平均温度:
29、(1)若没有积雪覆盖,即nsccurrent=0,则不模拟;
30、(2)若有积雪覆盖,即nsccurrent>0,则:
31、
32、第i层积雪平均温度(k);
33、第i层积雪顶部温度(k);
34、第i层积雪底部温度(k);
35、优选地,在步骤三中,模拟第1层积雪在当前时间段内各物理量状态,以下述方法确定积雪顶部温度和底部温度:
36、
37、
38、含义同上;
39、含义同上;
40、第1层积雪顶部温度(k),即积雪表层温度,为站点观测的地表温度;
41、最底层积雪底部温度(k),通过实测获得或者陆表模型模拟获得或者设为272.15k;
42、nsccurrent:含义同上;
43、第i层积雪深度(m)。
44、优选地,在步骤三中,模拟第1层积雪在当前时间段内各物理量状态,以下述方法确定积雪密度:
45、(1)若没有积雪覆盖,即nsccurrent=0,则不模拟;
46、(2)若有积雪覆盖,即nsccurrent>0,且pw>0,表示当前时间段内有新雪,第1层为新雪层,新雪密度为:
47、
48、ρns:新雪密度(kg/m3);
49、ta:含义同上;
50、tf:含义同上;
51、twb:湿球温度(k),可以通过ta,rh和大气压推算;
52、rh:含义同上。
53、(3)若有积雪覆盖,即nsccurrent>0,且pw=0,表示当前时间段内没有降雪,第1层为原有积雪雪层,则第1层积雪密度为:
54、
55、第i层积雪密度(kg/m3),这里i=1,为第1层积雪;
56、前一时间段第i层积雪密度(kg/m3);
57、第i层积雪压实作用导致的密度变化(kg/m3);
58、第i层积雪中液态水冻结导致的密度变化(kg/m3)。
59、优选地,在步骤三中,模拟第1层积雪在当前时间段内各物理量状态,以下述方法确定积雪压实作用导致的密度变化:
60、
61、含义同上;
62、含义同上;
63、积雪层上方积雪重量,用等效水厚度表示(m);
64、tf:含义同上;
65、含义同上;
66、n1,n2和n3是常数,例如为0.0013m-1s-1,0.08k-1和0.021m3 kg-1。
67、优选地,在步骤三中,模拟第1层积雪在当前时间段内各物理量状态,以下述方法确定积雪层上方积雪等效水厚度:
68、
69、含义同上;
70、第i层积雪的等效水厚度(m);
71、n:当前积雪层序号。
72、优选地,在步骤三中,模拟第1层积雪在当前时间段内各物理量状态,以下述方法确定第i层积雪的等效水厚度:
73、
74、含义同上;
75、含义同上;
76、第i层积雪深度(m);
77、ρw:纯水密度,1000kg/m3。
78、优选地,在步骤三中,模拟第1层积雪在当前时间段内各物理量状态,以下述方法确定积雪中液态水冻结导致的密度变化:
79、
80、含义同上;
81、qm:单位面积上用于积雪融化或冻结的能量(j/m2),大于0表示融化能量,小于0表示冻结能量;
82、第i层积雪被冻结的液态水深度(m);
83、ρw:含义同上;
84、第i层积雪前一时间段的雪深(m),这里i=1,指第1层;
85、优选地,在步骤三中,模拟第1层积雪在当前时间段内各物理量状态,以下述方法确定被冻结的液态水深度:
86、
87、含义同上;
88、qm:含义同上;
89、第i层积雪前序时刻的含水量(m);
90、lf:单位重量纯水冻结能量(j/kg),为常数,3.34×105j/kg;
91、第i层积雪最大潜在液态水冻结深度(m);
92、优选地,在步骤三中,模拟第1层积雪在当前时间段内各物理量状态,以下述方法确定最大潜在液态水冻结深度:
93、
94、含义同上;
95、ρsmax:最大积雪密度(kg/m3),一般设为550kg/m3,也可以根据研究需要自定义;
96、含义同上;
97、含义同上;
98、ρw:含义同上;
99、优选地,在步骤三中,模拟第1层积雪在当前时间段内各物理量状态,以下述方法确定积雪晶体平均直径:
100、(1)若没有积雪覆盖,即nsccurrent=0,则不模拟;
101、(2)若有积雪覆盖,即nsccurrent>0,且pw>0,表示当前时间段内有新雪,第1层为新雪层,新雪晶体平均直径确定方法为:
102、
103、gdns:新雪晶体平均直径(m);
104、g0:0.0005m;
105、a1;拟合参数,例如取0.0004m;
106、a2:拟合参数,例如取0.2s/m;
107、vwind:站点观测的风速(m/s)。
108、(3)若有积雪覆盖,即nsccurrent>0,且pw=0,表示当前时间段内没有降雪,第1层为原有积雪雪层,则第1层积雪晶体平均直径确定方法为:
109、
110、第i层积雪晶体平均直径(m);
111、第i层积雪前一时间段晶体平均直径(m);
112、g1:拟合参数,取值5.0×10-7m4 kg-1;
113、第i层积雪水蒸气通量(kg m-2s-1)。
114、优选地,在步骤三中,模拟第1层积雪在当前时间段内各物理量状态,以下述方法确定水蒸气通量确定方法为:
115、
116、uv:含义同上;
117、第i层积雪孔隙率(m3m-3);
118、第i层积雪有效扩散系数(m2s-1);
119、第i层积雪冰表面平衡蒸气压密度(kg m-3k-1),随温度变化;
120、积雪内部温度梯度(k/m)。
121、优选地,在步骤三中,模拟第1层积雪在当前时间段内各物理量状态,以下述方法确定孔隙率:
122、
123、含义同上;
124、含义同上;
125、ρi:冰密度,917kg/m3。
126、优选地,在步骤三中,模拟第1层积雪在当前时间段内各物理量状态,以下述方法确定积雪有效扩散系数:
127、
128、含义同上;
129、ee0s:1个标准大气压和0℃下有效扩散系数,9.2×10-5m2/s;
130、λa:大气压(pa),根据站点观测数据获得或者根据积雪所在位置的海拔推算;
131、含义同上;
132、tf:含义同上;
133、优选地,在步骤三中,模拟第1层积雪在当前时间段内各物理量状态,以下述方法确定冰表面平衡蒸气压密度:
134、
135、含义同上;
136、c1i:7.964×109kg k m-3;
137、lvi:冰的升华潜热,2.838×106j kg-1;
138、rw:水蒸气的气体常数,461.5j k-1kg-1;
139、含义同上;
140、优选地,在步骤三中,模拟第1层积雪在当前时间段内各物理量状态,以下述方法确定积雪表面导热系数:
141、(1)若没有积雪覆盖,即nsccurrent=0,则不模拟;
142、(2)若有积雪覆盖,即nsccurrent>0,积雪表面导热系数确定方法为:
143、(2.1)第1层积雪晶体平均直径介于0.1-0.5mm之间的:
144、
145、(2.2)第1层积雪晶体平均直径达到最大值5.0mm的:
146、
147、(2.3)第1层积雪晶体平均直径介于两者之间的,采用(2.1)和(2.3)的结果线性插值获得。
148、ks:积雪表面导热系数(w m-1k-1);
149、含义同上。
150、优选地,在步骤三中,模拟第1层积雪在当前时间段内各物理量状态,以下述方法确定积雪雪深:
151、(1)若没有积雪覆盖,即nsccurrent=0,则不模拟;
152、(2)若有积雪覆盖,即nsccurrent>0,且pw>0,表示当前时间段内有新雪,第1层为新雪层,新雪雪深确定方法为:
153、
154、含义同上,这里i=1,表示第1层积雪;
155、pw:含义同上;
156、ρns:含义同上;
157、(3)若有积雪覆盖,即nsccurrent>0,且pw=0,表示当前时间段内没有降雪,第1层为原有积雪雪层,则第1层积雪雪深确定方法为:
158、
159、含义同上,这里i=1,表示第1层积雪;
160、含义同上;
161、含义同上;
162、含义同上;
163、表层积雪前一时间段升华或凝结导致的积雪等效水厚度变化(m),负值表示升华,正值表示凝结;
164、表层积雪前一时间段消融导致的积雪等效水厚度变化(m)。
165、优选地,在步骤三中,模拟第1层积雪在当前时间段内各物理量状态,以下述方法确定表层积雪升华或凝结导致的积雪等效水厚度变化:
166、
167、ss:当前时间段内表层积雪升华或凝结导致的积雪等效水厚度变化(m),负值表示升华,正值表示凝结;
168、qe:当前时间段内积雪表面单位面积潜热交换乱流能量(j/m2);
169、含义同上;
170、含义同上;
171、lvi:含义同上。
172、优选地,在步骤三中,模拟第1层积雪在当前时间段内各物理量状态,以下述方法确定表层积雪消融导致的积雪等效水厚度变化:
173、
174、sm:当前时间段内表层积雪消融导致的积雪等效水厚度变化(m);
175、qm:含义同上;
176、含义同上;
177、含义同上;
178、lf:含义同上。
179、优选地,在步骤三中,模拟第1层积雪在当前时间段内各物理量状态,以下述方法确定第i层积雪的含水量:
180、(1)若i=1,即表层积雪:
181、
182、(2)若i>1:
183、
184、第i层积雪的含水量(m);
185、含义同上;
186、含义同上;
187、含义同上;
188、第(i-1)层积雪,即上一层积雪析出的液态水量(m)。
189、优选地,在步骤三中,模拟第1层积雪在当前时间段内各物理量状态,以下述方法确定第i层积雪析出的液态水量:
190、(1)若i=1,即表层积雪:
191、
192、(2)若i>1:
193、
194、含义同上;
195、含义同上;
196、含义同上;
197、含义同上。
198、优选地,在步骤三中,模拟第1层积雪在当前时间段内各物理量状态,以下述方法确定单位面积上用于积雪融化或冻结的能量:
199、qm=(1-α)qsi+qli+qle+qh+qe+qc
200、qm:含义同上;
201、α:当前时间段内积雪反照率(取值0-1之间);
202、qsi:当前时间段内入射到积雪表面的单位面积太阳短波辐射能量(j/m2);
203、qli:当前时间段内入射到积雪表面的单位面积大气长波辐射能量(j/m2);
204、qle:当前时间段内积雪表面单位面积出射长波辐射能量(j/m2);
205、qh:当前时间段内积雪表面单位面积显热交换乱流能量(j/m2);
206、qe:含义同上;
207、qc:当前时间段内积雪表面单位面积热传导能量(j/m2)。
208、优选地,在步骤三中,模拟第1层积雪在当前时间段内各物理量状态,以下述方法确定当前时间段内积雪反照率:
209、
210、α:含义同上;
211、αns:新雪反照率(取值0-1之间);
212、αbg:无积雪覆盖时地表反照率,草原地区取值0.2~0.3之间;
213、k:跟积雪粒径大小相关的参数;
214、积雪总深度(m);
215、t:新雪下完后到当前时间持续的天数(d);
216、a:拟合参数,例如取0.82;
217、b:拟合参数,例如取0.46。
218、优选地,在步骤三中,模拟第1层积雪在当前时间段内各物理量状态,以下述方法确定积雪总深度:
219、
220、含义同上;
221、含义同上。
222、优选地,在步骤三中,模拟第1层积雪在当前时间段内各物理量状态,以下述方法确定新雪反照率:
223、αns=c1×ln(gdns×106)+c2
224、αns:含义同上;
225、gdns:新雪晶体平均直径(m);
226、c1:拟合参数,例如取-0.059;
227、c2:拟合参数,例如取1.1132。
228、优选地,在步骤三中,模拟第1层积雪在当前时间段内各物理量状态,以下述方法确定跟积雪粒径大小相关的参数:
229、
230、k:含义同上;
231、gdns:含义同上;
232、c3:拟合参数,例如取23.939;
233、c4:拟合参数,例如取-0.629。
234、优选地,在步骤三中,模拟第1层积雪在当前时间段内各物理量状态,以下述方法确定当前时间段内入射到积雪表面的单位面积太阳短波辐射能量:
235、
236、qsi:含义同上;
237、stotal:当前时间段持续总秒数(s);
238、s*:太阳常数,1,370w/m2;
239、ψdir:直射率(0-1);
240、ψdif:漫反射率(0-1);
241、y:太阳直射光与积雪所在位置坡面夹角(rad);
242、z:太阳天顶角(rad)。
243、优选地,在步骤三中,模拟第1层积雪在当前时间段内各物理量状态,以下述方法确定太阳天顶角:
244、
245、z:含义同上;
246、积雪所在位置的纬度(rad);
247、τ:积雪所在位置和时刻的小时夹角(根据当地正午太阳测得)(rad);
248、δ:积雪所在位置和时刻的太阳高度角(rad)。
249、优选地,在步骤三中,模拟第1层积雪在当前时间段内各物理量状态,以下述方法确定积雪所在位置和时刻的小时夹角:
250、
251、τ:含义同上;
252、h:积雪所在位置的白天时长小时数(单位:小时),根据位置所在纬度和一年中的第几天查表或计算得到。
253、优选地,在步骤三中,模拟第1层积雪在当前时间段内各物理量状态,以下述方法确定积雪所在位置和时刻的太阳高度角:
254、
255、δ:含义同上;
256、0.409(rad),北回归线纬度处的弧度;
257、d:当前时间段处在一年中的第几天;
258、dr:173,北半球的夏至日;
259、dy:365.25,年均日数。
260、优选地,在步骤三中,模拟第1层积雪在当前时间段内各物理量状态,以下述方法确定太阳直射光与积雪所在位置坡面夹角:
261、cosy=cosβcosz+sinβsinzcos(μ-ξ)
262、y:含义同上;
263、z:含义同上;
264、β:积雪所在位置的坡度(rad);
265、ξ:积雪所在位置的坡向,正南方向为0;
266、μ:太阳方位角(rad),也是以正南方向为0。
267、优选地,在步骤三中,模拟第1层积雪在当前时间段内各物理量状态,以下述方法确定太阳方位角:
268、
269、μ:含义同上;
270、δ:含义同上;
271、τ:含义同上;
272、z:含义同上;
273、优选地,在步骤三中,模拟第1层积雪在当前时间段内各物理量状态,以下述方法确定直射率:
274、ψdir=(0.6-0.2cosz)(1.0-ρc)
275、ψdir:含义同上;
276、z:含义同上;
277、vc:云覆盖比例(0-1),由站点观测的云盖比例得到。
278、优选地,在步骤三中,模拟第1层积雪在当前时间段内各物理量状态,以下述方法确定漫反射率:
279、ψdif=(0.3-0.1cosz)σc
280、ψdif:含义同上;
281、z:含义同上;
282、σc:含义同上。
283、优选地,在步骤三中,模拟第1层积雪在当前时间段内各物理量状态,以下述方法确定当前时间段内入射到积雪表面的单位面积大气长波辐射能量:
284、
285、qli:含义同上;
286、stotal:含义同上;
287、ε:大气出射率(0-1);
288、σ:stefan-boltzmann常数,5.670×10-8w/m2 k4;
289、ta:含义同上。
290、优选地,在步骤三中,模拟第1层积雪在当前时间段内各物理量状态,以下述方法确定大气出射率:
291、
292、ε:含义同上;
293、ω:拟合参数,设为1.083;
294、xs:与海拔相关的拟合参数1;
295、ys:与海拔相关的拟合参数2;
296、zs:与海拔相关的拟合参数3;
297、σc:含义同上;
298、ta:含义同上;
299、ea:大气水气压(pa),采用《地面气象观测规范第6部分院空气温度和湿度观测(qx/t50-2007)》标准中的方法计算;
300、优选地,在步骤三中,模拟第1层积雪在当前时间段内各物理量状态,以下述方法确定与海拔相关的拟合参数:
301、
302、h:积雪所在位置的海拔高度(m);
303、式中i用x,y,z替代即可。x1=0.35,x2=0.51,y1=0.100k pa-1,y2=0.130k pa-1,z1=0.224,z2=1.100,h1=200m,h2=3000m。
304、优选地,在步骤三中,模拟第1层积雪在当前时间段内各物理量状态,以下述方法确定当前时间段内积雪表面单位面积出射长波辐射能量:
305、
306、qle:含义同上;
307、stotal:含义同上;
308、εs:雪表面长波发射率,设为0.98;
309、σ:含义同上;
310、含义同上。
311、优选地,在步骤三中,模拟第1层积雪在当前时间段内各物理量状态,以下述方法确定当前时间段内积雪表面单位面积显热交换乱流能量:
312、
313、qh:含义同上;
314、stotal:含义同上;
315、ρa:空气在0℃和1.0×105pa大气压下的密度,1.275kg/m3;
316、ta:含义同上;
317、含义同上。
318、cpa:0℃时空气热力系数,1004j/kg k;
319、γh:显热交换系数(m/s);
320、ξ:无量纲稳态方程或无量纲稳态系数。
321、优选地,在步骤三中,模拟第1层积雪在当前时间段内各物理量状态,以下述方法确定显热交换系数:
322、
323、γh:含义同上;
324、vwind:含义同上;
325、κ:0.4,无量纲常数;
326、zr:站点测得的风速所在高度(m),即测得vwind所在参考表面高度,一般为4m;
327、z0:动量粗糙长度,设为1.0×10-4m;
328、优选地,在步骤三中,模拟第1层积雪在当前时间段内各物理量状态,以下述方法确定无量纲稳态方程:
329、
330、ξ:含义同上;
331、ri:richardson参数,无量纲;
332、η:拟合参数,无量纲,设为9.4;
333、γ:拟合参数,无量纲;
334、优选地,在步骤三中,模拟第1层积雪在当前时间段内各物理量状态,以下述方法确定无量纲拟合参数γ:
335、
336、γ:含义同上;
337、ψ:拟合参数,无量纲,设为5.3;
338、η:含义同上;
339、γh:含义同上;
340、zr:含义同上;
341、z0:含义同上;
342、优选地,在步骤三中,模拟第1层积雪在当前时间段内各物理量状态,以下述方法确定richardson参数ri:
343、
344、ri:含义同上;
345、g:重力加速度,9.81m/s2;
346、温度梯度(k/m);
347、风速梯度(s-1);
348、优选地,在步骤三中,模拟第1层积雪在当前时间段内各物理量状态,以下述方法确定温度梯度:
349、
350、含义同上;
351、ta:含义同上;
352、含义同上;
353、含义同上;
354、zt:站点测得的气温所在高度(m),即测得ta所在参考表面高度,一般为1.5m;
355、优选地,在步骤三中,模拟第1层积雪在当前时间段内各物理量状态,以下述方法确定风速梯度:
356、
357、含义同上;
358、vwind:含义同上;
359、zr:含义同上;
360、含义同上。
361、优选地,在步骤三中,模拟第1层积雪在当前时间段内各物理量状态,以下述方法确定当前时间段内积雪表面单位面积潜热交换乱流能量:
362、
363、qe:含义同上;
364、stotal:含义同上;
365、ρa:含义同上;
366、lv:含义同上;
367、γe:潜热热交换系数(m/s),计算方法同γh;
368、ξ:含义同上;
369、ea:含义同上;
370、e0:积雪表面饱和蒸气压(pa),采用《地面气象观测规范第6部分院空气温度和湿度观测(qx/t 50-2007)》标准中的方法计算;
371、λa:含义同上。
372、优选地,在步骤三中,模拟第1层积雪在当前时间段内各物理量状态,以下述方法确定当前时间段内积雪表面单位面积热传导能量:
373、
374、qc:含义同上;
375、stotal:含义同上;
376、ks:积雪表面导热系数(w m-1k-1);
377、含义同上;
378、优选地,在步骤三中,模拟第1层积雪在当前时间段内各物理量状态,以下述方法确定积雪内部温度梯度:
379、
380、含义同上;
381、含义同上;
382、含义同上;
383、含义同上。
384、优选地,在步骤四中,模拟第2层到最底层积雪在当前时间段内各物理量状态,各层积雪平均温度确定方法同步骤三中的方法。
385、优选地,在步骤四中,模拟第2层到最底层积雪在当前时间段内各物理量状态,以下述方法确定各层积雪密度:
386、
387、含义同上,这里i>1;
388、含义同上;
389、含义同上。
390、优选地,在步骤四中,模拟第2层到最底层积雪在当前时间段内各物理量状态,各层积雪晶体平均直径确定方法同步骤三中的方法。
391、优选地,在步骤四中,模拟第2层到最底层积雪在当前时间段内各物理量状态,以下述方法确定各层积雪雪深确定方法为:
392、
393、含义同上,这里i>1;
394、含义同上;
395、含义同上;
396、含义同上。