黄靖雯,孟 钊,冯青郁,马伟伟,万修全,周 巍,李佳俊,冯 洋,2,3,*
1 中国科学院南海海洋研究所, 热带海洋环境国家重点实验室, 广州 510301 2 南方海洋科学与工程广东省实验室, 广州 511458 3 中国科学院南海海洋研究所,广东省海洋遥感重点实验室,广州 510301 4 中国科学院大学, 地球与行星科学学院, 北京 100049 5 中国科学院大学, 资源与环境学院, 北京 100049 6 中国科学院生态环境研究中心, 城市与区域生态国家重点实验室, 北京 100085 7 中国海洋大学, 海洋与大气学院, 青岛 266000 8 中国水产科学院,南海水产研究所, 广州 510300
河口海岸是地球岩石圈、大气圈、水圈和生物圈四大圈层交汇,物质和能量集散的重要地带,受人类陆地活动和全球气候变化的双重影响。尽管该区面积占不到全球的十分之一,但却承载了全球近三分之一的初级生产力[1]。人类陆地活动导致大量的无机氮磷营养盐排入河口海岸区域,而极端天气事件(如风暴潮、台风等过程)会对营养盐的再分布产生较大影响[2],二者复合积累,在短时间内对浮游植物生长繁殖产生剧烈影响。
我国珠江河口及邻近海域生物资源丰富、生态系统多样,承载着粤港澳大湾区最频繁,最活跃的社会经济活动,掌控着广东省乃至全国经济可持续发展命脉。近十年来,随着我国粤港澳大湾区的建设,珠江流域人口增长十分迅速,食品和能源需求急剧增加,这导致流域上游省份农业用地面积增加,下游省份城市化进程加速,随之导致含农业化肥的农田灌溉用水和生活污水排放大幅度增加,这些流域的面源污染物进入珠江河口,导致河口海域赤潮发生频率明显增高。此外,在过去的几十年间,随着全球暖化的加剧,台风活动越来越频繁。2000—2009年间,西北太平洋大概有239个台风爆发,其中79个发生在中国南海[3]。而珠江河口海岸线绵长,地势低平[4],面向大海呈高度开放格局,其独特的地理纬度和气候条件,使得南海和西太生成的台风在珠江口登陆频繁,每年给粤港澳大湾区带来巨大的经济损失[5—6]。
Qiu等[7]利用现场观测及卫星遥感数据确认了台风影响浮游植物时空分布的4个主要影响因子为:海水的垂直混合、流域的降水增加、河流淡水及营养盐的大量输入和风场变化。Zhao等[8]利用航次数据发现,这4个影响因子也是台风过境前后导致缺氧区破坏和恢复的重要因素。然而,基于观测数据的研究对机理认知还比较有限。因此,十分有必要通过数值模式来对近岸台风对生态灾害的影响深入探究。由Qiu等[7]的观测研究已知,近岸台风对浮游植物藻华的影响是一个综合流域、河口海洋、外部气象条件为一体的复杂过程。因而,模式必需集成陆面过程、海洋过程及大气过程为一体来对其进行模拟。本研究针对粤港澳大湾区,构建陆海气一体化集成模式系统(Great Bay Coupled Land-Ocean-Atmo Simulation System, GB-CLOASS),在利用观测资料对模式系统进行充分校正的基础上,一方面通过控制实验来对模式敏感性进行测试,另一方面对台风过境期间海表叶绿素时空变化机理进行探究。本研究所构建的集成模式系统对定量研究极端天气事件下河流面源污染对近岸生态灾害影响具有普适性意义,并为未来实现赤潮、缺氧等生态灾害预报提供可能性。
本研究所构建的评价体系是一个集成了陆地、海洋、大气、生源要素等多种因素,多尺度,多圈层的模式系统。其中陆面模式基于SWAT(Soil and Water Assessment Tool)构建、河口模式基于ROMS(Regional Ocean Modelling System)构建、海洋模式基于ECCO(Estimating the Circulation and Climate of the Ocean Simulation)构建、气象模式基于WRF(Weather Research and Forecasting Model)构建,不同模块之间物质或能量单向或双向交换(图1)。
图1 粤港澳大湾区陆海气集成模式示意图
1.1.1陆面模块
模式的陆面模块覆盖珠江流域,模拟了植被吸收、土壤渗透和径流输入等过程;将珠江各子流域输出的淡水、有机氮、硝氮、铵氮和无机悬浮泥沙作为驱动场单向输入河口海洋模式。
1.1.2河口海洋模块
基于ROMS构建的河口海洋模式是该评价体系的核心。该模式覆盖112.6°E—115.4°E,21.1°N—23.1°N的珠江口区域,包括珠江八大口门入海区域及珠江冲淡水影响的近岸海域,水平向为384 × 356的正交曲线网格,空间分辨率在河道内最高,约0.1 km,随离岸距离增加而减小,在河口外约3 km,垂向20层采用地形跟随坐标,每层厚度随地形而变化,其中表层和底层具有较高的分辨率。该模式在线耦合了Feng等[9]河口含碳生物地球化学循环模块,该模块包含11个生态变量:硝酸盐([NO3])、铵盐([NH4])、浮游植物(Phytoplankton,P)、浮游动物(Zooplankton,Z)、大小颗粒碎屑(Large Debris,DL和Small Debris,DS)、半易降解溶解有机氮(Semilabile Dissolved Organic Nitrogen,[DON]SL)、不易降解溶解有机氮(Refractory Dissolved Organic Nitrogen,[DON]RF)、无机悬浮物(Inorganic Suspended Solids,[ISS])、叶绿素(Chlorophyll,[Chl])和溶解氧(O2])。
其中,模式中叶绿素浓度的变化受控于浮游植物的生长与死亡过程,浮游植物生长、初级生产力增加会提高叶绿素浓度,而浮游植物代谢与死亡、浮游动物的摄食作用则会降低水体中的叶绿素浓度。此外,叶绿素浓度还因物理对流及混合过程改变,生态部分叶绿素浓度随时间变化如下:
(1)
式中,ρ[Chl]表示浮游植物生长中用于叶绿素合成的部分(无量纲);μ0为浮游植物生长速率(d-1),这里取为常数[10],受两种氮盐(LNH4和LNO3)限制和光(L1)限制;L1用以表征光合作用与光(Photosynthesis-light,P—I)的关系;LNH4和LNO3分别表示浮游植物吸收铵盐的限制和铵盐限制下吸收硝酸盐的限制(无量纲);[Chl]为叶绿素浓度(mg/m3)。E表征浮游植物代谢过程,G表征浮游动物的摄食作用,M,A和S分别表征浮游植物的死亡、聚集和沉降过程,详细的计算公式和研究中所涉及的参数具体取值参见Feng等[9]表A1和表A4。
浮游植物生长中用于叶绿素合成部分ρ[Chl]计算公式如下:
(2)
式中,θmax为叶绿素浓度与浮游植物的最大比值(mg Chl)/(mg C);P为浮游植物量(mmol N/m3);α为P—I关系线初始斜率(W-1m2d-1);I为光合作用可用辐射(W/m2)。
浮游植物铵盐限制和LNH4铵盐限制下硝酸盐限制LNO3计算公式如下:
(3)
(4)
式中,KNH4为吸收铵盐的半饱和常数(mmol N/m3),KNO3为吸收硝酸盐的半饱和常数(mmol N/m3)。NH4和NO3分别为铵盐和硝酸盐浓度(mmol/m3)。
浮游植物光合作用与光的关系L1计算公式如下:
(5)
式中, 光合作用可用辐射I,随水深z呈指数递减,计算公式如下:
I(z)=I0·PARfrac·ezKD
(6)
式中,I0为海表面之下的初始光强(W/m2); PARfrac为光强中可用于光合作用的部分(无量纲);KD为光的扩散衰减系数(m-1),在河口海域中受海表面盐度、叶绿素浓度、总悬浮物质和总溶解有机氮控制[11—12]。
1.1.3气象模式
区域气象模式向陆面模式单向输入降雨量、温度、风场、相对湿度、短波和长波辐射;同时,针对海洋与大气之间主要进行的动量、感热和潜热相互作用[13],基于COAWST(Coupled Ocean-Atmosphere-Wave-Sediment Transport)框架实现了海气耦合模式: 气象模式向海洋提供蒸发降水、海表感热、潜热通量及海面风应力和气压,海洋模式向气象模式反馈海表温度,交换频率为10 min一次。此外,大洋模式为河口模式提供开边界条件,实现河口和外部海域的物质互换。
为验证GB-CLOASS对珠江口模拟的有效性,本研究将观测站点与模式中的变量进行对比(图2),参与对比的观测量及时间跨度参见表1,其中温度和盐度代表了物理场变量,叶绿素和溶解无机氮代表了生态场变量。比较的量化指标包括均方根误差(Root Mean Square Error,RMSE)和相关系数r,计算方法如下:
表1 观测量时间
图2 模式模拟区域、模拟结果分析区域与观测站点示意图
(7)
(8)
式中,Xobs为观测值;Xmodel为模拟值;Cov(Xobs,Xmodel)表示两者协方差,Var[Xobs]和Var[Xmodel]分别表示观测与模拟值的方差。
台风“天鸽”于2017年8月19日(世界时,下同)前后在北太平洋西部(128°E,20.4°N)形成,向西移动过程中发展为热带风暴,于8月22日08时加强为强热带风暴向西北移动,22日15时强度达到台风级,并在23日07时升级为强台风[14],过境期间最大风速高达48 m/s,中心气压低至940 hPa[15]。最终台风“天鸽”于8月23日04时前后在广东省珠海市登陆,并在登陆后不久迅速减弱,于8月24日08时消散。
为了对台风过境后叶绿素变化机制进行研究,本研究设置了四组控制实验,其中,实验1包括单纯的珠江口物理-生化耦合模式,利用陆面模式SWAT模拟的气候态月平均河流径流作为输入;实验2在实验1的基础上采用SWAT模拟的日径流;实验3在实验1的基础上气象场采用WRF模拟;实验4在实验3的基础上采用SWAT日径流(表2)。本研究采用珠江口门内观测平均的温、盐及营养盐浓度作为模式的初始条件,实验1将模式从1980年1月1日开始转到2017年9月13日结束。其中,前28年为模式的加速阶段,2008—2016为模式模拟结果和观测比较的阶段,长时间段模拟足以保证台风“天鸽”过境前物理-生化耦合模式达到稳定状态。由于台风是短时间尺度内的天气事件,本研究将实验1中2017年8月15日模拟结果作为初始条件,实验2,3,4从该时间开始,分别对其进行28天的模拟。
表2 对照数值实验设置
其中,实验1与实验2、实验3与实验4采用相同的气象条件。模拟结果差异主要由径流输入差异引起;实验1与实验3、实验2与实验4径流相同、模拟结果差异主要由气象驱动条件差异引起。
珠江口以不规则半日潮为主[16],属潮优型河口[17],故集成模式对水位的模拟效果至关重要。通过对2015年5—7月赤湾、内伶仃和舢板洲三个站点处观测与模拟的海表面高度值的对比。集成模式系统的模拟结果能够体现出“朔望大潮,两弦小潮”的潮汐月不等现象以及珠江口的半日潮特征,除了对农历十五前后的大潮模拟结果较弱外,其余时间内观测与模拟的海表高度都较一致,模式几乎重现了珠江口的潮汐水位特征,高潮期间水位可达2 m左右。观测与模式之间的均方差RMSE小于0.5,相关系数r大于等于0.75(表3)。
此外,陆面模式SWAT对河道径流的模拟结果和八大口门的径流观测值表明,陆面模式的模拟结果能够再现珠江口丰枯水期的季节变化特征。观测值和模式模拟值之间的均方差RMSE在8000 m3/s左右。但是在部分月份,如2008年6月,2014年5月模式明显高估了观测的径流量。二者之间的相关系数达0.5(表3)。然而,现阶段由于施肥量的数据获取不足,SWAT模拟的氮磷浓度在入河口处和实际仍然有一定的差别,后续工作中将进一步改良。然而,由于模式采用的径流和观测具有相当可观的可比性,因此,模式模拟的进入珠江口总的面源污染物的量和实际相当。
表3 观测与模式间均方根误差(RMSE)及相关系数(r)比较
为了评估陆海气集成模式对珠江冲淡水海域物理及生态相关变量的模拟效果,对比了2014年8月、2015年1月、4月和10月月平均的观测和模拟的海表面物理场(温度、盐度)和生态场(叶绿素、无机氮)的空间分布特征,上述各海洋要素的观测时间段分别代表了珠江口夏、冬、春、秋季节性变化特征。图3中背景颜色为模式模拟值,圆点为观测值,圆点和背景颜色的融合程度反映了观测和模拟结果的一致程度。我们可以看到模式与观测的海表温度时空分布特征一致,模式再现了观测温度4个季节的变化特征。进一步地,我们将模式结果在观测的经纬度和时间点上取值,量化评估模式的模拟效果。均方差RMSE和相关系数r的计算为所有点先在每个季节进行空间排列,然后将不同季节排列好的数据组合在一起,因而模式和观测温度的RMSE在2 ℃左右,r值高达0.97(表3)。
图3 观测与模拟的海表温度、盐度、叶绿素及溶解无机氮浓度对比
在海表盐度的对比结果中,模式与观测的珠江冲淡水羽流在河口分布的形状特征基本一致,且能够体现珠江冲淡水离岸扩散较弱,外海盐度较高[18],夏季冲淡水范围最广的特点[19],但模式在夏季对珠江冲淡水的扩展有一定低估(图3)。量化比较结果显示,模式和观测盐度的RMSE在5.8 左右,r值高达0.83(表3)。
模式的模拟和观测结果均显示叶绿素浓度空间分布特征与冲淡水分布相似,体现了珠江口高浓度叶绿素夏季分布最广且主要受控于珠江径流量的河口海域特征[20—21]。而复杂的风场、流场和河口生态过程可能是导致模拟与观测存在差异性的原因。观测与模拟结果共同显示出溶解无机氮含量在春季最高,但在空间分布上,观测与模式仍存在一定差异(图3)。现阶段,河口生态模式部分是基于美国柴桑比克湾生物地球化学循环的条件构建的,而珠江口生态和柴桑比克湾有较大的不同,其中重要的一点就是珠江口大部分区域为磷限制[22],而模式中尚未考虑磷的循环过程。未来生态模式针对珠江口区域的改进,有可能大幅度改进模式模拟的叶绿素部分和营养盐部分。尽管模式模拟出的生态场和观测之间RMSE较大,但是由于生态变量的输运很大程度上被物理过程控制,因而二者取得了较好的相关性(r ≈0.9)(表3),共同体现了陆面径流对河口海洋要素分布的重要性。
由于上述模式量化RMSE和r值的计算糅合了时空分布的共同结果。为了进一步了解上述物理场和生态场观测和模式的季节变化特性。本研究进一步量化比较了模式与观测空间平均后的数理统计特性(图4)。我们发现,模拟的海表温度与观测具有高一致性,模式很好地再现了温度的季节变化,且二者之间量值相差不大。模式模拟的盐度略高于观测,但是季节趋势较为一致。海表叶绿素浓度在夏季最高,溶解无机氮浓度在春季最高。模式输出的海表面叶绿素和溶解无机氮含量与观测虽然存在一定的差异,但与观测捕捉到的季节性变化是一致的。
图4 模拟值在观测点的海表面温、盐、叶绿素、溶解无机氮对比
通过将模式输出的海表高度、径流,物理场中的温度和盐度,生态场中的海表面叶绿素浓度与溶解无机氮含量这6个变量与观测站点结果进行对比,我们发现模拟结果都较好地再现了与观测较为一致的时空分布规律,说明我们利用模式来进一步探索珠江口海表叶绿素对极端天气事件的反馈以评估陆地河流所致的面源污染对近岸水体富营养化的影响具有相当的可信性。
下面我们利用GB-CLOASS来探索河口生态环境对极端天气事件的响应。以2017年台风“天鸽”为例,WRF与ROMS耦合之后,模式模拟出的台风路径与观测显示出很好的一致性(图5),并且较好地再现了台风从8月22日由南海北部向西运移并于23日在广东省珠海市登陆的运动轨迹。
图5 模式模拟与观测的台风“天鸽”路径对比
本研究选取了WRF气象模式驱动以及SWAT陆面模式为日分辨率径流的实验4中的模拟结果与MODIS(Moderate-resolution Imaging Spectroradiometer)卫星观测数据作比较(图6)。受极端天气影响,台风“天鸽”过境期间数据缺少,过境后数据在河口区处缺失较多,因而本研究将台风过境一周(8月25日至8月31日)及过境两周(9月1日至9月10日)数据合成平均,和模式同时间段的平均做比较。近岸光学环境比较复杂,卫星反演的叶绿素一定程度上低估了实际的观测。尽管叶绿素在绝对值上和观测尚不能比较,但模拟结果和观测整体形态上显示出很强的一致性。即台风刚过境时口门外部叶绿素值偏低。而过境两周后叶绿素浓度上升的趋势,高值区面积向南海北部扩展。
图6 模拟与观测的台风“天鸽”过境后叶绿素浓度分布情况
为了更好地评估台风影响下陆地河流面源污染对近岸水体叶绿素分布的影响及模式设置给结果带来的差异,本研究对表2中的四组控制实验结果进行对比分析。首先,本研究在河口附近做区域平均,得到的时间序列如图7所示。这里实验1与实验3采用的均为SWAT陆面模式月分辨率的径流,而实验2与实验4采用的则为SWAT陆面模式日分辨率径流。实验1与实验2中的模式由观测的气象场驱动,而实验3与实验4中的模式由WRF气象模式驱动,径流分辨率驱动相同的实验模拟结果相近,相比于月分辨率下的模拟结果,采用日分辨率模拟的叶绿素浓度值在台风过境时更低,在过境后恢复与增长的幅度更剧烈,伴随着台风过境所带来的两次增水,日分辨率下盐度降低得更剧烈。不同设置下的模式结果对比说明径流是控制选择区域海表盐度及叶绿素变化的重要因素,现有的海气耦合模式工具如COAWST在近岸模拟中尚未整合陆面模块,这十分值得在未来模式耦合工具发展中考量。
图7 叶绿素和海表盐度对台风“天鸽”过境前后的响应
台风过境期间,四组实验都明显显示海表盐度迅速增加, 这是由于“天鸽”过境时向陆的风促进高盐海水进入口门内,同时抑制了低盐、高叶绿素浓度的河流冲淡水向外扩展。伴随着低盐水收缩,海表高叶绿素浓度的水也同时收缩,导致海表叶绿素浓度小幅下降。而台风过境后2—4天内,海表盐度迅速回落。其中,实验2、4的回落幅度明显大于实验1、3;这是由于实验2、4的径流由于台风带来的暴雨所致增加明显,因而盐度下降也十分明显。实验4中的叶绿素浓度在台风过境后并没有迅速回升,说明风致混合可能通过使沉积物再悬浮导致浮游植物生长受到光限制[23]。8月27日台风过境后四天,盐度出现先升高后降低的现象,升高可能是台风导致的海水近惯性运动所致,而大幅度降低则是因为台风导致的上游降雨通过径流进入珠江口具有一定的时间滞后性。海表盐度两次降低对应了实验4中叶绿素浓度的大幅度增加,叶绿素浓度在9月2号左右达到极大值,且日后仍出现再次增大趋势,这是因为珠江径流排放到河口携带的富含氮的营养盐,促使叶绿素浓度在台风过境后两周大幅度增长,持续性的增长或会促进珠江口水体的富营养化。
为了更清楚地显示台风过境前后及过境期间海表叶绿素浓度及盐度的变化,本研究选取更贴近台风实际影响下的日分辨率SWAT陆面模式驱动的实验2与实验4,分别画了两组实验在8月19日,8月23日和9月1日下海表叶绿素浓度和等盐线的空间分布图,代表台风过境前、过境期间和过境后二者的变化特征。如图8所示,台风“天鸽”过境期间,实验2与实验4中无论是叶绿素分布还是等盐线轮廓,都清楚地显示了台风对珠江冲淡水扩展的抑制作用。而在台风过境一周后即9月1日,海表叶绿素浓度最高,向口门外扩展明显。对比实验2和实验4,台风过境期间,实验4模拟的叶绿素浓度略低于实验2;而台风过境后略高于实验2。这可能是由于台风过境时强风所致的混合提高了悬浮泥沙的浓度,抑制了浮游植物生长,营养盐因未被利用而在水体中堆积。然而当台风到达松弛阶段,即过境后,前期未被利用的营养盐通过光合作用进一步促进浮游植物生长,导致叶绿素浓度在后期上升。
图8 台风“天鸽”过境前后海表面叶绿素浓度分布
认识和利用河口及其邻近海域海洋资源,保护海洋生态一直以来是我国建设海洋强国至关重要的问题。珠江河口作为粤港澳大湾区中重要一环,因特殊的地理纬度和气候条件,常年受台风侵袭,每年给粤港澳大湾区带来的经济损失惨重。在我国正将粤港澳大湾区打造成世界级城市群的时代背景下,加强台风对河口海岸生态灾害的影响机制的研究具有重大的战略和经济意义。
台风过境时的强迫阶段和松弛阶段会引发海洋水体物理和生态过程的剧烈变化[24]。一般而言,在热带及亚热带海域,例如,中国的南海和美国南部的墨西哥湾,海水温度季节性变化明显, 夏季较高的海温会促进水体分层[25],阻碍营养盐垂向输送,抑制浮游植物生长[26]。而台风过境时强风所致的混合和上升流会促进海洋的垂向交换,破坏夏季高度分层的水柱,加深混合层厚度,将低温、高营养盐的深层海水夹卷至上层海洋透光区,使海表面温度冷却[27]。而在台风过境后几天,被输送至上层海洋的营养盐通过光合作用促进浮游植物生长,导致叶绿素浓度上升[28—29]。然而,不同于南海海盆尺度的研究,近岸台风对初级生产力的影响是一个综合陆、海、气为一体的复杂过程[7—8,30]。本研究基于粤港澳大湾区综合陆地、河口、海洋及大气为一体的集成模式(GB-CLOASS),模式验证结果表明,水位、径流、温度和盐度均和观测具有良好的一致性;珠江口为典型的磷限制环境,由于模式尚未引入磷循环过程,模式模拟的叶绿素及营养盐和观测具有一定差异。以2017年8月台风“天鸽”事件为例,我们利用该集成模式研究了珠江河口海域叶绿素对极端天气事件的响应。我们发现,台风过境期间向岸的风将珠江冲淡水抑制在河口内,促进外海高盐水向口门内涌入,海表盐度增加,叶绿素浓度降低。而台风过境两周后,带给陆地的强降水通过径流流入河口,导致海表盐度大幅度下降;同时,径流携带丰富的营养盐及污染物质,加剧台风过境后近岸藻华的发生。
尽管现有的海气耦合模式工具如COAWST可以对台风影响下浮游植物动力过程进行模拟[31],然而该工具尚未包含陆面过程,因而对近岸台风影响机制研究具有很大的局限性。陆海气不同界面模式衔接,一直是全球和区域地球系统模式中具有挑战性的技术问题。目前,全球尺度地球系统模式,如Community Earth System Model(CESM)等已经在厄尔尼诺等一些大尺度海气相互作用的研究中取得了突出进展[32]。然而,现阶段大尺度的地球系统模式对陆—海相互作用的考虑较为欠缺,且模式分辨率尚不能实现对近岸海域海表盐度时空变化的有效辨认[33]。区域模式较全球尺度模式而言,陆海衔接相对成熟,如珠江口1D—3D等模式[34—35],将珠江河口海洋和流域河网模式相连接。然而,模式仅对径流进行模拟,尚未有对陆地面源污染等陆面过程的考量,且未利用统一的气象条件对陆面模式和海洋模式进行综合驱动,发展高分辨率的陆海气耦合模式是实现短时间尺度内极端天气事件对近岸物理生态过程影响的重要手段。目前针对近岸台风所涉及的陆、海、气等综合一体化过程,本研究所建立的珠江口陆海气集成模式既模拟了陆面径流输入,能够较为真实地反映陆面河流面源污染对近岸水体的影响,又充分考虑了海气相互作用,能够综合反映叶绿素浓度在台风影响下随海洋及陆地环境产生的双重变化,从而有效地评估近岸水体的富营养化情况,这为探究及预测极端天气下近岸水体富营养化,以及河口区域海洋的物理、生态响应情况提供了可能。此外,现有模式在陆面和大气衔接部分,陆面和海洋衔接部分仅采用单项输入,未来耦合器技术的引入,将最终实现陆气交换,陆海互换,并实现台风过境前后赤潮、缺氧等灾害的预警[36]。
致谢:感谢南海海洋研究所高性能计算部门提供数据支持。