郭云彤，周 妍，崔亞莉，邵景力

（中國地質(zhì)大學(xué)(北京)水資源與環(huán)境學(xué)院, 北京 100083）

我國西北干旱區(qū)水資源匱乏，隨著社會(huì)經(jīng)濟(jì)的發(fā)展，區(qū)域用水需求增大，生產(chǎn)生活用水不斷擠占生態(tài)用水，地下水水位下降導(dǎo)致植被退化、濕地面積衰減、生態(tài)系統(tǒng)退化等一系列生態(tài)問題[1]。為緩解區(qū)域各部門用水矛盾，遏制生態(tài)環(huán)境惡化的趨勢，近年來生態(tài)用水相關(guān)研究已成為當(dāng)前西北干旱區(qū)水資源開發(fā)利用的重點(diǎn)。

青土湖是我國西北河西走廊三大流域之一—石羊河流域的尾閭湖[2]，作為騰格里沙漠和巴丹吉林沙漠之間的生態(tài)屏障，青土湖濕地在防止沙漠合攏、遏制流域生態(tài)惡化趨勢上具有重要作用。由于氣候變化及人類活動(dòng)的影響，青土湖于1959年完全干涸，濕地消失[3]；自2010年開始，紅崖山水庫有計(jì)劃地沿渠道向下游進(jìn)行生態(tài)輸水，青土湖開始形成季節(jié)性水面，周邊生態(tài)環(huán)境明顯好轉(zhuǎn)[4]。因此，研究青土湖生態(tài)輸水量與湖水面面積的關(guān)系，對于確定合理的生態(tài)輸水量尤為重要。石羊河的地表來水及調(diào)水工程的水匯入紅崖山水庫后，通過渠道供給下游民勤盆地用水（主要為農(nóng)業(yè)灌溉用水）和青土湖生態(tài)輸水（圖1）。當(dāng)前針對生態(tài)輸水對青土湖生態(tài)環(huán)境的影響已有一些研究成果[3,5-7]，然而這些研究多是結(jié)合遙感解譯和定位觀測的方法，對獲得的數(shù)據(jù)進(jìn)行統(tǒng)計(jì)分析得到輸水量-湖水面積的規(guī)律，如果用這個(gè)統(tǒng)計(jì)規(guī)律外推，幅度偏大，其結(jié)果會(huì)存在著很大的不確定性。

圖1 紅崖山水庫、民勤盆地與青土湖Fig.1 Location of the Hongyashan Reservoir, Minqin Basin and Qingtu Lake

地表水-地下水耦合數(shù)值模擬方法是定量分析輸水與湖區(qū)面積關(guān)系最有效的方法。當(dāng)前地表水-地下水的耦合模擬模型大致分為獨(dú)立型、聯(lián)合型和集成型3 類。獨(dú)立型是在相對成熟的地下水模型的基礎(chǔ)上拓展部分地表水模擬的功能，如模塊化三維有限差分地下水流動(dòng)模型（Modular Three-dimensional Finite-difference Ground-water Flow Model，MODFLOW）的湖泊（LAK）、河流（RIV）、蒸散發(fā)（EVP、EVT）和入滲補(bǔ)給（RCH）模塊等，這些模塊已經(jīng)在實(shí)際地下水模擬研究中被普遍應(yīng)用。聯(lián)合型是將成熟的地表水和地下水模型通過一定手段進(jìn)行連接，共同構(gòu)建地表水-地下水的模擬系統(tǒng)，模型中各子系統(tǒng)按一定順序進(jìn)行獨(dú)立計(jì)算，而子系統(tǒng)之間只進(jìn)行單向或雙向傳輸[8-9]。集成型模型是一種完全耦合模型，該類模型將地表水和地下水作為一個(gè)系統(tǒng)，通過同時(shí)求解各個(gè)水文過程的控制方程描述地表水和地下水水分的交換和動(dòng)態(tài)變化過程，每一步都有其具體的物理意義，機(jī)理性強(qiáng)[10-11]。獨(dú)立型模型雖加入了部分地表水模擬功能，但對大區(qū)域地表水流動(dòng)過程的刻畫仍不足；聯(lián)合型模型較獨(dú)立型模型機(jī)理性更強(qiáng)的同時(shí)，又不需要集成型模型中過多的參數(shù)，是當(dāng)前地表水-地下水耦合模擬中最為實(shí)用的模型。

由美國地質(zhì)調(diào)查局2008年公布的用于模擬地表水-地下水相互作用的模型（Coupled Groundwater and Surface-Water Flow Model，GSFLOW）是近年來應(yīng)用較多的聯(lián)合型耦合模型，該模型將降雨徑流模型系統(tǒng)（Precipitation Runoff Modeling System，PRMS）和MODFLOW 進(jìn)行耦合，能夠同時(shí)模擬氣候、地表徑流、地下潛流以及融雪、湖泊、溪流和濕地等與地下水之間的相互作用，已在國內(nèi)外部分地區(qū)的地表水-地下水相互作用研究中得以成功應(yīng)用[12-14]。當(dāng)前GSFLOW 的模擬技術(shù)多應(yīng)用于流域尺度的地表水和地下水關(guān)系的模擬中，多研究河流與地下水的相互作用。GSFLOW 很好地將MODFLOW 的湖泊（LAK）模塊和地表水模塊進(jìn)行耦合，適用于本次針對青土湖的模擬研究。本次研究利用地表水-地下水耦合數(shù)值模擬的方法，定量化研究了流域中上游生態(tài)輸水對下游尾閭湖水域面積的影響。

本文基于高分辨率的DEM 數(shù)據(jù)、利用ArcGIS 分析功能確定了水面面積與水位、蓄水量的轉(zhuǎn)換關(guān)系；進(jìn)而運(yùn)用GSFLOW 軟件構(gòu)建地表水-地下水耦合數(shù)值模型，通過模擬得到不同生態(tài)輸水量情景下青土湖水面面積，論證青土湖的適宜生態(tài)輸水量。

1 數(shù)據(jù)來源及分析

1.1 數(shù)據(jù)來源

本次研究針對青土湖生態(tài)輸水量-湖水面積的關(guān)系，以數(shù)值模型為基礎(chǔ)，其中涉及的數(shù)據(jù)主要包括鉆孔資料、補(bǔ)給項(xiàng)（降雨、生態(tài)輸水?dāng)?shù)據(jù)）、排泄項(xiàng)（蒸發(fā)數(shù)據(jù)）以及模型校準(zhǔn)所需的驗(yàn)證數(shù)據(jù)（水位數(shù)據(jù)）等。具體數(shù)據(jù)如下：

（1）研究區(qū)不同分辨率的高程數(shù)據(jù)，包括90 m×90 m 的DEM 數(shù)據(jù)及1 m×1 m 的無人機(jī)遙感地形數(shù)據(jù)[15]，低分辨率數(shù)據(jù)作為數(shù)值模型的頂板標(biāo)高，高分辨率數(shù)據(jù)用于后續(xù)確定不同輸水量下水面面積的實(shí)際分布。

（2）研究區(qū)2010—2019年多期Lansat5 及l(fā)ansat8 及國產(chǎn)高分2 號（GF-2）遙感影像數(shù)據(jù)，用于識別植被及水體面積。參考NDVI 及MNDWI 進(jìn)行分類[16-17]，利用ArcGIS 對其分布進(jìn)行矢量化處理，獲得研究區(qū)年內(nèi)較為連續(xù)的湖面及蘆葦分布及面積信息[18]。

（3）本次研究主要采用民勤站2010—2019年日尺度的蒸發(fā)和降水?dāng)?shù)據(jù)，氣象數(shù)據(jù)來自中國氣象數(shù)據(jù)網(wǎng)的中國地面國際交換站氣候資料日值數(shù)據(jù)集（data.cma.cn），經(jīng)處理后主要作為數(shù)值模型中的源匯項(xiàng)。

（4）生態(tài)輸水?dāng)?shù)據(jù)來自民勤縣水務(wù)局2010—2019年的統(tǒng)計(jì)資料；地下水水位數(shù)據(jù)來自位于青土湖的地下水水位觀測井的記錄值。

（5）含水介質(zhì)參數(shù)和水文地質(zhì)參數(shù)初值根據(jù)研究區(qū)地層巖性及前人研究成果確定[19]。

綜合以上資料，青土湖水面面積和蘆葦面積的年內(nèi)變化特征見圖2。從數(shù)據(jù)較完整的2013—2018年看，水面面積增加趨勢與輸水時(shí)段基本重合，隨著輸水的進(jìn)行，水面面積持續(xù)增大，直至輸水結(jié)束，在水面蒸發(fā)的作用下水面面積隨即減小。蘆葦面積逐年上升，且隨著多年生態(tài)輸水的進(jìn)行，蘆葦面積趨于穩(wěn)定。綜上，青土湖區(qū)水面面積與蘆葦面積的年內(nèi)變化特征為：水面面積變化在生態(tài)輸水期主要受生態(tài)輸水影響，非生態(tài)輸水期主要受蒸發(fā)影響；生態(tài)輸水后的水面面積增大對蘆葦面積增大有促進(jìn)作用。

圖2 水體、蘆葦面積及入湖水量變化Fig.2 Variations in the surface water area, reed area and water inflow into the lake

1.2 水面面積與水位、蓄水量的轉(zhuǎn)換方法

建立水面面積與水位、蓄水量的轉(zhuǎn)換關(guān)系是后續(xù)濕地地表-地下水耦合模擬驗(yàn)證的重要參考依據(jù)，轉(zhuǎn)換方法具體如下：

（1）利用高分辨率的DEM 數(shù)據(jù)識別湖泊區(qū)域；

（2）利用ArcGIS 的Surface Analysis 工具以0.1 m 的水位差為單位，計(jì)算不同湖水位時(shí)對應(yīng)的淹沒面積及淹沒面以下的體積；

（3）對得到的水位與其對應(yīng)的面積及體積數(shù)據(jù)進(jìn)行趨勢分析，即可得到水面面積與水位、蓄水量的轉(zhuǎn)換關(guān)系公式。

遙感解譯結(jié)果為湖面面積的變化，需要建立水面面積與水位、蓄水量的轉(zhuǎn)換關(guān)系，并將湖面面積轉(zhuǎn)換為湖水水位，作為驗(yàn)證湖泊模擬結(jié)果合理性的驗(yàn)證數(shù)據(jù)。水面面積與水位、蓄水量的轉(zhuǎn)換關(guān)系見圖3。

圖3 青土湖水面面積與水位關(guān)系Fig.3 Relationship between the surface water area and water level of the Qingtu Lake

如圖3 所示，水位與水面面積的關(guān)系式可表示為：

式中：S—水面面積/km2；

hl—湖水水位/m。

式（1）中R2分別達(dá)0.996 7 和0.999 9。

2 數(shù)值模擬方法

本次研究采用GSFLOW 進(jìn)行青土湖地下水-地表水的耦合模擬。GSFLOW 將MODFLOW 中使用的有限差分單元和PRMS 中使用的水文響應(yīng)單元（Hydrological Response Unit，HRU）進(jìn)行空間鏈接是PRMS 和MODFLOW 模型耦合中的一個(gè)關(guān)鍵步驟。這個(gè)過程是通過生成重力儲(chǔ)層（GRV）實(shí)現(xiàn)的，GRV 的作用是進(jìn)行HRU 和有限差分單元間的水量傳輸。由于HRU 和有限差分單元具有不同的空間范圍，所以每個(gè)GRV的空間范圍由HRU 和有限差分單元的交叉部分確定。每一個(gè)GRV 被分配一個(gè)唯一的標(biāo)識號，每個(gè)URU和有限差分單元內(nèi)可以有多個(gè)GRV。如圖4 中A 部分所示，第一層為MODFLOW 中的有限差分單元，第二層為PRMS 的HRU，第三層顯示了GRV 分布。具體關(guān)系如圖4 中B 部分所示，包括4 個(gè) 水文響應(yīng)單元、6 個(gè)有限差分單元和9 個(gè)GRV。在實(shí)際應(yīng)用中，需要為每個(gè)GRV 指定拓?fù)鋮?shù)，將每個(gè)GRV 與相應(yīng)的HRU 和有限差分單元聯(lián)系。每個(gè)GRV 的重力排水通過GSFLOW 模塊gsflow_prms2mf 被添加到對應(yīng)的有限差分單元中。類似地，每個(gè)有限差分單元的地下水排出量也可通過GSFLOW 模塊gsflow_mf2prms 被添加到GRV 中。

圖4 PRMS 與MODFLOW 耦合模式[20]Fig.4 Coupling model between PRMS and MODFLOW[20]

利用GSFLOW 中湖泊（LAK）模擬功能對湖泊隨生態(tài)輸水的變化以及與地下水的交互進(jìn)行模擬，GSFLOW 中的湖泊在PRMS 中表示為湖泊HRU，在MODFLOW-2005 中表示為一組有限差分單元。在MODFLOW中直接在湖泊模塊中輸入作用于湖泊的降水量、蒸發(fā)量和地表徑流數(shù)據(jù)[21]。然而在GSFLOW 中，這些過程以及通過土壤帶的壤中流是在PRMS 中計(jì)算的，因此LAK 模塊的這些輸入變量應(yīng)設(shè)置為0。流入湖泊HRU 的地表徑流及壤中流計(jì)算公式為：

式中：—第m時(shí)間步長第n次迭代中從貢獻(xiàn)HRU 到湖泊HRU 的地表徑流和壤中流體積/m3；

FJ,lakeHRU—HRUJ 中為湖泊HRU 提供地表徑流和壤中流的面積占總面積的十進(jìn)制分?jǐn)?shù)，在GSFLOW 的參數(shù)hru_pct_up中進(jìn)行定義；

J—HRU 的計(jì)數(shù)；

JJ—某一段河段貢獻(xiàn)地表徑流和壤中流的HRU 的總數(shù)。

進(jìn)入MODFLOW 中湖泊的量表示為：

式中：—第m時(shí)間步長第n次迭代中從湖泊HRU 到MODFLOW-2005定義的湖的體積流量/（m3·d-1）；

C'prms2mf—prms 的單位到modflow 中單位的換算系數(shù)；

—第m時(shí)間步長的湖泊HRU 上的降雨量/m；

—第m時(shí)間步長的湖泊HRU 上的蒸發(fā)量/m；

AlakeHRU—湖泊HRU 的面積/m2。

假設(shè)湖床存在于湖泊單元和含水層有限差分單元之間，并且具有與底層有限差分單元不同的特性，則湖泊和地下水之間交換量的計(jì)算公式為：

式中：—第m時(shí)間步長第n次迭代中穿過湖床至含水層有限差分單元中心的體積流量/（m3·d-1）；Klkbd—湖床的水力傳導(dǎo)系數(shù)/（m·d-1）；

thicklkbd—湖床厚度/m；

thickaq—含水層厚度/m；

—湖床覆蓋有限差分單元的面積/m2；

Kaq—靠近湖單元的含水層有限差分單元的水平或垂向滲透系數(shù)/（m·d-1）；

—第m時(shí)間步長第n次迭代中的湖水位/m；

—第m時(shí)間步長第n次迭代中有限差分單元靠近單元節(jié)點(diǎn)處湖泊的地下水水頭/m。

3 結(jié)果與分析

3.1 模型建立與識別

依據(jù)研究區(qū)的地質(zhì)和水文地質(zhì)條件等，確定本次模型模擬范圍包括青土湖及其周邊部分沙漠地區(qū)，南部以隱伏斷層為界，設(shè)定為隔水的零流量邊界；東北側(cè)以沙漠邊緣為界，設(shè)定為通用水頭邊界；考慮到裸土區(qū)地下水的極限蒸發(fā)埋深，以地下水埋深2.5 m 為界，西北部以等水位線1 309 m 和1 311.5 m 為界，西南部以1 307.5 m 為界，設(shè)定為定水頭邊界；總面積約159 km2。模擬區(qū)的補(bǔ)給項(xiàng)主要為降水補(bǔ)給和人工生態(tài)輸水補(bǔ)給，其中上游紅崖山水庫對青土湖的人工輸水是整個(gè)區(qū)域的主要補(bǔ)給來源。模型的主要排泄項(xiàng)為湖水及地下水的蒸散發(fā)。模擬期為2015年1月1日—2019年12月31日，1 d 為一個(gè)應(yīng)力期進(jìn)行模擬。

地下水模型區(qū)剖分為500 m×500 m 的網(wǎng)格，如圖5（a）。研究區(qū)地表標(biāo)高采用90 m×90 m 的DEM 數(shù)據(jù)及1 m×1 m 的無人機(jī)遙感地形數(shù)據(jù)拼接作為研究區(qū)的地表高程使用。

圖5 模型范圍與空間離散Fig.5 Simulation area and spatial discretization

含水層結(jié)構(gòu)如圖6 所示，湖區(qū)西部、東南部及南部均有黏粒含量高的沉積物分布，透水性差，東北部主要以沙土為主。為刻畫含水層水文地質(zhì)條件，反映湖水與地下水間的補(bǔ)排特征，將整個(gè)含水系統(tǒng)劃分為一層，含水層厚度定為15 m。

圖6 湖區(qū)鉆孔分布及地層剖面[19]Fig.6 Borehole distribution and stratigraphic profile of the Qingtu Lake[19]

地表水模型部分，對研究區(qū)的地表空間、地下空間及模型接口進(jìn)行劃分。研究區(qū)內(nèi)地形相對平坦，降水稀少，地表產(chǎn)匯流較少，因此將研究區(qū)整體劃分為陸地HRU 和湖泊HRU，圖5（b）；利用GSFLOW 模型的GRV 接口將地表HUR 和地下水模型的有限差分單元網(wǎng)格連接起來，共劃分769 個(gè)GRV，圖5（c）。

利用民勤縣氣象站監(jiān)測的2015—2019年降雨量數(shù)據(jù)、氣溫?cái)?shù)據(jù)作PRMS 的數(shù)據(jù)文件，對研究區(qū)的降雨入滲和蒸散發(fā)進(jìn)行計(jì)算；利用輸水渠道模塊SFR2 模擬向青土湖的渠道輸水，利用湖泊模塊LAK模擬湖泊的變化。

青土湖水位變化擬合效果見圖7。模型模擬水位變化與遙感解譯水位變化趨勢基本一致，模型模擬結(jié)果基本上能夠反映湖水實(shí)際的水位變化情況，納什指數(shù)NSE 為0.766、決定系數(shù)R2為0.772。然而，由于青土湖區(qū)實(shí)際的地形相對平坦，但也存在微小的起伏，網(wǎng)格剖分的空間分辨率不足以反映某些微小的地形起伏變化，因此模型模擬的湖泊水位與實(shí)際水位存在一定的偏差，特別是當(dāng)水位較低時(shí)，偏差相對較大。

圖7 青土湖模擬水位與實(shí)際水位變化Fig.7 Changes in the simulated water level and actual water level of the Qingtu Lake

模擬區(qū)觀測孔位置如圖8 所示。為了驗(yàn)證模型模擬效果，2018年和2019年觀測的實(shí)際水位與模擬水位進(jìn)行對比（圖9），結(jié)果表明，觀測孔的模擬水位變化與實(shí)測值擬合較好。

圖8 觀測孔分布圖Fig.8 Distribution of the observation wells

圖9 觀測孔擬合圖Fig.9 Fitting of the predicted water level and measured water level

3.2 生態(tài)輸水方案預(yù)測與結(jié)果分析

生態(tài)輸水量是青土湖湖面變化最主要的影響因素，為研究生態(tài)輸水量對湖面面積變化及地下水水位變化的影響，分別設(shè)置不同生態(tài)輸水量：2 000×104～4 500×104m3/a 每隔100×104m3/a 設(shè)置1 個(gè)輸水方案，4 500×104～6 000×104m3/a 每隔500×104m3/a 設(shè)置1 個(gè)輸水方案。輸水時(shí)段均在8—10月，共設(shè)置29 個(gè)輸水方案。運(yùn)用已建立的地表水-地下水流耦合模型，預(yù)測未來20年模擬區(qū)內(nèi)地表水和地下水的變化情況。模型的結(jié)構(gòu)、水文地質(zhì)參數(shù)及邊界條件均保持不變，模型模擬期設(shè)置為2020年1月1日—2040年12月31日。模型的降雨蒸發(fā)項(xiàng)需輸入逐日數(shù)據(jù)。本次研究搜集到的逐日降雨蒸發(fā)數(shù)據(jù)為2010—2019年。從蒸發(fā)數(shù)據(jù)看，2000年后研究區(qū)年蒸發(fā)量在1 876～2 924 mm范圍內(nèi)波動(dòng)（圖10）。故在預(yù)測模型中降雨蒸發(fā)項(xiàng)取對應(yīng)2010—2019年每日的多年平均值，年總蒸發(fā)量約2 300 mm。

圖10 民勤盆地1960—2019年降水、蒸發(fā)量變化Fig.10 Variation of rainfall and evaporation from 1960 to 2019

不同輸水量情景下預(yù)測湖水水位變化情況見圖11。隨著生態(tài)輸水量的增大，水位變化呈現(xiàn)由大到小的過程。預(yù)測初期5～6年，湖水水位快速增高，變化幅度與生態(tài)輸水量的大小有關(guān)；6年后湖水水位緩慢變化。當(dāng)生態(tài)輸水量為2 000×104m3/a 時(shí)，湖水水位逐漸下降，最高水位保持在1 210.4 m（平均水位1 210.0 m），對應(yīng)水面面積6.68 km2；當(dāng)年輸水量保持在3 100×104m3/a時(shí)，與現(xiàn)狀模型2019年相比，水位變化不大，最高水位保持在約1 212.2 m（平均水位1 211.7 m），說明現(xiàn)狀3 100×104m3/a 的生態(tài)輸水量可以保證青土湖維持當(dāng)前水面面積的生態(tài)需水。隨著輸水量的增大，到2040年，生態(tài)輸水量達(dá)6 000×104m3/a 時(shí)，湖泊最高水位可達(dá)1 217 m（平均水位1 215.9 m），年內(nèi)變化幅度達(dá)2.45 m。

圖11 不同輸水方案下湖水水位變化圖Fig.11 Water level variations of the lake under different water conveyance schemes

由于各觀測孔水位變化趨勢相近，故以V01 孔為例（圖12）顯示不同方案下地下水水位的變化情況，該位置地表高程約為1 211 m，見圖12。與湖泊水位變化相似，隨著生態(tài)輸水量的增大，預(yù)測初期5～6年，地下水水位快速增高， 6年后地下水水位緩慢變化。當(dāng)輸水量保持在3 100×104m3/a 時(shí)，地下水水位的年際變化不大，與現(xiàn)狀年相比基本保持穩(wěn)定，最高水位維持在約1 210.8 m，最低水位維持在約1 210.6 m；隨著生態(tài)輸水量的增大，湖泊對地下水的補(bǔ)給量逐漸增大，地下水的穩(wěn)定水位逐漸提升，當(dāng)水量達(dá)4 000×104m3/a 時(shí)，地下水水位達(dá)到地表高程；當(dāng)輸水量達(dá)5 000×104m3/a 時(shí)，水位常年高于地表以上；當(dāng)輸水量達(dá)6 000×104m3/a 時(shí)，最高水位可達(dá)1 211.3 m，高出地表面約0.3 m。

圖12 不同方案下V01 觀測孔地下水水位變化圖Fig.12 Water level variations of observation well v01 under different water conveyance schemes

生態(tài)輸水量分別為3 100×104，4 500×104，6 000×104m3/a 時(shí)，預(yù)測2039年湖泊水均衡情況見表1。

如表1 所示，當(dāng)生態(tài)輸水量小于4 500×104m3/a時(shí)，隨著生態(tài)輸水量的增加，水面面積不斷增大，水面蒸發(fā)量增大幅度明顯，是湖泊的主要排泄項(xiàng)。當(dāng)生態(tài)輸水量大于4 500×104m3/a 時(shí)，隨著湖泊水位的提高，地下水和地表水位差增大，湖泊向地下水的排泄量增大，湖面面積增大幅度很小，水面蒸發(fā)量增量有限。

表1 不同生態(tài)輸水方案下2039年預(yù)測湖泊水均衡情況Table 1 Predicted lake water balance in 2039 under different ecological water conveyance schemes /104 m3

3.3 青土湖生態(tài)輸水適宜量確定

隨著生態(tài)輸水量的增大，湖面面積也逐漸增大，但增大的幅度有所變化。將水位轉(zhuǎn)化成湖水的水面面積，則生態(tài)輸水量與水面面積及水面面積變化率的關(guān)系見圖13。隨著生態(tài)輸水量的增大，水面面積增大情況大致可以分為3 個(gè)階段：當(dāng)輸水量為2 000×104～3 700×104m3/a 時(shí)，水面面積隨輸水量增大而增大，基本為線性關(guān)系，面積變化率相對穩(wěn)定，保持在8 122～10 796 m2/（104m3·a）；輸水量為3 700×104～4 500×104m3/a時(shí)，水面面積增大幅度減緩，面積變化率逐漸大幅度減小至2 000 m2/（104m3·a）；當(dāng)輸水量大于4 500×104m3/a時(shí)，隨生態(tài)輸水的增大水面面積增大幅度很小，特別是當(dāng)生態(tài)輸水量大于5 500×104m/a 時(shí)，面積變化率趨近于0，水面面積幾乎穩(wěn)定在26 km2。從維持當(dāng)前生態(tài)水面面積不至減小的前提看，生態(tài)輸水量不宜低于3 100×104m3/a。另外從生態(tài)輸水的效益考慮，生態(tài)輸水不宜高于4 500×104m3/a。因此，青土湖生態(tài)輸水適宜量為3 100×104～4 500×104m3/a，維持湖面面積大致為16.27～26.60 km2。

圖13 不同輸水量條件下最大水面面積變化圖Fig.13 Maximum surface water area under different water conveyance schemes

針對生態(tài)輸水與青土湖水面面積的變化研究已有一些研究成果，如利用統(tǒng)計(jì)分析和水均衡分析方法對2010年生態(tài)輸水以來的生態(tài)輸水?dāng)?shù)據(jù)、氣象數(shù)據(jù)和水面面積變化數(shù)據(jù)進(jìn)行分析，結(jié)果顯示生態(tài)輸水量在3 145×104m3/a 時(shí)可以基本維持當(dāng)前的最大水面面積13.43 km2[22-23]。而本次研究的結(jié)果顯示生態(tài)輸水量3 100×104m3/a 時(shí)水面面積可維持16.27 km2，結(jié)果與統(tǒng)計(jì)分析及水均衡分結(jié)果相比，相同生態(tài)輸水狀況下，維持的水面面積偏大。這可能是由于數(shù)值模擬得到的水面面積是生態(tài)輸水多年累計(jì)達(dá)到的較為穩(wěn)定的水面面積，略大于短期生態(tài)輸水形成的水面面積。

4 結(jié)論

（1）依據(jù)不同生態(tài)輸水方案的模擬結(jié)果，當(dāng)前3 100×104m3/a 的生態(tài)輸水量可以保證青土湖維持當(dāng)前水面面積的生態(tài)需水。

（2）生態(tài)輸水量與水面面積變化的關(guān)系大致可分為3 個(gè)階段。結(jié)合該關(guān)系，同時(shí)考慮向青土湖生態(tài)輸水的效益及保證當(dāng)前水面面積的需求，紅崖山水庫向青土湖的生態(tài)輸水適宜量為3 100×104～4 500×104m3/a，湖面面積大致可維持在16.27～26.60 km2。

致謝：感謝清華大學(xué)胡宏昌教授及首都師范大學(xué)李浩乾碩士在地表水-地下水耦合數(shù)值模擬技術(shù)方面提供的支持，感謝馬瑞教授、陳喜教授共享青土湖地形數(shù)據(jù)、地下水水位數(shù)據(jù)，感謝民勤縣水務(wù)局提供的詳細(xì)生態(tài)輸水?dāng)?shù)據(jù)！