全球变暖对海洋浮游生物生长–产氧的影响与应对策略

Effects of Global Warming on the Growth-Oxygen Production of Marine Plankton and Coping Strategies

摘要: 本文研究在全球变暖的背景下，海洋浮游生物生长–产氧情况的变化趋势。基于浮游动植物与氧气的关系，用函数刻画相关物质流的循环与传递，得到恒温状态下三者的基础模型，接下来引入温度上升对其中各化学反应的影响，同时在原有模型的基础上加上一个空间扩散维度，并考虑生态系统对气候的负反馈作用，最终得到浮游生物–氧气–气温变化模型。用有限差分法将微分方程进行离散，并结合python编程对浮游植物、浮游动物、氧气浓度和上升温度这四者的关系进行求解，绘制数量变化曲线。模拟结果表明，近500年内，浮游生物与氧气的波动将趋于平缓，振幅逐渐减小，且整体有略微下降的趋势；与此同时，温度上升的幅度也趋于平缓，浮游生物与气温二者互相制约，使生态系统更加稳定。

Abstract: We consider the changing trends of marine plankton growth-oxygen production in the context of global warming in this paper. Based on the relationship between phytoplankton and oxygen, the cycle and transport of related material flows are characterized by functions, and the basic model of the three is obtained under constant temperature, and then the effect of temperature rise on each of the chemical reactions is introduced, and a spatial diffusion dimension is added on the basis of the original model, and the negative feedback effect of the ecosystem on the climate is considered, and the plankton-oxygen-temperature change model is finally obtained. Using the differential differentiation method of differential equations and combined with python programming, the relationship between phytoplankton, zooplankton, oxygen concentration and rising temperature is solved, and the quantitative change curve is plotted. The simulation results show that in the past 500 years, the fluctuations of plankton and oxygen will tend to be flattened, the amplitude will gradually decrease, and the overall trend will be slightly decreasing; at the same time, the amplitude of temperature rise will also tend to be flat, and the two will restrict each other, making the ecosystem more stable.

1. 引言

随着经济全球化的发展，人类对生态环境的影响也越来越显著，其中以气候变化对海洋生物多样性的改变最为深刻且影响面较广，人们也愈加关注这方面的环境问题。海洋浮游生物作为海洋生态系统中的生产者、消费者，对海洋生态系统的各项循环起着至关重要的调控作用，与此同时，也受到海洋环境(如海流、温度、营养盐等)的影响，二者之间具有强烈的相互制约关系 [1] [2]。

而全球变暖带来的海洋温度上升，一方面促使藻类释放出更多的二甲基硫(DMS) [3]，形成大量云层并带来酸雨 [4]，使温度降低的同时制约藻类等浮游生物的生长；另一方面，主要温室气体二氧化碳含量的升高将强化“海洋表层混合层”的分层 [5]，从而阻止水面生长的浮游生物吸收关键营养物质，最终导致其光合作用下降。

海洋的浮游生物供给了地球一大部分的氧气 [6]，因而在海洋生态环境问题 [7] 逐渐严峻的今天，对氧气–浮游生物模型和全球变暖影响问题的研究也越来越迫切 [8]。

2. 模型的建立

我们需要建立模型，量化分析全球变暖对海洋浮游生物生长和产氧速率的影响关系，故问题的重点在于浮游生物(包含浮游动物和浮游植物)与氧气浓度的关系以及全球变暖带来海洋温度升高对这种关系的影响程度 [9] [10]。

2.1. 浮游生物与氧气的制约与循环

浮游生物根据生态系统中的角色分类，有浮游动物和浮游植物两种。浮游植物光合作用导致氧气浓度增加，这些产生的氧气在浮游植物和浮游动物的呼吸、代谢活动中被消耗，同时浮游动物捕食浮游植物，控制其生长密度。三者关系 [11] 如图1所示。

Figure 1. Schematic diagram of the relationship between plankton and oxygen

图1. 浮游生物与氧气关系示意图

图1中c表示海洋中的氧气浓度，u表示浮游植物密度，v表示浮游动物密度，箭头表示在该系统中相关物质流的传递，其速率如箭头上所附参数所示。

其中 $f\left(c\right)$ 表示溶解氧从浮游植物细胞输送到周围水体的浓度增加速率，参数A量化了细胞内叶绿体的产氧率，其主要与阳光利用率、二氧化碳利用率和水温等因素有关，总的来看， $Af\left(c\right)$ 描述了每单位质量浮游植物的产氧率； $u_R$ 、 $v_R$ 分别表示浮游植物与浮游动物呼吸作用消耗氧气的速率； $mc$ 表示水中氧气自然消耗的情形(比如一些生化反应)，其中m为该消耗的速率； $g\left(c,u\right)$ 表示浮游植物的平均生长率，我们认为其与浮游植物的密度与水体含氧量有关，故定义为关于二者的函数； $p\left(u,v\right)$ 表示浮游动物对浮游植物的捕食速率，与二者的密度密切相关； $\sigma u$ 表示浮游植物自然死亡率， $\mu v$ 表示浮游动物自然死亡率，二者均与自身密度呈正比关系。

由此，海洋中氧气的增长率可以表示为：浮游植物产氧率减去浮游生物呼吸作用消耗氧气效率，再减去水体自身氧气耗损率 [12]，即：

$\frac{\text{d}c}{\text{d}t}=Af\left(c\right)u-u_R\left(c,u\right)-v_R\left(c,v\right)-mc$ (1)

浮游植物的增长率 [12] 可以表示为自身平均生长速率减去被浮游动物捕食的速率，再减去自然状态下的死亡速率，即：

$\frac{\text{d}u}{\text{d}t}=g\left(c,u\right)u-p\left(u,v\right)-\sigma u$ (2)

浮游动物的增长率 [12] 可以表示为自身平均生长速率减去自然状态下的死亡速率，并且此时自身增长可以看作是对浮游植物捕食的一种转化，效率为 $k\left(c\right)$，即：

$\frac{\text{d}v}{\text{d}t}=k\left(c\right)p\left(u,v\right)-\mu v$ (3)

最终三者的关系可以表示为：

$\{\begin{array}{l}\frac{\text{d}c}{\text{d}t}=Af\left(c\right)u-u_R\left(c,u\right)-v_R\left(c,v\right)-mc\hfill \\ \frac{\text{d}u}{\text{d}t}=g\left(c,u\right)u-p\left(u,v\right)-\sigma u\hfill \\ \frac{\text{d}v}{\text{d}t}=k\left(c\right)p\left(u,v\right)-\mu v\hfill \end{array}$ (4)

由文献 [13] 可知，公式(4)中相关函数如表1所示。

Table 1. Correlation function expressions

表1. 相关函数表达式

其中 $k_i\left(i=1,2,3,4\right)$ 为半饱和常数；B为浮游植物线性增长率； $\beta$ 为浮游动物捕食率， $\eta$ 为浮游动物最大摄食效率； $\gamma$ 为当种群密度较大时出现的内部竞争速率。需要注意的是，所有系数均是无量纲的。

于是公式(4)可以转化为：

$\{\begin{array}{l}\frac{\text{d}c}{\text{d}t}=A\left(1-\frac{c}{c-k_0}\right)u-\frac{uc}{c+k_2}-\frac{\eta cu}{c+k_3}-c\hfill \\ \frac{\text{d}u}{\text{d}t}=\left(\frac{Bc}{c+k_1}-\gamma u\right)u-\frac{\beta uv}{u+h}-\sigma u\hfill \\ \frac{\text{d}v}{\text{d}t}=\frac{\eta c^2}{c^2+k_4^2}\cdot \frac{\beta uv}{u+h}-\mu v\hfill \end{array}$ (5)

公式(5)为温度恒定情况下的基础模型。

2.2. 全球变暖对模型的影响 [14] [15] [16]

基础模型揭示了氧气、浮游动物、浮游植物间的相互制约、转化关系，其中涉及到的化学反应如光合作用、呼吸作用等均与环境温度有着密切的联系 [17]，在全球变暖的背景下，需对该模型进行进一步明确、改进。

前文曾提及：参数A量化了细胞内叶绿体的产氧率，其主要与阳光利用率、二氧化碳利用率和水温等因素有关，在这里将其看做光合作用的剧烈程度，又根据已有的资料显示，温度每上升10℃，光合作用的强度 [18] 提升至原来的2倍，即：

$A=A_0\cdot 2^{\frac{\Delta T}{10}}$ (6)

其中 $A_0$ 为初始的光合作用的强度， $\Delta T$ 为上涨的温度。

同时，对于浮游植物的自然增长率 $g\left(c,u\right)$，我们认为其会受到温度的影响且这种影响是线性的，即在其前面乘上一个函数 ${\alpha }_1\left(T\right)$ ：

${\alpha }_1\left(T\right)=1-{\rho }_1\Delta T$ (7)

同理，对于浮游生物的呼吸作用耗氧量 $u_R\left(c,u\right)$ 、 $v_R\left(c,v\right)$，我们同样认为其受温度的线性影响，即：

$u_R\left(c,u\right)\to {\alpha }_2\left(T\right)u_R\left(c,u\right)$ (8)

$v_R\left(c,v\right)\to {\alpha }_3\left(T\right)v_R\left(c,v\right)$ (9)

其中，两类浮游生物所受线性影响视作一致，即：

${\alpha }_2\left(T\right)={\alpha }_3\left(T\right)=1+{\rho }_2\Delta T$ (10)

代入公式(5)，得到全球变暖背景下，海洋浮游生物的生长–产氧模型：

$\{\begin{array}{l}\frac{\text{d}c}{\text{d}t}=A_0\cdot 2^{\frac{\Delta T}{10}}\left(1-\frac{c}{c-k_0}\right)u-\left(1+{\rho }_2\Delta T\right)\frac{uc}{c+k_2}-\left(1+{\rho }_2\Delta T\right)\frac{\eta cu}{c+k_3}-c\hfill \\ \frac{\text{d}u}{\text{d}t}=\left(1-{\rho }_1\Delta T\right)\left(\frac{Bc}{c+k_1}-\gamma u\right)u-\frac{\beta uv}{u+h}-\sigma u\hfill \\ \frac{\text{d}v}{\text{d}t}=\frac{\eta c^2}{c^2+k_4^2}\cdot \frac{\beta uv}{u+h}-\mu v\hfill \end{array}$ (11)

2.3. 增加移动维度

在实际情况中，浮游生物可以选择离开原有位置以摆脱灭绝 [19]，氧气也有着由高浓度向低浓度扩散的趋势，与此同时，即使是对于同一片海域，由于水层并没有充分混合，各个点的浮游生物密度、初始含氧量、营养物质等量均有差异，即在空间上存在梯度，现对该梯度的存在予以考虑。

为简化问题，仅考虑一维情形，即仅沿着x方向变化。

2.3.1. 氧气

此时氧气浓度c不仅与时间有关，还与空间x有关，故等式左端写作 $\frac{\partial c\left(x,t\right)}{\partial t}$，同时在等式右端加上扩散项 ${\nabla }^{2}c$，得到现有的关于氧气浓度c的微分方程：

$\frac{\partial c\left(x,t\right)}{\partial t}={\nabla }^{2}c+A\left(1-\frac{c}{c-k_0}\right)u-{\alpha }_2\left(T\right)\frac{uc}{c+k_2}-{\alpha }_3\left(T\right)\frac{\eta cu}{c+k_3}-c$ (12)

其中：

$A=A_0\cdot 2^{\frac{\Delta T}{10}}$ (13)

${\alpha }_2\left(T\right)={\alpha }_3\left(T\right)=1+{\rho }_2\Delta T$ (14)

2.3.2. 浮游植物

同理得到现有的关于氧气浓度u的微分方程：

$\frac{\partial u\left(x,t\right)}{\partial t}={\nabla }^{2}u+{\alpha }_1\left(T\right)\left(\frac{Bc}{c+k_1}-\gamma u\right)u-\frac{\beta uv}{u+h}-\sigma u$ (15)

其中：

${\alpha }_1\left(T\right)=1-{\rho }_1\Delta T$ (16)

2.3.3. 浮游动物

同理得到现有的关于氧气浓度v的微分方程：

$\frac{\partial v\left(x,t\right)}{\partial t}={\nabla }^{2}v+\frac{\eta c^2}{c^2+k_4^2}\cdot \frac{\beta uv}{u+h}-\mu v$ (17)

2.3.4. 增加空间梯度后的氧气–浮游生物模型

此时的模型变为：

$\{\begin{array}{l}\frac{\partial c\left(x,t\right)}{\partial t}={\nabla }^{2}c+A\left(1-\frac{c}{c-k_0}\right)u-{\alpha }_2\left(T\right)\frac{uc}{c+k_2}-{\alpha }_3\left(T\right)\frac{\eta cu}{c+k_3}-c\hfill \\ \frac{\partial u\left(x,t\right)}{\partial t}={\nabla }^{2}u+{\alpha }_1\left(T\right)\left(\frac{Bc}{c+k_1}-\gamma u\right)u-\frac{\beta uv}{u+h}-\sigma u\hfill \\ \frac{\partial v\left(x,t\right)}{\partial t}={\nabla }^{2}v+\frac{\eta c^2}{c^2+k_4^2}\cdot \frac{\beta uv}{u+h}-\mu v\hfill \end{array}$ (18)

其中：

$\{\begin{array}{l}A=A_0\cdot 2^{\frac{\Delta T}{10}}\hfill \\ {\alpha }_1\left(T\right)=1-{\rho }_1\Delta T\hfill \\ {\alpha }_2\left(T\right)={\alpha }_3\left(T\right)=1+{\rho }_2\Delta T\hfill \end{array}$ (19)

此时浮游生物可选择在 方向上移动、迁徙以摆脱灭绝，氧气也会由高浓度向低浓度渗透。

2.4. 浮游生物对气温的反馈

2.4.1. DMS对气温上升的影响

浮游生物对于气温的影响在于，当气候变暖引起海洋温度上升，某些藻类的数量因此剧烈增长而在全体浮游植物中的占比增大，这些藻类释放出更多的二甲基硫(DMS)，促进产生大量云层，起到减少太阳热辐射、降低气温的作用。

因而二甲基硫的产量可用下式简化表示：

$P\_DMS=\alpha u\cdot k$ (20)

其中 $\alpha$ 表示产生二甲基硫的浮游植物在整体浮游植物中的占比，k表示这些浮游植物的单位产量。

对于 $\alpha$，假定其增长速度可用下式表示：

$\alpha =1-{\text{e}}^{\frac{\tau }{\Delta T}}$ (21)

即随着温度的上升，该类浮游植物的占比逐渐增大且无限接近于1。

我们认为二甲基硫的升温速度的影响是线性的，即：

$\omega ={\omega }_0-\gamma P\_DMS$ (22)

其中 $\omega$ 表示经过浮游植物调节后的温度上升速率， ${\omega }_0$ 表示初始的温度上升速率，即地球变暖速率。

则温度T可用下式表示：

$T={\displaystyle {\int }_{t_0}^t\omega \text{d}\tau }+T_0$ (23)

其中 $T_0$ 表示初始温度。

2.4.2. 浮游生物–氧气–气温变化模型

对公式(20)~(23)进行整理化简，得到温度T的最终表达式为：

$T={\displaystyle {\int }_{t_0}^t{\omega }_0-\left(1-\frac{\tau }{{\text{e}}^{\Delta T}}\right)}\cdot \gamma u\cdot k\text{d}\tau +T_0$ (24)

结合公式(18)，可得到最终的表达式为：

$\{\begin{array}{l}\frac{\partial c\left(x,t\right)}{\partial t}={\nabla }^{2}c+A\left(1-\frac{c}{c-k_0}\right)u-{\alpha }_2\left(T\right)\frac{uc}{c+k_2}-{\alpha }_3\left(T\right)\frac{\eta cu}{c+k_3}-c\hfill \\ \frac{\partial u\left(x,t\right)}{\partial t}={\nabla }^{2}u+{\alpha }_1\left(T\right)\left(\frac{Bc}{c+k_1}-\gamma u\right)u-\frac{\beta uv}{u+h}-\sigma u\hfill \\ \frac{\partial v\left(x,t\right)}{\partial t}={\nabla }^{2}v+\frac{\eta c^2}{c^2+k_4^2}\cdot \frac{\beta uv}{u+h}-\mu v\hfill \\ T={\displaystyle {\int }_{t_0}^t{\omega }_0-\left(1-\frac{\tau }{{\text{e}}^{\Delta T}}\right)}\cdot \gamma u\cdot k\text{d}\tau +T_0\hfill \end{array}$ (25)

其中：

$\{\begin{array}{l}A=A_0\cdot 2^{\frac{\Delta T}{10}}\hfill \\ {\alpha }_1\left(T\right)=1-{\rho }_1\Delta T\hfill \\ {\alpha }_2\left(T\right)={\alpha }_3\left(T\right)=1+{\rho }_2\Delta T\hfill \end{array}$ (26)

3. 模型的求解

3.1. 有限差分法

对于模型求解，我们采用简单显格式的有限差分法 [20]。

将公式(20)~(23)所示微分方程差分化，得到如下差分方程：

$\{\begin{array}{l}P\_DMS\left[i\right]=averag{e}_j\left[u\left(i,j\right)\right]\cdot k\cdot \left(1-{\text{e}}^{\frac{\tau }{\Delta T}}\right)\hfill \\ \omega \left[i\right]={\omega }_0-\gamma P\_DMS\left[i\right]\hfill \\ T\left[i+1\right]=T\left[i\right]+\omega \left[i\right]\cdot \Delta t\hfill \\ \Delta T=T\left[i\right]-T\left[0\right]\hfill \end{array}$ (27)

对公式(18)进行差分化，得：

$\{\begin{array}{l}{c}_{\left(t+1,x\right)}=\Delta t\left[c_k\left(t+1,x\right)+A\frac{k_0}{c+k_0}-{\alpha }_1\frac{u_n{c}_n}{c_n+k_2}-{\alpha }_2\frac{v_n{c}_n{u}_n}{c_n+k_3}-c_n\right]+c_n\hfill \\ u_{\left(t+1,x\right)}=\Delta t\left[u_k\left(t+1,x\right)+\left(\frac{{\alpha }_{3}B{c}_n}{c_n+k_1}-\gamma u_n\right)\cdot u_n-\frac{\beta v_n{u}_n}{u_n+h}-\sigma u_n\right]+u_n\hfill \\ v_{\left(t+1,x\right)}=\Delta t\left[v_k\left(t+1,x\right)+\frac{\eta c_n^2}{c_n^2+k_4^2}\cdot \frac{u_n{v}_n}{u_n+h}-\mu v_n\right]+v_n\hfill \end{array}$ (28)

其中 $c_k\left(t+1,x\right)$ 、 $u_k\left(t+1,x\right)$ 、 $v_k\left(t+1,x\right)$ 分别表示氧气、浮游植物和浮游动物在x轴方向上的扩散项，即：

$\{\begin{array}{l}{c}_k\left(t+1,x\right)=D_T\frac{c_{\left(t,x-1\right)}+c_{\left(t,x+1\right)}-2c_{\left(t,x\right)}}{\Delta x^2}\hfill \\ u_k\left(t+1,x\right)=D_T\frac{u_{\left(t,x-1\right)}+u_{\left(t,x+1\right)}-2u_{\left(t,x\right)}}{\Delta x^2}\hfill \\ v_k\left(t+1,x\right)=D_T\frac{v_{\left(t,x-1\right)}+v_{\left(t,x+1\right)}-2v_{\left(t,x\right)}}{\Delta x^2}\hfill \end{array}$ (29)

得到最终的差分方程为：

$\{\begin{array}{l}{c}_{\left(t+1,x\right)}=\Delta t\left[c_k\left(t+1,x\right)+A\frac{k_0}{c+k_0}-{\alpha }_1\frac{u_n{c}_n}{c_n+k_2}-{\alpha }_2\frac{v_n{c}_n{u}_n}{c_n+k_3}-c_n\right]+c_n\hfill \\ u_{\left(t+1,x\right)}=\Delta t\left[u_k\left(t+1,x\right)+\left(\frac{{\alpha }_{3}B{c}_n}{c_n+k_1}-\gamma u_n\right)\cdot u_n-\frac{\beta v_n{u}_n}{u_n+h}-\sigma u_n\right]+u_n\hfill \\ v_{\left(t+1,x\right)}=\Delta t\left[v_k\left(t+1,x\right)+\frac{\eta c_n^2}{c_n^2+k_4^2}\cdot \frac{u_n{v}_n}{u_n+h}-\mu v_n\right]+v_n\hfill \\ P\_DMS\left[i\right]=averag{e}_j\left[u\left(i,j\right)\right]\cdot k\cdot \left(1-{\text{e}}^{\frac{\tau }{\Delta T}}\right)\hfill \\ \omega \left[i\right]={\omega }_0-\gamma P\_DMS\left[i\right]\hfill \\ T\left[i+1\right]=T\left[i\right]+\omega \left[i\right]\cdot \Delta t\hfill \\ \Delta T=T\left[i\right]-T\left[0\right]\hfill \end{array}$ (30)

3.2. 有关参数求解

3.2.1. 恒温基础模型求解

2.1中给出了在温度恒定的基础情形下，氧气与浮游生物的变化关系，现对其求解以获取常态化情形时三者的密度大小。

同理对公式(5)采用简单显格式的有限差分法，不再赘述。对于初始条件以及相关参数，我们查阅相关文献 [13] 后采取以下取值，如表2所示。

Table 2. Correlation function expressions

表2. 相关函数表达式

代入对应的差分方程，使用python对其迭代求解，得到稳定状态下浮游生物密度与氧气浓度的关系图，如图2所示。

Figure 2. Graph of phytoplankton density versus oxygen concentration in a stable temperature state

图2. 气温稳定状态下浮游动植物密度与氧气浓度的关系图

观察图2可以看到，浮游植物密度与海水内氧浓度呈同相位正弦变化，而比其浮游动物的正弦波则滞后半个周期，符合捕食关系的种群密度变化，因而基础模型是可靠的。

3.2.2. 最终模型求解

由于基础模型已经给出了温度恒定的情况下相对稳定、合理的解，故最终模型的求解中，c、u、v三者的初始值可以更正为基础模型中三者是随机抽样取值，其余参数保持不变，即如表3所示。

Table 3. Correlation function expressions

表3. 相关函数表达式

代入式(30)，使用python对其迭代求解，得到如图3所示，最终模型下温度、浮游植物密度、浮游动物密度、氧气浓度四者之间的关系图。

Figure 3. Plankton-oxygen-temperature variation model

图3. 浮游生物–氧气–气温变化模型

4. 模型灵敏度分析

对初始的光合作用的强度 $A_0$ 做灵敏度分析，分别得到 $A_0=2.08,2.11,2.14,2.16$ 时四者的关系图，如图4~7所示。

Figure 4. Diagram at $A_0=2.08$

图4. $A_0=2.08$ 时的关系图

Figure 5. Diagram at $A_0=2.11$

图5. $A_0=2.11$ 时的关系图

Figure 6. Diagram at $A_0=2.14$

图6. $A_0=2.14$ 时的关系图

Figure 7. Diagram at $A_0=2.16$

图7. $A_0=2.16$ 时的关系图

图4~7可以看出光合作用强度轻微的扰动反应在曲线上时将造成巨大的变化，尤其当其增长到2.16时，浮游生物将在近百年内灭绝，海洋中的氧气也将消耗殆尽。

5. 结论

由图3可以看出，当浮游植物产氧对气候变化的影响被考虑后，浮游生物与氧气的波动趋于平缓，振幅变得越来越小，并且整体有略微下降的趋势。与此同时，温度上升幅度也趋于平缓，由原本的直线变为斜率逐渐减小的曲线，可以看出二者之间形成相互制约，从而使生态整体更加稳定。在此背景下，未来的500年内该生态系统未见明显的崩溃。

由图4~7可以看出，光合作用强度对生态系统稳定性的影响是巨大的，持续的高强度光照将对生态系统造成毁灭性打击。

6. 建议

出于保护生态的目的，在全球变暖的大背景下，我们应当大力保护生态系统的多样性，避免单一物种充当单一角色的情形，同时积极保护植被，净化二氧化碳、释放氧气，最后还应克制化石燃料的燃烧、使用，减少有毒物质的产生和释放。

需要指出的是，在氧气、生物扩散维度方面，模型还值得进一步改进，即可增加y或(y, z)等维度对模型进行补充，而二氧化碳浓度 [21] 对光合作用以及全球变暖的影响也值得被考虑。