礦渣堆積對小流域行洪能力影響的數(shù)值模擬研究

李桂伊， 侯精明， 高 波， 韓 浩， 馬利平

(1. 西安理工大學(xué)西北旱區(qū)生態(tài)水利工程國家重點(diǎn)實(shí)驗(yàn)室， 陜西 西安 710048；2.中國地質(zhì)調(diào)查局西安市地質(zhì)調(diào)查中心， 陜西 西安 710054)

1 研究背景

在礦山開采過程中會(huì)產(chǎn)生大量的廢棄礦渣，若隨意堆積在溝谷中容易產(chǎn)生滑坡、泥石流等地質(zhì)災(zāi)害，但全部運(yùn)輸出山谷則成本過高，因而礦渣的處理成為制約礦產(chǎn)業(yè)健康發(fā)展的重要問題。目前對于礦渣治理的研究更多是從防止堆積礦渣產(chǎn)生地質(zhì)災(zāi)害的角度，采用實(shí)地調(diào)查、模型實(shí)驗(yàn)、土力學(xué)模型和經(jīng)驗(yàn)公式進(jìn)行礦渣啟動(dòng)機(jī)理分析，如：徐友寧等[1]對潼關(guān)金礦區(qū)的礦渣堆數(shù)目、方量、穩(wěn)定性以及泥石流隱患進(jìn)行實(shí)地調(diào)研，并提出對應(yīng)的工程治理措施；常鳴等[2]運(yùn)用經(jīng)驗(yàn)公式從地形、水文、靜動(dòng)力及物源條件的角度定性分析泥石流運(yùn)動(dòng)模式；林玫玲等[3]采用PFC2D顆粒流軟件模擬分析了礦渣型泥石流從啟動(dòng)到暴發(fā)過程中顆粒內(nèi)部的力學(xué)特征。然而從流域產(chǎn)匯流角度研究分析粗顆粒礦渣堆積引起的地形改變對于下游流域防洪的影響的研究尚屬少見。

本文采用由西安理工大學(xué)與英國紐卡斯?fàn)柎髮W(xué)聯(lián)合開發(fā)的二維水動(dòng)力模型[4]GAST(GPU accelerated surface water flow and associated transport)建立大西岔流域雨洪徑流模型，對比分析不同重現(xiàn)期降雨條件下，礦渣堆積前后溝谷各個(gè)斷面的流量變化情況，研究礦渣堆積引起的地形變化對小流域行洪能力的影響，以期為礦渣的分類利用提供參考。

2 模型基本原理簡介

2.1 控制方程

控制方程服從二維淺水方程[8]，考慮降雨和下滲，忽略黏性項(xiàng)、風(fēng)應(yīng)力和科氏力，其矢量守恒格式如下：

(1)

(2)

式中：x、y、t分別為笛卡爾空間坐標(biāo)、時(shí)間坐標(biāo),m、s；q表征每個(gè)網(wǎng)格上的變量矢量;f和g分別為每個(gè)網(wǎng)格x、y方向上的通量矢量，包括水深h、兩個(gè)方向上的單寬流量qx和qy，m、m2/s；u、v分別表示x、y方向上的流速,m/s；S為源項(xiàng)矢量,包括摩阻力源項(xiàng)Sb和底坡源項(xiàng)Sf;i為降雨強(qiáng)度，mm/h；zb為河床底面高程,m，表征實(shí)際地形的高低起伏，作為邊界條件參與方程求解；謝才系數(shù)Cf=(gn2)/h1/3，m1/2/s，其中n為曼寧系數(shù)，是反映地形邊壁粗糙度、邊界形狀以及其他因素影響的綜合性經(jīng)驗(yàn)系數(shù)，g為重力加速度，m/s2。

2.2 數(shù)值方法

本文建立了耦合分布式水文和水動(dòng)力過程的雨洪數(shù)值模型，采用Godunov格式有限體積法[6]進(jìn)行空間離散，為穩(wěn)健地處理非連續(xù)問題，通量格式采用HLLC格式近似黎曼求解器[7]。干濕動(dòng)邊界處理在采用精度格式自適應(yīng)方法[7]的基礎(chǔ)上對動(dòng)邊界處網(wǎng)格的水深、單寬流量和底坡變量進(jìn)行靜水重構(gòu)[5]，使得在復(fù)雜地形上滿足全穩(wěn)條件。時(shí)間積分采用二階顯式Runge Kutta方法，從而構(gòu)造具有二階時(shí)間和空間精度的穩(wěn)健算法，可有效提高模型計(jì)算的精度和穩(wěn)定性。同時(shí)引入GPU并行技術(shù)[8]來加速計(jì)算，以提高模型計(jì)算效率，使得在單機(jī)上可進(jìn)行大規(guī)模計(jì)算。在文獻(xiàn)[6]中，該模型應(yīng)用于某小流域的計(jì)算結(jié)果與實(shí)測結(jié)果比對，驗(yàn)證了其計(jì)算精度。

3 模型構(gòu)建

3.1 研究區(qū)域概況

大西岔小流域位于陜西、河南兩省交接的小秦嶺金礦區(qū)西段北坡(圖1)，為我國著名的四大黃金產(chǎn)區(qū)之一[1]。大西岔溝為典型的粗顆粒礦渣堆積小流域，在地貌上，大西岔溝溝谷堆積著大量礦山開挖的廢棄礦渣，95%以上為粒徑大于1 mm的粗顆粒；水文條件上，區(qū)域內(nèi)多年平均降水量為645.8 mm，降雨量較大；地形條件上，主溝床平均縱比降為32.7%，地形陡峭。因此選取大西岔流域作為研究不同降雨條件下礦渣堆積對溝道行洪能力的影響規(guī)律具有代表性。

3.2 模型參數(shù)設(shè)置

模型邊界條件設(shè)置為自由出流的開邊界，CFL數(shù)為0.5，輸入地形、糙率、入滲和降雨數(shù)據(jù)，模擬歷時(shí)為6 h，每隔0.5 h輸出一次結(jié)果文件。

3.2.1 地形數(shù)據(jù) 研究區(qū)域流域面積為1.89 km2,地形采用實(shí)測3 m精度的DEM地形數(shù)據(jù)，分為無堆渣地形數(shù)據(jù)和有堆渣地形數(shù)據(jù)(圖2)。依據(jù)衛(wèi)星影像(圖1)將研究區(qū)域分為山坡林地和溝道碎石兩種土地利用類型，并依據(jù)相關(guān)標(biāo)準(zhǔn)及文獻(xiàn)[11]，選取林地和碎石土的曼寧系數(shù)分別為0.4、0.03，穩(wěn)定下滲率分別為6和2 mm/h。

3.2.2 降雨數(shù)據(jù) 模型設(shè)置降雨輸入采用不同重現(xiàn)期設(shè)計(jì)暴雨。根據(jù)當(dāng)?shù)厮氖謨?渭南地區(qū)暨銅川市實(shí)用水文手冊)，流域面積小于50 km2的區(qū)域設(shè)計(jì)暴雨歷時(shí)為1～3 h，本文計(jì)算了設(shè)計(jì)歷時(shí)為1～3 h的設(shè)計(jì)面暴雨量，并按照渭銅地區(qū)設(shè)計(jì)暴雨時(shí)程分配雨型表進(jìn)行暴雨時(shí)程分配，得到5種重現(xiàn)期下(5、10、20、50、100年)降雨歷時(shí)為3 h的設(shè)計(jì)暴雨過程(圖3)。

4 模擬結(jié)果對比及討論

地質(zhì)勘察報(bào)告表明：大西岔礦渣95%以上為粒徑大于1 mm的粗顆粒，不易參與泥石流過程。經(jīng)過近幾年的治理，各個(gè)礦渣堆坡度基本均小于30°，穩(wěn)定性較好。因此本文基于堆積體穩(wěn)定的假設(shè)下，模擬分析礦渣堆積引起的地形變化對大西岔流域產(chǎn)匯流過程的影響。

圖4給出了5、20、100年一遇3種典型重現(xiàn)期暴雨下大西岔流域有、無堆渣的最大水深分布情況，對比不同重現(xiàn)期的水深分布，可以看出隨著降雨量的增加，相應(yīng)的地表徑流量增加，研究區(qū)域的下游最大水深明顯增大；對比相同重現(xiàn)期原始地形與堆渣地形的水深分布，可見在3種重現(xiàn)期降雨下，4個(gè)礦渣堆積區(qū)水深均相應(yīng)有所增大，對水流產(chǎn)生一定程度的阻滯作用。

圖1大西岔流域圖圖2研究區(qū)域數(shù)字地形高程

圖3 不同重現(xiàn)期設(shè)計(jì)暴雨過程

4.1 流域出口斷面徑流量分析

如圖5所示，給出了礦渣堆積前后5種不同重現(xiàn)期的流域出口斷面流量過程。研究發(fā)現(xiàn)：不同重現(xiàn)期降雨下，礦渣堆積后洪水起漲時(shí)間延遲約1 h，洪水退水過程也有所延遲，礦渣堆積具有滯洪效果，分析認(rèn)為礦渣減小了行洪通道面積，增大了行洪阻力，從而顯示出滯洪效果。5、10、20、50、100年一遇降雨情況下，洪峰流量分別減少了5.09、7.52、8.65、8.94 、8.16 m3/s，洪峰削減率為30.94%、32.04%、28.09%、23.33%、17.30%，10年一遇降雨時(shí),洪峰削減效果最為顯著，在小于10年一遇的雨強(qiáng)時(shí)地形對于洪峰大小的影響逐漸增大，當(dāng)遇到超過10年一遇的暴雨時(shí)，地形對于洪峰大小的影響逐漸減小。5、10年一遇降雨的洪峰滯后約30 min，其他洪峰到達(dá)時(shí)間并無區(qū)別。當(dāng)雨強(qiáng)較小時(shí)，礦渣堆積引起的地形改變對洪峰到達(dá)時(shí)間起主要影響作用，隨著雨強(qiáng)增大，降雨成為影響洪峰到達(dá)時(shí)間的主要因素。分析認(rèn)為由于在溝谷堆放的礦渣束窄了行洪通道，增大了行洪阻力，對于不易參與泥石流過程的粗顆粒礦渣堆，在采取工程措施保證堆渣穩(wěn)定的前提下，堆渣的阻滯作用對徑流產(chǎn)生滯洪削峰的效果，對于減輕下游流域洪水風(fēng)險(xiǎn)具有一定積極意義。

4.2 礦渣斷面流量過程分析

從上游到下游，在長度近1 000 m的谷底分布著大量礦渣，經(jīng)過數(shù)年的工程治理[14]，主要對溝內(nèi)方量較大的礦渣堆進(jìn)行就地封存，對溝內(nèi)方量較大、堆積較高的礦渣堆進(jìn)行削坡，將堆積坡度超過30°的渣堆全部削減到坡度不超過30°，提高礦渣堆的穩(wěn)定性和抗水沖擊能力。在4個(gè)具有代表性的主要礦渣堆上依次截取4個(gè)斷面(圖1)，分析各個(gè)堆渣對流域徑流過程的影響。以50年一遇降雨為例，礦渣堆積后1、2、3、4號(hào)斷面的洪峰流量分別減少了1.77、3.00、5.36、7.51 m3/s(表1)，洪峰削減率分別為13.30%、16.30%、18.20%、20.33%，洪峰削減率隨著斷面向下游移動(dòng)逐漸增大，各個(gè)堆渣削減效果層層疊加。同時(shí)5、10、20、100年一遇降雨下的4各斷面徑流均反映以上規(guī)律。分析認(rèn)為從上游到下游堆渣數(shù)目越來越多，對水流的阻滯作用逐步累積，相應(yīng)的洪峰削減效果逐步累加。

圖4 典型重現(xiàn)期暴雨時(shí)研究區(qū)域最大水深分布圖

圖5 礦渣堆積前后5種不同重現(xiàn)期下流域出口斷面流量過程

表1 各個(gè)堆渣斷面洪峰峰值情況 m3/s

注：Qa為無堆渣時(shí)斷面峰值流量，Qb為有堆渣時(shí)斷面峰值流量。

5 結(jié) 論

本文采用基于水動(dòng)力學(xué)的雨洪數(shù)值模型，以礦渣堆積典型小流域大西岔為研究對象，計(jì)算了5種不同重現(xiàn)期暴雨下，礦渣堆積引起的地形改變對流域行洪能力的影響，得出以下結(jié)論：

(1)礦渣堆積后，流域出口洪水起漲時(shí)間和洪水退水過程均有所延遲，顯示出一定的滯洪效果。

(2)5、10、20、50、100年一遇降雨情況下，流域出口斷面洪峰削減率分別為30.94%、32.04%、28.09%、23.33%、17.30%，10年一遇降雨時(shí)洪峰削減效果最為顯著。

(3)當(dāng)小于10年一遇降雨時(shí)，地形對洪峰大小、到達(dá)時(shí)間起主要影響作用，反之，降雨成為主要影響因素。

(4)洪峰削減率隨著斷面向下游移動(dòng)而逐步增大，各個(gè)堆渣對洪峰削減效果層層疊加。

鑒于粗顆粒礦渣堆積引起的地形改變對下游流域能夠產(chǎn)生滯洪削峰效果，對減輕下游流域洪水風(fēng)險(xiǎn)具有一定積極意義。這一結(jié)論給類似的礦渣治理工程提供了相應(yīng)的借鑒，將不易參與泥石流的粗顆粒礦渣堆積在溝谷，并采用削坡、鉛絲籠加固等工程措施保證堆渣體穩(wěn)定[15]，可起到一定的滯洪削峰作用。后續(xù)將對不同溝道比降的流域中的不同礦渣堆堆放位置、不同暴雨雨型及重現(xiàn)期的滯洪削峰效果展開進(jìn)一步研究。

參考文獻(xiàn)：

[1] 徐友寧,陳社斌,何 芳,等. 潼關(guān)金礦區(qū)礦渣型泥石流災(zāi)害及防治對策[J]. 山地學(xué)報(bào),2006,24(6):667-671.

[2] 常 鳴,唐 川. 基于水動(dòng)力的典型礦山泥石流運(yùn)動(dòng)模式研究[J]. 水利學(xué)報(bào),2014,45(11):1318-1326.

[3] 林玫玲,簡文彬,胡海瑞,等. 基于離散元的礦渣泥石流運(yùn)動(dòng)過程研究[J]. 中國地質(zhì)災(zāi)害與防治學(xué)報(bào),2017,28(2):10-14+120.

[4] HOU Jingming,LIANG Qiuhua,SIMONS F,et al. A stable 2D unstructured shallow flow model for simulations of wetting and drying over rough terrains[J]. Computers & Fluids,2013,82(17):132-147.

[5] HOU Jingming,LIANG Qiuhua,ZHANG Hongbing,et al. An efficient unstructured MUSCL scheme for solving the 2D shallow water equations[J]. Environmental Modelling & Software,2015,66(C):131-152.

[6] HOU Jingming,?ZGEN I. A model for overland flow and associated processes within the hydroinformatics modelling system[J]. Journal of Hydroinformatics,2014,16(2):375-391.

[7] HOU Jingming,SIMONS F,MAHGOUB M,et al. A robust well-balanced model on unstructured grids for shallow water flows with wetting and drying over complex topography[J]. Computer Methods in Applied Mechanics & Engineering,2013,257(15):126-149.

[8] LIANG Qiuhua,XIA Xilin,HOU Jingming. Catchment-scale high-resolution flash flood simulation using the GPU-based technology[J]. Procedia Engineering,2016,154:975-981.

[9] SMITH L S,LIANG Qiuhua. Towards a generalised GPU/CPU shallow-flow modelling tool[J]. Computers & Fluids,2013,88(12):334-343.

[10] LIANG Qiuhua,MARCHE F. Numerical resolution of well-balanced shallow water equations with complex source terms[J]. Advances in Water Resources,2009,32(6):873-884.

[11] ENGMAN E T. Roughness coefficients for routing surface runoff[J]. Journal of Irrigation & Drainage Engineering,1986,112(1):39-53.

[12] 李貴玉. 黃土丘陵區(qū)不同土地利用類型下土壤入滲性能對比研究[D]. 楊凌: 西北農(nóng)林科技大學(xué). 2007.

[13] 渭南地區(qū)水利水土保持局. 渭南地區(qū)暨銅川市實(shí)用水文手冊[M].渭南: 渭南地區(qū)暨銅川水利局水土保持局,1973.

[14] 方茜娟,任軍旗,郭書科. 礦山泥石流形成條件綜述[J]. 巖土工程界,2008，11(5):66-67.

[15] 孫新國,彭 勇,周惠成. 基于SWAT分布式流域水文模型的下墊面變化和水利工程對徑流影響分析[J]. 水資源與水工程學(xué)報(bào),2016,27(1):33-39.

[16] 黃 沛,張秋文. 幾種典型流域水文模型類比分析[J]. 水資源與水工程學(xué)報(bào),2006,17(5):27-30.