深层叶绿素最大值(Deep Chlorophyll Maximum, DCM)现象是全球海洋普遍存在的特征之一^[1-2], 一般出现在浮游植物水华期间稳定海洋水体的近表层和真光层底部附近。事实上, 决定海洋上层水体中浮游植物浓度大小的主要因素除了直接参与浮游植物生长所需要的营养盐与来自更高级食物链的摄食作用外, 水体的物理条件, 包括光照, 海温, 混合层深度, 三维对流和湍流扩散运输等都会影响浮游植物在垂向上的分布^[3]。在用数值模式进行模拟研究时, 模式中往往既包含代表生物化学过程的模式参数, 也包含描述物理过程和机制的模式参数。其中, 代表湍流扩散运输过程的垂向湍流扩散系数是最重要的物理参数之一。

海洋湍流混合动力过程是由海面风场, 海底摩擦, 内波破碎等触发生成的, 是海洋中动量, 质量, 能量的输送, 并通过海气相互作用影响着气候变化, 同时它也是控制海洋生态环境的关键因素。小尺度的湍流扩散对浮游植物的接触率, 摄食率, 生长率都有影响。在用数值模式进行研究时, 根据研究的问题和时间尺度, 往往有很多不同的取值办法。早在20世纪60年代, MUNK就曾指出为了维持大洋层化结构的稳定, 大洋的平均扩散系数至少需要达到10^–4 m²/s的量级。但在随后的观测实验中发现, 在远离边界的大洋区, 其扩散系数仅为10^–6~10^–5 m²/s量级。直至20世纪90年代以后, 一些观测证实了强烈的跨密度混合发生在海底边界, 特别是海槛, 海山, 海脊和海底峡谷, 其扩散系数可达10^–4 m²/s, 甚至更强^[4]。Chai等^[4]在利用一个复杂的海洋生态模式进行数值模拟研究时, 只在表层的75 m取了量级为10^–4 m²/s且递减的垂向湍流扩散系数, 75 m以下的区域为0。Fennel等^[5]取了量级为10^–5 m²/s左右的随深度变化的垂向湍流扩散系数, 进行海洋生态模式表层稳态解的研究。Liccardo等^[6]在数值模式中引入了随时间变化的垂向湍流扩散系数, 并比较了垂向湍流扩散与光的重要性。由此可见, 在用数值模式对DCM进行模拟时, 垂向湍流扩散系数的取值存在很大差异。那么垂向湍流扩散系数的不确定性到底给海洋数值模式的模拟, 尤其是DCM现象的数值模拟, 带来了怎样的影响呢？

本文将利用条件非线性最优扰动(Conditional nonlinear optimal perturbation related to parameter, CNOP)方法研究一维海洋生态模式中垂向湍流扩散系数的不确定性对DCM模拟的影响。

1 方法与模式介绍 1.1 一维海洋生态模式

本文使用的一维海洋生态模式(Nutrient-Phytoplankton model)在许多DCM的研究中都有广泛应用^[7-8]。该模式包含如下的原始方程:

$\begin{gathered} \frac{{\partial P}}{{\partial t}} = {\rm{growth - loss - sinking + mixing}} \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \;\;{\kern 1pt} {\kern 1pt} {\kern 1pt} = \frac{N}{{{H_{\rm{N}}} + N}}\frac{I}{{{H_{\rm{I}}} + I}}P - mP - \nu \frac{{\partial P}}{{\partial z}} + \kappa \frac{{{\partial ^2}P}}{{\partial {z^2}}} \\ \end{gathered}, $

(1)

$\begin{gathered} \frac{{\partial N}}{{\partial t}} = - {\rm{uptake + recycling + mixing}} \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = - {\rm{ \mathsf{ α}}} \frac{N}{{{H_{\rm{N}}} + N}}\frac{I}{{{H_{\rm{I}}} + I}}P + \varepsilon \alpha mP + \kappa \frac{{{\partial ^2}N}}{{\partial {z^2}}} \\ \end{gathered}, $

(2)

其中, $I(z, t) = {I_{\rm{in}}}{{\rm{e}}^{ - {K_{\rm{bg}}}z - k\int_0^z {P(\xi, t){\rm{d}}\xi } }}$表示入射光的光强。

垂直一维生态动力学模型不考虑水平输运对生态系统的影响, 因此我们取从海洋表面到海面下一定深度处的垂直水柱为研究对象, z和t分别表示水柱的深度和时间。P和N代表浮游植物(Phytoplankton)和营养盐(Nutrients)的浓度。κ为湍流扩散系数, 其他参数的意义和取值见表 1。

在水柱的上边界(海面: z=0)和底边界(海面下一定深度z=Z_B)处, 满足下面的边界条件:

$\left[{\kappa \frac{{\partial P}}{{\partial z}}-\nu P} \right]\left| {_{z = 0, {Z_{\rm{B}}}}} \right. = 0, $

(3)

$\frac{{\partial N}}{{\partial z}}\left| {_{z = 0}} \right., {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} N\left| {_{z = {Z_{\rm{B}}}}} \right. = {N_{\rm{B}}}, $

(4)

其中, 水柱底部的营养盐浓度用常数N_B表示。

首先我们保持各参数取值如表 1(此时湍流扩散系数在垂向上为常值), 模式的时间步长设为0.01 h, 将上面的一维海洋生态模式运行20 a。结果显示模式从第3年开始逐渐达到稳定态, 且呈现出显著的DCM特征(图 1), 50~100 m是浮游植物比较集中的区域(Productivity Layer, P-L), 叶绿素浓度取得最大值的位置位于水面下84.5 m。这与其他文献中的研究结果是类似的^[9-10]。在接下来的研究中, 我们将以此稳定态作为扰动的参考态。

1.2 条件非线性最优扰动方法(CNOP)

Mu等^[11]在2003年提出CNOP后, 又在2010年把该方法拓展到初始误差与模式参数误差同时存在的情形^[12]。其中, 与初始条件有关的参数误差称为条件非线性最优初始扰动(CNOP-I), 与模式参数有关的误差称为条件非线性最优参数扰动(CNOP-P)^[13-15]。本文中使用的是CNOP-P, 下面我们对此方法进行详细说明。

假设下面是一个以U为状态变量的动力系统:

$\left\{ {\begin{array}{*{20}{l}} {\frac{{\partial U}}{{\partial t}} = F(U, \mathit{\boldsymbol{p}}){\rm{ }}U \in {\mathit{\boldsymbol{R}}^n}, {\rm{ }}t \in \left[{{\rm{0, }}T} \right]}\\ {U\left| {_{t = 0} = {U_0}} \right.} \end{array}} \right., $

(5)

其中, F为非线性函数, t表示时间, p是参数向量, U₀为状态变量在初始时刻的值, Rⁿ表示n维向量空间。从初始时刻0到时刻T, 上述问题又可以写成:

$\left\{ {\begin{array}{*{20}{l}} {U(T;{U_{\rm{0}}}, \mathit{\boldsymbol{p}}) = {M_T}({U_{\rm{0}}}, \mathit{\boldsymbol{p}})}\\ {U(T;{U_{\rm{0}}}, \mathit{\boldsymbol{p}}) + u(T;{U_{\rm{0}}}, \mathit{\boldsymbol{p}}') = {M_{\rm{T}}}({U_{\rm{0}}}, \mathit{\boldsymbol{p}} + \mathit{\boldsymbol{p'}})} \end{array}} \right., $

(6)

其中, ${M_{\rm{T}}}$表示含有参数向量p的动力系统公式(5)0~ T时刻的传播算子, 向量p′为参数向量p的扰动, $u(T;{U_{\rm{0}}}, \mathit{\boldsymbol{p}}')$为该扰动引起的状态变量U的变化量。

给定一个范数$\left\| \bullet \right\|$, 定义如下的目标函数:

$J(\mathit{\boldsymbol{p'}}) = \left\| {{M_{\rm{T}}}({U_{\rm{0}}}, \mathit{\boldsymbol{p}} + \mathit{\boldsymbol{p'}})-{M_{\rm{T}}}({U_{\rm{0}}}, \mathit{\boldsymbol{p}})} \right\|, $

(7)

则参数向量p的扰动${\mathit{\boldsymbol{p}}_{\rm{ \mathsf{ δ} }}}$称为CNOP-P当且仅当

$J({\mathit{\boldsymbol{p}}_{\rm{ \mathsf{ δ} }}}) = \mathop {\max }\limits_{\mathit{\boldsymbol{p'}} \in \Omega } J(\mathit{\boldsymbol{p'}}).$

(8)

所以, CNOP-P即代表了在给定的扰动约束下, 引起模式模拟结果最大变化的一类参数误差。

1.3 差异进化算法(Differential Evolution)

在试验中, 我们使用目标函数:

$J(\mathit{\boldsymbol{p'}}) = \frac{1}{2}\int_{\;0}^{\;{Z_{\rm{B}}}} {P'{{(z, t)}^2}{\rm{d}}z}, $

(9)

其中, P′表示模式中浮游植物生物量的变化, p′则表示模式中参数向量的扰动。

另外, 经常用来求解CNOP-P即最优化问题公式(8)的算法, 有SPG算法^[16]、SQP算法^[17]、DE算法^[18]等等。这里我们选择了DE算法^[19], 因为DE算法是一种随机的并行直接搜索算法, 可对非线性不可微连续空间函数进行最小化, 以其易用性, 稳健性和强大的全局寻优能力在多个领域取得成功, 例如DE算法曾成功用于陆面模式中CNOP-P的计算^[19]。它的计算流程如图 2。

2 实验设计与结果分析

由于不同深度的海洋内波破碎, 涡旋, 上升流等对垂向湍流扩散系数带来的扰动是不同的, 且扰动的幅度也是未知的, 所以我们假设湍流扩散系数的扰动是随深度变化的, 为此我们把海洋上层300 m的水柱按照每50 m一个区间, 分成了6个区间, 分别记为: Λ₁=[0, 50], Λ₂=[51, 100], Λ₃=[101, 150], Λ₄=[151, 200], Λ₅=[201, 250], Λ₆=[251, 300]。

2.1 实验一:各层的不确定性比较

我们取绝对值范数, 扰动约束δ_i(i=1, 2, 3)为湍流扩散系数κ在表 1中的参考值(0.4)的20%, 40%和60%, 即: ${\delta _1} = 0.08, {\delta _2} = 0.16, {\delta _3} = 0.24$。此外, 模式中其他参数的值保持不变。且当在区间Λ_i上对湍流扩散系数κ加扰动${p_{{\Lambda _i}}}$时, 其它区间的湍流扩散系数仍保持参考值(见表 1)。则:

$J({p_{{\Lambda _i}}}) = \left\| {{M_{\rm{T}}}({U_0}, \kappa + {p_{{\Lambda _i}}}) - {M_{\rm{T}}}({U_0}, \kappa )} \right\|, $

(10)

其中, $\left| {{p_{{\Lambda _i}}}\left( z \right)} \right| \le \delta, {\kern 1pt} {\kern 1pt} z \in {\Lambda _i};{p_{{\Lambda _i}}}(z){\rm{ = }}0, {\kern 1pt} {\kern 1pt} z \notin {\Lambda _i}, {\kern 1pt} 1 \le i \le 6$,

$J\left( {p_{{\Lambda _i}}^{{\rm{CNOP}}}} \right) = \max J\left( {{p_{{\Lambda _i}}}} \right), {\kern 1pt} \;\;\;\;\;\;\;\;1 \le i \le 6.$

(11)

结果显示(图 3), 当扰动约束为参考值的20%和40%时, 使得目标函数取最大值的湍流扩散系数的扰动位于区间Λ₂(50~100 m)内, 即生产力层中DCM所在的深度以上的区域。而当振幅增加到60%时, 由区间Λ₃(100~150 m)内的湍流扩散系数的扰动导致的目标函数值迅速增加, Λ₃成为扰动对模拟结果产生最大不确定性的区间, 即生产力层中DCM以下的区域。另外, 最接近海面的区间Λ₁(0~50 m)内的湍流扩散系数的扰动对模式模拟结果的影响始终最小。相比而言, 水柱的底部边界处的扰动显然比水柱上边界的扰动对模式模拟的影响大很多。事实上, 湍流混合作为海洋混合中重要而普遍的形式, 在营养盐的垂向输送中起着重要作用, 对浮游植物的垂直分布也有着显著的影响。海面近表层(Λ₁), 风浪搅拌混合通常占有主导地位, 较强的混合作用将导致混合层内的浮游植物生物量在垂向上的分布趋于“均一”。而临近跃层附近(Λ₂)湍流扩散的强度加强, 向上输送的下层营养盐增多, 上层营养盐缺乏的水层将变浅, 导致浮游植物得以在更大范围内生长^[20]。

2.2 实验二:去掉各层的CNOP型误差

在接下来的试验中, 我们就扰动约束取参考值的20%(即δ₁=0.08)时的结果, 设计数值实验, 讨论在某个区间(${\Lambda _i}, 1 \le i \le 6$)内, 去掉一定程度的垂向湍流扩散系数的误差, 模式模拟结果的改进程度。同时, 定义湍流扩散系数在整个水柱上的不确定性引起的模式误差如下:

$J({p_{{Z_{\rm{B}}}}}) = \left\| {{M_{\rm{T}}}({U_0}, \kappa + {p_{{Z_{\rm{B}}}}}) - {M_{\rm{T}}}({U_0}, \kappa )} \right\|, $

(12)

$ 其中, \;{p_{{Z_{\rm{B}}}}}\left( z \right) \le \delta, {\kern 1pt} {\kern 1pt} z \in [0, {Z_{\rm{B}}}{\kern 1pt}], $

(13)

$J(p_{{Z_{\rm{B}}}}^{{\rm{CNOP}}}) = \max J({p_{{Z_{\rm{B}}}}}).$

(14)

另外, 我们采用下面的度量公式τ来分析减少$p_{{Z_{\rm{B}}}}^{{\rm{CNOP}}}$在某个区间(${\Lambda _i}, 1 \le i \le 6$)内一定程度的值, 模式模拟结果的改进程度:

$ \begin{array}{l} \tau = \left| {\frac{{J(p_{{Z_{\rm{B}}}}^{{\rm{CNOP}}}) - J(p')}}{{J(p_{{Z_{\rm{B}}}}^{{\rm{CNOP}}})}}} \right|, \\ p'(z) = \left\{ {\begin{array}{*{20}{c}} {\left( {1 - \mu } \right)p_{{Z_{\rm{B}}}}^{{\rm{CNOP}}}, {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} z \in {\Lambda _i}}\\ {p_{{Z_{\rm{B}}}}^{{\rm{CNOP}}}, {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} z \notin {\Lambda _i}} \end{array}} \right., {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} 1 \le i \le 6, \end{array}$

(15)

其中, 误差的减小系数μ=0.2, 0.4, 0.6, 0.8, 1.0。

从模式的输出结果(图 4)可以看出, 除了对模拟结果中浮游植物生物量影响不大的区间Λ₁外, 去掉其他区间内的不同程度的CNOP型扰动, 模式的模拟程度都得到了一定幅度的提高。当去掉CNOP型扰动所在的区间Λ₂内的误差时, 模式的改善程度最高, 最高达到了接近50%。所以, 除了接近海面的位置, 改善模式的CNOP型误差均能很好的改进模式的模拟误差, 提高模式的模拟能力, 其中, 改进生产力层上半部分, 模式误差的改进程度最高。这说明在垂向上, 区间Λ₂内的湍流扩散系数的误差对数值模式模拟DCM有着最为显著的影响。事实上, 这个位置位于混合层与温跃层交界处, 随着垂向上的湍流扩散带来的营养盐的继续补充以及合适的光照, 浮游植物在该区间内得到了较快生长, 产生了一个浮游植物相对聚集的生产力层。

2.3 实验三:去掉P-L层的CNOP型误差

浮游植物相对集中的初级生产力层(Productivity Layer, P-L)对整个水体的贡献很大, 它所含的初级生产力往往占到整个水体的30%~70%, 也具有最大的新生产力比例。而新生产力反映了海洋真光层从大气中净吸收CO₂的能力, 可见这一层对估算海洋初级生产力及认识全球碳循环都具有重要意义^[20]。下面我们在一定程度上减小实验二中得到的$p_{{Z_{\rm{B}}}}^{{\rm{CNOP}}}$在生产力层中的值, 观察模式的提高程度(表 2)。

$ \begin{array}{l} \tau = \left| {\frac{{J(p_{{Z_{\rm{B}}}}^{{\rm{CNOP}}})-J(p')}}{{J(p_{{Z_{\rm{B}}}}^{{\rm{CNOP}}})}}} \right|, \\ p'(z) = \left\{ {\begin{array}{*{20}{c}} {\left( {1-\mu } \right)p_{{Z_{\rm{B}}}}^{{\rm{CNOP}}}, {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} z \in P-L}\\ {p_{{Z_{\rm{B}}}}^{{\rm{CNOP}}}, {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} z \notin P - L} \end{array}} \right., \end{array} $

(16)

其中, $\mu = 0.2, 0.4, 0.6, 0.8, 1.0$。

通过表 2的数据可以看出, 与整个水柱上都叠加CNOP-P型扰动相比, 改善生产力层上的CNOP-P型误差确实能有效提高模式的模拟能力。当全部去掉生产力层的CNOP-P型误差后, 模式的模拟能力最大能提高了接近80%, 即只改进生产力层的垂向湍流扩散系数的误差, 就已经能够很好的改进模式的模拟结果。所以, 生产力层不但是浮游植物聚集的区域, 也是湍流扩散扰动对深层叶绿素最大值的模拟产生最大影响的区域。

3 结论与讨论

研究深层叶绿素最大值现象的分布规律和影响机制, 对提高真光层内叶绿素的含量与初级生产力的估算精度及增进对海洋碳循环过程的认识都具有重要意义, 也为校正海洋生物地球化学模式和地球系统模式参数提供重要依据^[20]。本文使用CNOP-P方法研究了一个一维海洋生态模式中, 垂向湍流扩散的不确定性对数值模式模拟误差的影响。通过计算垂向湍流扩散系数的CNOP-P型误差, 我们发现对浮游植物的生长和发展产生最大影响的垂向扩散扰动发生在生产力层的上部, 并且随着扰动幅度的增大, 逐渐过渡到生产力层的下半部分。这可能是由于此区域内对湍流扩散产生扰动的因素主要来源于海洋层结和切变。在稳定的水体中, 浮游植物的生长与营养盐的含量密切相关, 温跃层掌控着营养盐的向上补充。因此垂向湍流扩散在混合层与温跃层交接的位置(也是生产力层所在的区域)发生的扰动, 对深层叶绿素的模拟影响最大。这对我们提高DCM的模拟能力和开展进一步的研究都具有非常重要的指导意义。