Python手敲耦合地球系统模式

介绍

马普学会气象研究所[Max-Planck-Institut für Meteorologie]拥有先进的地球系统模式，其所属IMPRS-ESM（博士研究生培养单位）每年在全球招收十多名博士。感兴趣的同学可以查看官网。在今年的新生导论课程An introduction to Earth System Modelling中，我们用python手敲了一个地球系统模式，模型虽然简单（因此称为box model），但是对于理解耦合地球系统模式是怎样运行的非常有帮助，比如，你会知道，原来模型真的就是把初值设置为0，然后让它跑。

本文不再赘述理论部分，参考资料见链接https://owncloud.gwdg.de/index.php/s/P53YvEggfM5ezbj，实际上课程上也很少讲这个资料中的部分。1-5纯背景，从第六章开始介绍课程使用的模型。

此处简单介绍一下要敲的模型。

The model we will develop, which for lack of a better word we call the Earth-system Box Model, consists of a coupled set of differential equations describing how changes in the heat and carbon content reservoirs influence one another, and Earth’s surface temperature.

模型基于能量和质量守恒，分别建立 Energy Balance Box Model (EBBM) 和 Mass (carbon) Balance Box Model (MBBM) （Fig 5.1）. 前者仅两层，分别是上层或表层海洋，用字母表示，以及深海，用字母表示。课上bjoren介绍了为什么这个简单的box model不包含大气的合理性。虽然他是atmosphere department的director。海洋的相当深度（ equivalent depth）为2620m，表层占比用参数表示，其余为深海。后者仅考虑碳守恒，包含四层：表示大气，表示陆地，以及表层和深层海洋。我们也可以从下图看出，改模型包含两个外部输入，一个是来自外太空的能量，以及通过燃烧化石能源注入的碳.

IMPRS-BOX MODEL

整个模型由以下六个微分方程构成：

其中，带点的符号代表对时间的微分。表示焓，表示热量（传输），表示碳储（carbon reservoir），表示碳吸收。

类似表示能量（质量）从表层输送到深层。比如，式描述，海表的焓变，等于外部输入的能量的变化（比如返照率变化引起的吸收太阳短波辐射变化）减去从海表传输到深海的能量。式描述，大气中的碳变化，等于化石燃料燃烧引起的碳注入，减去传输到海表的碳，减去传输到陆地的碳。其他公式一样。

上述公式只是一个表示。参考材料第六七章中有比较详细的公式推导，即，怎样从公式得到下列公式，本文不再推导：

以上公式中，斜体表示相对于工业革命前的变化量，比如表示大气中的碳相对于工业革命前的变化量，而花里胡哨的表示大气中碳含量的绝对初值。比如第一个公式表示，海表气温的变化率，由四项决定，第一项表示向外的长波辐射，第二项表示温室效应，第三项表示表层向深层的热量传输，第四项是一个什么辐射项（我也忘了）。虽然此处因为篇幅和时间原因没有详细介绍以上公式描述的过程，但是如果真正想通过这个project增进对模型运作原理的理解，十分建议阅读参考材料，大致掌握以上公式说了什么。

从上面的公式可以看出，该模型是一个耦合的模型，比如表层温度随时间的变化率（对时间的微分，第一行公式）的变化，依赖于海表的温度本身，表层海和深海的温差，以及大气中的碳含量。而大气中的碳含量（最后一行公式）与其他碳储有关，进而与大气温度相关。

Python 手敲耦合模式

实际上，run模型也很简单，我们的任务就是把以上的公式转换成python代码。我们此次用到的一个重要的包是sympy。我的理解，sympy可以将公式表示为我们能理解的，常见的公式的样子（就像是上文展示的公式那样）的同时，可以将公式转化为可以运行的函数。

Import

import numpy as np
import sympy
from matplotlib import pyplot as plt
from scipy import integrate

sympy define symbols and functions

简单理解，就是定义，键盘上输入'tau_l0'，sympy在展示公式的时候展示为。

t,tau_l0,Chi = sympy.symbols("t, τ_l^0, χ")
Pi0,beta_Pi,C_a0, C_l0 = sympy.symbols("Π^0, β_Π, C_a^0, C_l^0")
gam,k_a,kappa_0,EC,delta = sympy.symbols ("γ,k_a,κ_0,ηC,δ")
tau_l = sympy.symbols("τ_l")
R_h = sympy.symbols("R_h")
P = sympy.symbols("P")

c_star,lam,beta,E_H,lam_star= sympy.symbols("c_*,λ,β,η_H,λ_*")

C_a = sympy.Function("C_a")
C_s = sympy.Function("C_s")
C_d = sympy.Function("C_d")
C_l = sympy.Function("C_l")
J   = sympy.Function("J")
T_s = sympy.Function("T_s")
T_d = sympy.Function("T_d")

Define the six equations

Using the symbols we defined above, represent our equations using sympy.

1. Energy

yr2s = 1   # change year to seconds, but here we don't use it.

ode_Ts = 1 / (delta*c_star) * (-lam*T_s(t*yr2s) + beta*sympy.log(C_a(t)/C_a0 + 1) 
                               - E_H*(T_s(t*yr2s)-T_d(t*yr2s)) 
                               - lam_star * (T_s(t*yr2s) - T_d(t*yr2s)))

ode_Ts   # here the sympy output even the same equations as the Latex. but it's not the excutable equation yet.

2. Energy

ode_Td = E_H*(T_s(t*yr2s)-T_d(t*yr2s))/((1-delta) * c_star)

ode_Td

3. function of

function τ

tau_l = tau_l0 * Chi ** (-T_s(t)/10)




    
tau_l

function

R_h = (C_l0+C_l(t))/tau_l

R_h

function

P = Pi0 * (1+beta_Pi*sympy.ln(C_a(t)/C_a0 +1))

function

ode_Cl = P - R_h

ode_Cl

4. function of

ode_Cs = gam * (C_a(t)/k_a - C_s(t)/kappa_0) - EC*(C_s(t)/delta - C_d(t)/(1-delta))

ode_Cs

5. function of

ode_Cd = EC*(C_s(t)/delta - C_d(t)/(1-delta))

ode_Cd

6. function of

To define , we need the emission first. The logistic equation that emulate the historical and future high-end emission sce- nario SSP5-8.5 is as follows (Archer and Brovkin, 2008):

where , . By definition , which can be derived analytically.

A_tot,t_opt = sympy.symbols ("A_tot,t_opt")

A = sympy.Function("A")

A = A_tot*( 1/(1+2.5*sympy.exp((t_opt-t)/50)) - 1/(1+2.5*sympy.exp(t_opt/50))) # rcp8.5 emission

A_ = []
for i in range(600):
    A_.append(A.subs({t:i,t_opt:250,A_tot:5000}))

plt.plot(A_)

Fig. 2 historical and future high-end emission sce- nario SSP5-8.5

J=A.diff(t)

ode_Ca = J - ode_Cl-ode_Cs-ode_Cd

ode_Ca

Equations To functions

上述代码定义了6个方程，接下来，我们要将以上六个方程转变成可以执行的函数，sympy函数sympy.lambdify可以直接调用。

# Create the ODE
ode_sys = [ode_Ts, ode_Td,ode_Cl, ode_Cs, ode_Cd, ode_Ca]




    
ode_sys_np = sympy.lambdify(([T_s(t), T_d(t), C_l(t), C_s(t), C_d(t), C_a(t)],
                            t, 
                             c_star, 
                             lam, 
                             beta, 
                             E_H, 
                             lam_star,
                             tau_l0, 
                             Chi, 
                             Pi0, 
                             beta_Pi, 
                             C_a0, 
                             C_l0, 
                             gam, 
                             k_a, 
                             kappa_0, 
                             EC, 
                             delta, 
                             A_tot,
                             t_opt),
                            ode_sys)

Solve the equations

We use numeric solution here, firstly define some initial conditions. Most of the values for the parameters can be found in the lecture note, for example, the table below:

1. initial condition

ics = { 
    'TC0': [0,0,0,0,0,0], # ode_Ts,ode_Td,ode_Cl, ode_Cs, ode_Cd, ode_Ca
                              
    # var     # values.             # unit (origin)
    c_star:   10.8*1e9,             # JK**-1m**-2 
    lam:      1.77*365*24*3600,     # W m**-2K**-1
    beta:     5.77*365*24*3600,     # W m**-2
    E_H:      0.73*(365*24*3600),   # None
    lam_star: 0,                    # None

    tau_l0:   41,                   # yr
    Chi:      1.8,                  # None     
    Pi0:      60 ,                  # GtC/yr 
    beta_Pi:  0.4,                  # None 
    C_a0:     589,                  # GtC
    C_l0:     2500,                 # GtC
    gam:      0.005,                # GtC/yr/ppm
    k_a:      2.12,                 # GtC/ppm
    kappa_0:  1.98,                 # GtC/ppm
    EC:       60*1e-12*365*24*3600, # s ** -1 
    delta:    0.015 ,               # None
    A_tot:    5000,                 # GtC
    t_opt:    250 ,                 # yr   
    
}
t_ = np.arange(0,5000,1)   # year

2. Solve the ode

xy_t = integrate.odeint(ode_sys_np, ics['TC0'],
                        t_, 
                        args=(
                        ics[c_star],
                        ics[lam],
                        ics[beta],
                        ics[E_H],
                        ics[lam_star],
                        ics[tau_l0],
                        ics[Chi],
                        ics[Pi0],
                        ics[beta_Pi],
                        ics[C_a0],
                        ics[C_l0],                        
                        ics[gam],
                        ics[k_a],
                        ics[kappa_0],  
                        ics[EC],
                        ics[delta],
                        ics[A_tot],
                        ics[t_opt]
                        ))




    
xy_t.shape

(5000, 6)

Results

fig,axes = plt.subplots(3,2,figsize = (12,10),dpi = 100)

colors = ['--k','--b','g','k','b','r']
labels = [r'$T_s$',r'$T_d$',r'$C_l$',r'$C_s$',r'$C_d$',r'$C_a$']

for i,ax in enumerate(axes.flat):
    ax.plot(xy_t[:,i],colors[i],label = labels[i])
    ax.legend(loc = "upper left")
    ax.set_xlabel("year")
    
for ax in axes.flat[:2]:
    ax.set_ylabel("T anomaly (K)")
for ax in axes.flat[2:]:
    ax.set_ylabel("C anomaly (Gt)")

至此，我们就完成了一个简单的耦合地球系统模式的建模。我们可以简单的解读这个模型，图中虚线代表温度的变化，实线代表碳的变化，值得一提的是，此处纵坐标代表的是相对于工业革命前的变化量。结合之前的(我们标为Fig. 2)中的排放模型，排放增加，首先海表温度快速升高（左1），深海温度变化较慢（右1）。海表最高升温4.7度左右，深海小于4度。最终大概在4000年之后再次达到稳定状态。从碳的角度来看，化石燃料燃烧，陆地和大气中的碳都是先升后降（左2和右3），大部分的碳最中被海洋吸收（右2和左3）。

虽然这只是一个简单的box model，有大量的参数化过程，所以也算不上是物理模型。但是对于帮助我理解耦合模型是如何运作的帮助很大。最后全班被分成每两人的小组，每个小组都展示了他们的project，要么改造了模型，比如我们小组将海表分为了冷热两部分模拟赤道和极地，有小组添加了海冰等，要么利用该模型做了一些应用。非常有创造力。