天目湖流域茶園非點源污染多尺度模擬

賴正清，肖雪純，李碩*，費國松，朱立國，閆浩，夏玉林

（1.南京師范大學地理科學學院，教育部虛擬地理環境重點實驗室，南京 210023；2.南京師范大學海洋科學與工程學院，南京210023；3.江蘇省水文水資源勘測局常州分局，江蘇 常州 213000；4.江蘇省水文水資源勘測局無錫分局，江蘇 無錫 214000）

土地利用變化是非點源污染的重要影響動因，會對流域的水量與水質產生強烈的影響[1-3]?？茖W地分析土地利用變化對非點源污染的影響，對提高流域水環境質量以及實現流域生態環境持續發展具有重要的研究意義[4]。我國東南丘陵山區優質耕地缺乏，適宜建設的土地資源具有局限性，地方政府出于緩解人地矛盾的考慮，相繼出臺了引導和鼓勵丘陵山區農業開發的政策[5-6]。江蘇溧陽天目湖流域是溧陽市近80萬人口的飲用水源地和國家級旅游度假區，自21 世紀初開始，當地政府實施天目湖5A級景區的打造，在流域內部實施了以鄉村旅游和鼓勵茶果園種植生產為主的丘陵農業開發，迄今逐漸形成了集供水、旅游、農業于一體的區域發展模式[7-10]。天目湖流域丘陵山區茶園的規?；_發顯著促進了該區域的社會經濟發展，但大范圍茶園種植伴隨著的水土流失加劇、施肥量劇增形成的非點源污染對下游水源地水質安全造成極大威脅[11-12]。2015—2018 年周監測資料顯示，沙河水庫和大溪水庫共5 條入庫支流的總氮（TN）都達到劣Ⅴ類水質標準，大溪水庫水質TN 指標僅為Ⅳ類水質，而沙河水庫在2017 年和2019 年均有超過Ⅳ類水質標準的情況出現[13]。兩個水庫現狀水質與《常州市天目湖保護條例》提出的全湖Ⅱ類水質目標有較大差距。2015 年以來，地方政府逐漸認識到現有開發模式的環境負面效應，因此陸續出臺了一系列的政策，實施了退耕還林、茶果園收儲等相關措施[14]。從前期的大力支持開發，到現階段的補救措施，政策必然會對當地經濟發展產生很大程度的影響。如何對規?；鑸@開發產生的長期環境效應進行科學評價，水源地保護、茶園開發和流域生態文明建設如何協同發展，成為該區域及類似發展模式地區生態文明建設中迫切需要解決的現實問題。

得益于計算機技術和“3S”技術的發展，國內外學者將地理過程模型與遙感技術、地理信息技術、計算機技術相結合，研發出了多種流域分布式非點源污染過程模型[15-18]，這些模型成為流域非點源污染時空過程定量評估的重要工具。各類非點源污染過程模型開發的背景和應用對象不同，導致其適用性和功能存在一定差異。SWAT 模型（Soil and Water Assessment Tool model）是由美國農業部農業研究中心開發的流域尺度非點源污染模型，主要用于預測較大流域復雜多變的土壤類型、土地利用方式和管理措施條件下土地管理對水分、泥沙和化學物質的長期影響[15]。該模型在計算中采用流域-子流域-水文響應單元的空間離散方式，子流域內部的水文響應單元采用統計疊加的方式生成，因此沒有明確的地塊劃分。該模型的模擬產出是將各個水文響應單元計算結果進行統計匯總后作為子流域的單獨產出，進行河道傳輸的演算，因此在子流域內部無法進行地塊尺度的空間確定性模擬輸出。APEX 模型（Agricultural Policy/Environ?mental eXtender model）由美國德州農工大學在EPIC模型（Environmental Policy Integrated Climate model）的基礎上開發，具備對中小型流域或復雜結構農場的土地管理影響的模擬能力，可配置灌溉、排水、溝渠、緩沖帶、梯田、施肥、水庫、輪作、農藥施用、放牧和耕作等多種土地管理策略[19-20]。APEX 模型以子區即同一性質的地塊單元為最小離散單元，不同子區通過流域內的水流方向相互聯系，可從子流域層面對地塊產出進行詳盡模擬[21]，但該模型缺乏子流域內部和子流域之間河道過程的模擬，因而不具備對較大流域處理的能力。

模型集成是實現不同非點源污染模型之間的優勢互補、提高模擬能力與效率的有效手段[22-25]。因此，本研究首先利用遙感手段揭示天目湖流域土地利用的時空變化特征，解析土地利用的數量結構和空間格局，再通過集成SWAT 模型和APEX 模型，構建基于“源-匯”的一系列要素和過程的多尺度集成模擬方法，定量解析天目湖流域低山丘陵區以茶園為主要土地利用類型的非點源污染產出對水環境的影響，探討并制訂丘陵山地茶園開發影響下的水源地水生態保護調控方案，推動區域經濟和生態河湖建設的和諧發展。

1 材料與方法

1.1 研究區概況與基礎數據

天目湖流域（圖1）屬江蘇省常州市溧陽市，位于蘇浙皖三省交界，屬天目山余脈的丘陵地區，介于31°07′~31°20′N、119°21′~119°29′E 之間，流域面積約245 km2，屬北亞熱帶季風氣候，雨量充沛，四季分明，年平均氣溫16.6 ℃，最高氣溫38.7 ℃，最低氣溫-4.6 ℃，年平均降水量1 251 mm。天目湖流域主要由沙河水庫流域和大溪水庫流域構成，擁有沙河、大溪兩座大（二）型水庫，是太湖流域上游重要的水源涵養區。

圖1 天目湖流域概況Figure 1 An overview of the Lake Tianmu basin′s topography with stream networks and reservoirs

天目湖流域兩大水庫上游的低山丘陵主要來源于天目山余脈，丘陵區從東、南、西3 個方向環繞水庫，具有階梯式地貌特征。天目湖流域的地帶性土壤為黃棕壤，其呈微酸性，主要分布在低山山頂以及緩坡，保水性較差，易造成水土流失現象。為了改善當地經濟，緩解人地矛盾，20 世紀末當地政府采取了鼓勵大規模種植茶葉、果樹和以水庫為重點的旅游景點建設等政策，流域內越來越多的丘陵坡地被改造用于農業生產。大規模的茶園開發以及茶葉種植中大量的氮肥施用，導致天目湖流域兩大水庫的TN 始終處于Ⅳ類至劣Ⅴ類標準，難以達標，嚴重威脅了流域生態環境的可持續發展和居民飲用水的安全。

本研究使用的基礎數據主要包括地形、土地利用和土壤等，具體見表1。根據模型輸入要求分別構建天目湖流域氣象、土壤類型、土地利用類型和農業管理措施參數庫。

表1 模型輸入數據Table 1 Model input data

1.2 天目湖流域土地利用變化分析

對天目湖流域 2000、2005、2010、2015 年以及2019年5期土地利用類型數據進行統計，獲得各期不同土地利用類型面積、占比和變化面積，見表2。從表2 中可以看出，2000 年到2019 年以來天目湖流域的土地利用面積呈現出較大的變化。2000 年時，溧陽市政府尚未對茶、果園等經濟作物進行規?；拈_發種植，兩者僅呈零星狀分布，分布在盆地邊緣的林地占據了天目湖流域較大的土地利用面積（110.86 km2），其次為農田（104.01 km2），兩者的面積占比分別為47.10%和44.19%。2005 年開始，該流域的土地利用發生了較為顯著的變化，林地、農田開始向茶、果園轉變，林地面積減少了15.39 km2，農田更是大面積減少至55.39 km2，原來分布在水庫周邊的部分林地也已經被逐漸開發轉變為可持續帶來經濟收益的茶、果園，茶園和果園面積分別達到了9.91 km2和47.10 km2。自2005 年以后，天目湖流域內的農田、林地、草地面積呈現出不斷減少的趨勢，分別由2005 年的55.39、95.47、1.75 km2下降至 2019 年的 21.84、85.18、0.32 km2。果園和茶園的面積持續增加，果園面積由2005 年的 47.10 km2增加至 2019 年的 57.74 km2，茶園面積也由9.91 km2擴張至21.72 km2，果園以及茶園15 a 的面積擴張速度分別達到了每年近0.7 km2和0.8 km2。此外，隨著人口激增帶來的用地需求的增加，使居民區也呈現出上升的趨勢?？v觀2000—2019 年的土地利用變化數據，農田、林地呈現出下降的趨勢，面積占比分別由2000 年的44.19%和47.10%下降至2019 年的8.95%和34.90%；而果園和茶園呈現出不斷上升的趨勢，面積占比分別由2000 年的0.03%和0.16%上升至2019年的23.66%和8.90%。

表2 天目湖流域土地利用結構變化Table 2 Land use change in Lake Tianmu basin

縱觀2000—2019 年天目湖流域近20 a 茶園面積變化的統計分析，發現其時間變化特征呈現出由快速擴張向緩慢增長轉變的趨勢，開發模式由規?；_發逐漸向限制保護性開發轉變。

從空間尺度上看（圖2），2000 年時天目湖流域的茶園僅零星分布在洙漕河附近，開始規?；_發建設茶園后，茶園主要分布在大溪水庫沿湖東岸、洙漕河流域、沙河水庫沿湖西岸以及平橋河流域。2005—2019 年，由于政策的驅動和茶園高效益產出所吸引的企業投資，位于沙河西岸、洙漕河流域、平橋河流域的茶園在原有基礎上迅猛擴張。洙漕河流域的茶園面積由 2005 年的 3.13 km2增長至 2019 年的 8.03 km2，沙河西岸的茶園面積由2005 年的2.84 km2增長至2019年的4.59 km2，平橋河流域的茶園面積由2005年的1.85 km2增長至2019年的4.39 km2。這三大茶園的聚集區均為居民區，道路較為集中，可達性較高，可見茶園在天目湖流域空間分布以及變化上的差異受到了社會經濟以及基礎設施完備性的極大影響。

圖2 2000—2019年天目湖流域茶園空間分布圖Figure 2 Spatial distribution of tea plantations in Lake Tianmu basin

1.3 SWAT-APEX多尺度集成模擬

1.3.1 模型選擇與集成方法

本研究采取SWAT 模型和APEX 模型集成模擬的方法，即：在SWAT 模型劃分的子流域中選擇茶園分布的子流域，利用APEX 模型進行空間子區細分，將子流域模擬結果導入SWAT 模型，進行大流域的河道演算。通過集成方式分別發揮兩個模型的優勢，完成天目湖流域以茶園為主要土地利用類型的非點源污染產出對水環境影響的評價目標。

在本研究構建的SWAT 模型中，根據地形、土地利用類型、土壤類型和河道位置，將天目湖流域劃分為359 個子流域。為實現對流域內各茶園地塊更高的空間分辨率定量模擬，在359個子流域中篩選出82個茶園聚集的子流域，針對這些子流域進行子流域-子區的空間離散化，作為APEX模型的基本模擬單元。為保證子區劃分的科學性與合理性，子流域-子區的空間離散化對地形、土地利用以及土壤類型進行了綜合考慮，共劃分出139 個茶園子區。最終的天目湖流域SWAT-APEX模型空間離散化結果如圖3所示。

圖3 天目湖流域空間離散化結果Figure 3 Subbasin and subarea of the Lake Tianmu basin

1.3.2 模型率定與驗證

針對SWAT-APEX 集成模型的模擬需要分別完成SWAT 模型和APEX 模型相關參數的更新和校正。SWAT模型和APEX模型提供的土地覆蓋類型參數數據庫中不包含茶園類型的參數，本研究基于SWAT 模型已有作物參數庫，根據文獻在作物參數庫中添加了茶園類型的參數[26]。在此基礎上，首先根據參數的敏感性完成SWAT 模型參數的校正，參數校正結果見文獻[27]；再針對APEX 模型控制數據庫、子區數據庫、參數數據庫中相關參數，如對徑流量模擬影響較大的SCS 曲線系數、對氮磷產出模擬影響較大的氮磷富集率系數等進行率定[20，28-29]。2010—2012年作為模型預熱期，利用鰱魚橋站2013—2018年月均徑流觀測數據對模型徑流模擬效果進行評價，利用鰱魚橋站、洙漕壩橋站、潘村橋站和徐家園橋站2013—2017 年月均總氮、總磷監測數據對模型水質模擬效果進行評價。精度評價指標選取決定系數（R2）和納什系數（Ens）[28，30]，評價函數的計算公式如下：

式中：Oi和Om分別代表實測值及其平均值；Si和Sm分別表示模擬值及其平均值；n為實測或模擬數據的個數。

2 結果與討論

2.1 SWAT-APEX模型的率定與驗證結果

鰱魚橋站的月均徑流模擬和實測結果對比如圖4所示。2013—2015年為模型徑流模擬校正期，2016—2018 年為模型徑流模擬驗證期，校正期和驗證期的R2和Ens均大于0.6，表明SWAT-APEX 集成模型得到了較好的徑流模擬效果。

圖4 鰱魚橋站月均徑流量模擬值與實測值對比Figure 4 Calibration and validation results of monthly average flow from 2013 to 2018 at the Lianyuqiao station

洙漕壩橋站、潘村橋站、鰱魚橋站和徐家園橋站的月均總氮、總磷模擬和實測結果對比以及精度評價如圖5 和圖6 所示。洙漕壩橋站、潘村橋站和徐家園橋站3 個站點在校正期（2013—2015 年）和驗證期（2016—2017 年）的總氮和總磷模擬結果R2都大于0.5、Ens都大于 0.4，且其中大部分R2大于 0.6、Ens大于0.5，模擬精度滿足要求[26-27]。鰱魚橋站在2013—2015年的模擬精度偏低，但是在2016—2017年的模擬精度較高。原因是當地在2013—2015 年期間對中田舍河實施了大量河道清淤和整治工程，位于中田舍河的鰱魚橋斷面該階段的水質監測數據受到強烈的人為影響，超出了模型的模擬能力。而在后期不受河道工程影響時，模型則較好地模擬了河流水質過程?？傮w來看，SWAT-APEX集成模型獲得了較好的氮、磷模擬精度，其中總磷模擬精度整體優于總氮模擬精度。

圖5 月均總氮模擬值與實測值對比Figure 5 Calibration and validation results of monthly TN from 2013 to 2017

圖6 月均總磷模擬值與實測值對比Figure 6 Calibration and validation results of monthly TP from 2013 to 2017

2.2 茶園規模開發前后氮磷產出變化模擬與分析

基于天目湖流域SWAT-APEX 集成模型參數率定結果，在SWAT 模型中分別輸入2000、2005、2010、2015、2019 年 5 期土地利用類型數據，在 APEX 模型中分別設定不同的情景，對天目湖流域2000—2019年的總氮、總磷產出分別進行流域整體、子流域、茶園子區的多尺度分析。

在整個流域層面，根據模型模擬結果對2000、2005、2010、2015 年和 2019 年不同用地類型總氮和總磷的貢獻率進行計算，如圖7所示，從2000—2019年，茶園用地類型對天目湖流域的氮磷負荷貢獻持續上升。至2019年，茶園的年均總氮和總磷負荷貢獻率分別達到了30%和20%，是天目湖流域年均總氮負荷的第一貢獻者、年均總磷負荷的第二貢獻者。

圖7 天目湖流域不同用地類型的總氮與總磷產出占比Figure 7 Percentage of TN and TP contribution from different land use types in Lake Tianmu basin

在子流域層面，利用2019 年與2000 年天目湖流域兩期土地利用模擬的子流域多年平均總氮、總磷產量相減，獲得了2000—2019 年總氮、總磷產量變化的空間分布圖，如圖8 所示。年均總氮、總磷產量增加的子流域面積分別占55%和77%，子流域總氮產出最大增加了3 891.7 kg，子流域總磷產出最大增加了326 kg。經分析，2000—2019 年總氮、總磷產出明顯增加的子流域與其內部茶園用地類型面積的增加具有正相關關系。

圖8 天目湖流域2000—2019年總氮與總磷產量變化空間分布Figure 8 Spatial distribution of changes in TN and TP production from 2000 to 2019 in Lake Tianmu basin

在茶園子區層面，結合SWAT-APEX 集成模型的模擬結果制作出2019年茶園子區和2000年對應范圍子流域的年均總氮與總磷產出空間分布圖，如圖9 所示。2019年天目湖流域的茶園子區在2000年時的土地利用類型主要為耕地和林地。根據土地利用變化分析，近20 a中，近56%的茶園是由林地轉換而來，其余茶園則是由耕地轉換而來。在2000 年，茶園子區對應范圍子流域的年均總氮產出為0.87~710.63 kg，年均總磷產出為0.19~151.62 kg。至2019 年，茶園子區的年均總氮產出增加到3.10~3 982.81 kg，年均總磷產出增加到2.98~316.07 kg。從2000 年到2019 年，經過茶園的規?；_發，茶園子區的年均總氮與總磷產量均得到大幅增加。

圖9 2000年與2019年各茶園子區年均總氮與總磷產出空間分布Figure 9 Spatial distribution of TN and TP production from tea plantation subareas in 2000 and 2019

通過茶園規?；_發后的2019 年茶園子區總氮、總磷產出空間分布圖可以看出，茶園子區年均氮磷產出較為嚴重的區域主要聚集在洙漕河流域和平橋河流域，這與茶園集中種植區相吻合。其中，洙漕河流域既屬于大溪水庫西部的水源涵養區，也屬于十里長山-南山山水生態保護帶，平橋河流域主要屬于沙河水庫東部的水源涵養區。2000 年部分位于大溪水庫西部水源涵養區的洙漕河流域以及沙河水庫東部水源涵養區平橋河流域的水源涵養林，到2019 年時已被聚集成片的茶園所替代。茶園植株高度遠低于涵養林，且在采摘茶葉時需經常對其除草，這必然導致兩大重要水源涵養區的入湖污染攔截能力大幅降低。一旦該區域發生降雨，將產生嚴重的氮磷流失，繼而引發水庫富營養化問題，此種風險在茶園開發初期時尤為顯著。因此，就茶園子區現狀的污染程度而言，有必要對水源涵養區的茶園采取優化種植管理、設置攔污設施，甚至退茶還林等措施。

2.3 針對茶園管理措施的污染削減效果評估

根據天目湖流域不同用地類型對養分負荷的貢獻，明確了茶園是自然用地中造成氮磷污染的最重要用地類型。參考《天目湖保護規劃（2018—2035）》草案，分別從農業生產措施管理（施肥削減）、工程管理（過濾帶設置）和土地資源管理（退茶還林）3 方面因素設計3 種情景來研究不同管理措施下茶園污染產出的削減效果。情景1 設計為將茶園化肥施用量削減至原來的80%；情景2 設計為在茶園的集中種植區設置10 m的植被過濾帶；情景3設計為將所有茶園轉換為林地。各情景下茶園子區的氮磷削減效果如表3所示。

表3 3種茶園管理措施下平均總氮、總磷削減率Table 3 Average TN and TP reduction rate of three scenarios of tea plantation managements

情景1 中，各茶園子區的年均總氮、總磷產出削減率分別為6.23%~24.31%、1.75%~13.22%，平均削減率分別為14.47%和6.77%。在化肥施用量適當削減的情況下，茶園子區的年均氮磷產出削減效果并不顯著。茶園的土壤氮素呈現出輸入量高、利用率低、流失量大的特點[31]，在茶園施肥過程中也常采用高氮含量的肥料。盡管削減施肥量對茶園年均總氮、總磷產出的削減效果有限，但仍可通過適當減少施肥量、避免雨前施肥或者施用稀釋肥料等措施來減少氮磷流失。

情景2 中，各茶園子區的年均總氮、總磷產出削減率分別為25.18%~53.85%、31.58%~84.69%，平均削減率分別為29.46%和42.37%，年均總磷產量的削減效果優于總氮。整體來說，設置10 m 植被過濾帶的效果要優于控制化肥施用量的效果，然而也需要較高成本。

情景3 中，各茶園子區的年均總氮、總磷產出削減效率分別為14.57%~98.68%、17%~96.76%，平均削減效率分別為71.97%和74.42%。將茶園轉換為林地是保護流域水環境最直接、最高效的方法，然而，這種方式將大幅降低丘陵山區農民的經濟收益，與地區經濟發展相矛盾。根據前文分析，天目湖流域的水源涵養區聚集了大量茶園，是氮磷流失的高風險區。因此，可適當將位于水源涵養區的茶園轉換為涵養林，加強水源涵養區的污染物入湖攔截能力。

綜合看來，3 種情景中，退茶還林的污染削減效果最優，設置10 m 植被過濾帶的污染削減效果次之，控制施肥量的污染削減效果較為有限。然而考慮到現實情況與經濟成本，可將3 種方法適當結合，如將水源涵養區中造成高污染負荷的茶園區用林地替代，在其余茶樹種植聚集區設置5~10 m 的植被過濾帶，同時避免雨前施肥。

3 結論

（1）2000 年，茶園僅零星分布在洙漕河流域。2005—2019 年規?；_發茶園后，沙河水庫西岸、洙漕河流域、平橋河流域的茶園迅速擴張，成為天目湖流域主要的茶園聚集區。茶園開發模式呈現出水庫水源涵養區集中成片、水庫周邊向上游持續蔓延的特征。

（2）經過APEX-SWAT 集成模型的參數率定與結果驗證，模型月尺度徑流產出的確定性系數和納什效率系數在校正期達到0.6 以上，驗證期達到0.8 以上；月尺度的總氮和總磷產出的確定性系數在率定期與驗證期基本達到0.6以上，納什效率系數基本達到0.5以上。SWAT-APEX 集成模型在天目湖流域具有較好的適用性。

（3）通過模型模擬對比了茶園規?；_發前后流域總氮、總磷的時空產出特征，即：茶園用地類型對天目湖流域的氮磷負荷貢獻持續上升，至2019 年，茶園總氮產出貢獻超出其他土地利用類型成為第一貢獻者，總磷產出僅次于耕地成為第二貢獻者，洙漕河流域和平橋河流域成為兩大重點氮磷產出的茶園聚集區。

（4）設定施肥削減（將茶園化肥施用量削減至原來的80%）、過濾帶設置（在茶園的集中種植區設置10 m 植被過濾帶）和退茶還林（將所有茶園轉換為林地）3 種不同情景對茶園氮磷產出的削減效果進行評估，3 種情景的平均總氮削減率分別為14.47%、29.46%、71.97%，平均總磷削減率分別為6.77%、42.37%、74.42%。建議根據當地實際情況與經濟效益，將水源涵養區高氮磷流失的茶園區轉換為林地，結合過濾帶設置與避免雨前施肥等措施，實現丘陵山地茶園開發影響下的水源地水生態保護調控。