The two-layer grey radiation model (aka the leaky greenhouse)

ATM 623: Climate Modeling

As of IPython 2.0, the user interface has changed significantly
Enter edit mode by pressing `Enter`
Don't try to type into a cell in command mode
  • Surface temperature is $T_s$
  • Atm. temperatures are $T_1, T_2$ where $T_1$ is closest to the surface.
  • absorptivity of atm layers is $\epsilon_1, \epsilon_2$
  • Surface emission is $\sigma T_s^4$
  • Atm emission is $\epsilon_1 \sigma T_1^4, \epsilon_2 \sigma T_2^4$ (up and down)
  • Abs = emissivity for atm layers
  • Transmissivity for atm layers is $\tau_1, \tau_2$ where $\tau_i = (1-\epsilon_i)$

Emission

$$ E_s = \sigma T_s^4 $$$$ E_1 = \epsilon_1 \sigma T_1^4 $$$$ E_2 = \epsilon_2 \sigma T_2^4 $$

Incident radiation

$$ F_s = \tau_1 E_2 + E_1 $$$$ F_1 = E_s + E_2 $$$$ F_2 = \tau_1 E_s + E_1 $$

Net radiation

(absorptivity) * incident - emission

$$ R_s = F_s - E_s $$$$ R_1 = \epsilon_1 F_1 - 2 E_1 $$$$ R_2 = \epsilon_2 F_2 - 2 E_2 $$

OLR

$$ OLR = \tau_2 F_2 + E_2 $$$$ = \tau_1 \tau_2 E_s + \tau_2 E_1 + E_2 $$$$ = \tau_1 \tau_2 \sigma T_s^4 + \tau_2 \epsilon_1 \sigma T_1^4 + \epsilon_2 \sigma T_2^4 $$

Net radiation in terms of emissions

$$ R_s = \tau_1 E_2 + E_1 - E_s $$$$ R_1 = \epsilon_1 (E_s + E_2) - 2 E_1 $$$$ R_2 = \epsilon_2 (\tau_1 E_s + E_1) - 2 E_2 $$

Net radiation in terms of temperatures

$$ R_s = \tau_1 \epsilon_2 \sigma T_2^4 + \epsilon_1 \sigma T_1^4 - \sigma T_s^4 $$$$ R_1 = \epsilon_1 (\sigma T_s^4 + \epsilon_2 \sigma T_2^4) - 2 \epsilon_1 \sigma T_1^4 $$$$ R_2 = \epsilon_2 (\tau_1 \sigma T_s^4 + \epsilon_1 \sigma T_1^4) - 2 \epsilon_2 \sigma T_2^4 $$

Net radiation in terms of temperatures and absorptivities

$$ R_s = (1-\epsilon_1) \epsilon_2 \sigma T_2^4 + \epsilon_1 \sigma T_1^4 - \sigma T_s^4 $$$$ R_1 = \epsilon_1 (\sigma T_s^4 + \epsilon_2 \sigma T_2^4) - 2 \epsilon_1 \sigma T_1^4 $$$$ R_2 = \epsilon_2 ((1-\epsilon_1) \sigma T_s^4 + \epsilon_1 \sigma T_1^4) - 2 \epsilon_2 \sigma T_2^4 $$

Solve for radiative equilibrium

Need to add the solar energy source. We assume atm is transparent, solar is all absorbed at the surface.

$$ R_1 = R_2 = 0$$$$ R_s = - (1-\alpha) Q $$

Introduce useful notation shorthand:

$$ (1-\alpha) Q = \sigma T_e^4 $$

This gives a 3x3 system which is linear in $T^4$ (divide through by $\sigma$)

$$ - T_s^4 + \epsilon_1 T_1^4 + (1-\epsilon_1) \epsilon_2 T_2^4 + T_e^4 = 0 $$$$ \epsilon_1 T_s^4 - 2 \epsilon_1 T_1^4 + \epsilon_1 \epsilon_2 T_2^4 = 0$$$$ \epsilon_2 (1-\epsilon_1) T_s^4 + \epsilon_1 \epsilon_2 T_1^4 - 2 \epsilon_2 T_2^4 = 0$$

Here we use the sympy module to solve the algebraic system symbolically.


In [1]:
import sympy
sympy.init_printing()
T_s, T_1, T_2, T_e, e_1, e_2 = sympy.symbols('T_s, T_1, T_2, T_e, e_1, e_2', positive=True )
system = [-T_s**4 + e_1*T_1**4 + e_2*(1-e_1)*T_2**4 + T_e**4,
              e_1*T_s**4 - 2*e_1*T_1**4 + e_1*e_2*T_2**4,
              e_2*(1-e_1)*T_s**4 + e_1*e_2*T_1**4 - 2*e_2*T_2**4]
out1 = sympy.solve( system, [T_s**4, T_1**4, T_2**4])
out1


Out[1]:
$$\left \{ T_{1}^{4} : \frac{T_{e}^{4} \left(- e_{1} e_{2} + e_{2} + 2\right)}{e_{1} e_{2} - 2 e_{1} - 2 e_{2} + 4}, \quad T_{2}^{4} : - \frac{T_{e}^{4}}{e_{2} - 2}, \quad T_{s}^{4} : \frac{T_{e}^{4} \left(- e_{1} e_{2} + 4\right)}{e_{1} e_{2} - 2 e_{1} - 2 e_{2} + 4}\right \}$$

In [2]:
quarter = sympy.Rational(1,4)
out2 = {}
for var4, formula in out1.iteritems():
    var = (var4)**quarter
    out2[var] = sympy.simplify(formula**quarter)
out2


Out[2]:
$$\left \{ T_{1} : T_{e} \sqrt[4]{\frac{- e_{1} e_{2} + e_{2} + 2}{e_{1} e_{2} - 2 e_{1} - 2 e_{2} + 4}}, \quad T_{2} : T_{e} \sqrt[4]{- \frac{1}{e_{2} - 2}}, \quad T_{s} : T_{e} \sqrt[4]{\frac{- e_{1} e_{2} + 4}{e_{1} e_{2} - 2 e_{1} - 2 e_{2} + 4}}\right \}$$

In [3]:
#  The special case of equal absorptivities
e = sympy.symbols('e')
out3 = {}
for var4, formula in out1.iteritems():
    var = (var4)**quarter
    simple_formula = sympy.cancel(formula.subs([(e_2, e),(e_1, e)]))
    out3[var] = sympy.simplify( simple_formula**quarter )
out3


Out[3]:
$$\left \{ T_{1} : T_{e} \sqrt[4]{\frac{- e - 1}{e - 2}}, \quad T_{2} : T_{e} \sqrt[4]{- \frac{1}{e - 2}}, \quad T_{s} : T_{e} \sqrt[4]{\frac{- e - 2}{e - 2}}\right \}$$

The solution is

\begin{align} T_s^4 &= T_e^4 \frac{4 - \epsilon_1 \epsilon_2}{4 + \epsilon_1 \epsilon_2 - 2 \epsilon_1 - 2 \epsilon_2} \\ T_1^4 &= T_e^4 \frac{2 -\epsilon_1 \epsilon_2 + \epsilon_2}{4 + \epsilon_1 \epsilon_2 - 2 \epsilon_1 - 2 \epsilon_2} \\ T_2^4 &= T_e^4 \frac{ 1}{2 - \epsilon_2} \end{align}

In the special case $\epsilon_1 = \epsilon_2$ this reduces to

\begin{align} T_s^4 &= T_e^4 \frac{2+\epsilon}{2-\epsilon} \\ T_1^4 &= T_e^4 \frac{1+\epsilon}{2-\epsilon} \\ T_2^4 &= T_e^4 \frac{ 1}{2 - \epsilon} \end{align}

In [4]:
out2[T_s].subs([(T_e, 255), (e_1, 0.4), (e_2, 0.4)])


Out[4]:
$$282.203889523582$$

In [5]:
for var, formula in out2.iteritems():
    print formula.subs([(T_e, 255), (e_1, 0.4), (e_2, 0.4)])


246.627893584128
282.203889523582
226.730624779963

In [6]:
# Coding up the analytical solutions for radiative equilibrium
#  These use the analytical results returned by sympy and wrap them in callable functions

def Ts(Te, e1, e2):
    #return Te*((4-e1*e2)/(4+e1*e2-2*(e1+e2)))**0.25
    return out2[T_s].subs([(T_e, Te), (e_1, e1), (e_2, e2)])
def T1(Te, e1, e2):
    #return Te*((2+e2-e1*e2)/(4+e1*e2-2*(e1+e2)))**0.25
    return out2[T_1].subs([(T_e, Te), (e_1, e1), (e_2, e2)])
def T2(Te, e1, e2):
    #return Te*(1/(2-e2))**0.25
    return out2[T_2].subs([(T_e, Te), (e_1, e1), (e_2, e2)])

In [7]:
%matplotlib inline
import numpy as np
from climlab import constants as const
from climlab.model import column

In [8]:
mycolumn = column.GreyRadiationModel( num_lev=2 )

In [9]:
print mycolumn


climlab Process of type <class 'climlab.model.column.GreyRadiationModel'>. 
State variables and domain shapes: 
  Tatm: (2,) 
  Ts: (1,) 
The subprocess tree: 
top: <class 'climlab.model.column.GreyRadiationModel'>
   LW: <class 'climlab.radiation.greygas.GreyGas'>
   SW: <class 'climlab.radiation.greygas.GreyGasSW'>
   insolation: <class 'climlab.radiation.insolation.FixedInsolation'>


In [10]:
mycolumn.integrate_years(10.)


Integrating for 3652 steps, 3652.422 days, or 10.0 years.
Total elapsed time is 9.99884460229 years.

In [11]:
print mycolumn.Ts
print mycolumn.Tatm


[ 287.84605967]
[ 229.43643402  252.95019305]

In [12]:
(e1, e2)= mycolumn.subprocess['LW'].absorptivity
print e1, e2


0.477374247427 0.477374247427

In [13]:
ASR = (1-mycolumn.param['albedo_sfc'])*mycolumn.param['Q']
Te = (ASR/const.sigma)**0.25
print Te


254.865280431

Check numerical versus analytical results

Use a tolerance value to test if the results are the same.


In [14]:
tol = 0.01

def test_2level(col):
    (e1, e2)= col.subprocess['LW'].absorptivity
    ASR = (1-col.param['albedo_sfc'])*col.param['Q']
    Te = (ASR/const.sigma)**0.25
    print 'Surface:'
    num = col.Ts
    anal = Ts(Te,e1,e2)
    print '  Numerical: %.2f   Analytical: %.2f    Same:' %(num, anal) , abs(num - anal)<tol
    print 'Level 1'
    num = col.Tatm[0]
    anal = T1(Te,e1,e2)
    print '  Numerical: %.2f   Analytical: %.2f    Same:' %(num, anal) , abs(num - anal)<tol
    print 'Level 2'
    num = col.Tatm[1]
    anal = T2(Te,e1,e2)
    print '  Numerical: %.2f   Analytical: %.2f    Same:' %(num, anal) , abs(num - anal)<tol

In [15]:
test_2level(mycolumn)


Surface:
  Numerical: 287.85   Analytical: 287.85    Same: [ True]
Level 1
  Numerical: 229.44   Analytical: 252.95    Same: False
Level 2
  Numerical: 252.95   Analytical: 229.44    Same: False

In [16]:
e1 = 0.3
e2 = 0.4
mycolumn.subprocess['LW'].absorptivity = np.array([e1,e2])
mycolumn.integrate_years(10.)


Integrating for 3652 steps, 3652.422 days, or 10.0 years.
Total elapsed time is 19.9976892046 years.

In [17]:
test_2level(mycolumn)


Surface:
  Numerical: 278.53   Analytical: 278.53    Same: [ True]
Level 1
  Numerical: 223.20   Analytical: 243.87    Same: False
Level 2
  Numerical: 241.15   Analytical: 226.61    Same: False

In [18]:
e1 = 0.6
e2 = 0.6
mycolumn.subprocess['LW'].absorptivity = np.array([e1,e2])
mycolumn.integrate_years(10.)
test_2level(mycolumn)


Integrating for 3652 steps, 3652.422 days, or 10.0 years.
Total elapsed time is 29.9965338069 years.
Surface:
  Numerical: 297.52   Analytical: 297.52    Same: [ True]
Level 1
  Numerical: 234.30   Analytical: 263.52    Same: False
Level 2
  Numerical: 263.52   Analytical: 234.30    Same: False

In [19]:
col1 = column.GreyRadiationModel(num_lev=2, abs_coeff=1.9E-4)
(e1, e2) = col1.subprocess['LW'].absorptivity
print e1, e2
col1.integrate_years(10.)
test_2level(col1)


0.652920962199 0.652920962199
Integrating for 3652 steps, 3652.422 days, or 10.0 years.
Total elapsed time is 9.99884460229 years.
Surface:
  Numerical: 301.92   Analytical: 301.92    Same: [ True]
Level 1
  Numerical: 236.57   Analytical: 268.24    Same: False
Level 2
  Numerical: 268.24   Analytical: 236.57    Same: False

In [20]:
col1 = column.GreyRadiationModel(num_lev=2, abs_coeff=1.9E-4)
lw = col1.subprocess['LW']
(e1, e2) = lw.absorptivity
e1 *= 1.2
lw.absorptivity = np.array([e1, e2])
col1.integrate_years(10.)
test_2level(col1)


Integrating for 3652 steps, 3652.422 days, or 10.0 years.
Total elapsed time is 9.99884460229 years.
Surface:
  Numerical: 307.85   Analytical: 307.85    Same: [ True]
Level 1
  Numerical: 242.68   Analytical: 272.49    Same: False
Level 2
  Numerical: 276.56   Analytical: 236.57    Same: False

Conclusion: The Two-level model works

Three-layer model

Extend the analysis to three layers. Start numbering the layers from 0 to be consistent with array indexing.


In [21]:
T_s, T_0, T_1, T_2, T_e, epsilon_0, epsilon_1, epsilon_2, sigma = \
    sympy.symbols('T_s, T_0, T_1, T_2, T_e, epsilon_0, epsilon_1, epsilon_2, sigma', positive=True )

Define the transmissivities $\tau_i$ for layers $i=0, 1, \dots, N-1$


In [22]:
tau_0 = (1-epsilon_0)
tau_1 = (1-epsilon_1)
tau_2 = (1-epsilon_2)
tau_0, tau_1, tau_2


Out[22]:
$$\left ( - \epsilon_{0} + 1, \quad - \epsilon_{1} + 1, \quad - \epsilon_{2} + 1\right )$$

Note that if the atmosphere has $N$ layers then $\epsilon_N = 0$

Define the emissions for each layer:


In [23]:
E_s = sigma*T_s**4
E_0 = epsilon_0*sigma*T_0**4
E_1 = epsilon_1*sigma*T_1**4
E_2 = epsilon_2*sigma*T_2**4
E_s, E_0, E_1, E_2


Out[23]:
$$\left ( T_{s}^{4} \sigma, \quad T_{0}^{4} \epsilon_{0} \sigma, \quad T_{1}^{4} \epsilon_{1} \sigma, \quad T_{2}^{4} \epsilon_{2} \sigma\right )$$

Define the longwave fluxes incident on each layer, $F_{i}$

Note that if the atmosphere has $N$ layers then $F_{N}$ is the OLR (emission to space)


In [24]:
F_s = E_0 + tau_0*E_1 + tau_0*tau_1*E_2
F_0 = E_s + E_1 + tau_1*E_2
F_1 = tau_0*E_s + E_0 + E_2
F_2 = tau_1*tau_0*E_s + tau_1*E_0 + E_1
F_3 = tau_2*tau_1*tau_0*E_s + tau_2*tau_1*E_0 + tau_2*E_1 + E_2

F_s, F_0, F_1, F_2, F_3


Out[24]:
$$\left ( T_{0}^{4} \epsilon_{0} \sigma + T_{1}^{4} \epsilon_{1} \sigma \left(- \epsilon_{0} + 1\right) + T_{2}^{4} \epsilon_{2} \sigma \left(- \epsilon_{0} + 1\right) \left(- \epsilon_{1} + 1\right), \quad T_{1}^{4} \epsilon_{1} \sigma + T_{2}^{4} \epsilon_{2} \sigma \left(- \epsilon_{1} + 1\right) + T_{s}^{4} \sigma, \quad T_{0}^{4} \epsilon_{0} \sigma + T_{2}^{4} \epsilon_{2} \sigma + T_{s}^{4} \sigma \left(- \epsilon_{0} + 1\right), \quad T_{0}^{4} \epsilon_{0} \sigma \left(- \epsilon_{1} + 1\right) + T_{1}^{4} \epsilon_{1} \sigma + T_{s}^{4} \sigma \left(- \epsilon_{0} + 1\right) \left(- \epsilon_{1} + 1\right), \quad T_{0}^{4} \epsilon_{0} \sigma \left(- \epsilon_{1} + 1\right) \left(- \epsilon_{2} + 1\right) + T_{1}^{4} \epsilon_{1} \sigma \left(- \epsilon_{2} + 1\right) + T_{2}^{4} \epsilon_{2} \sigma + T_{s}^{4} \sigma \left(- \epsilon_{0} + 1\right) \left(- \epsilon_{1} + 1\right) \left(- \epsilon_{2} + 1\right)\right )$$

Now define the net absorbed longwave radiation (flux divergence) in each layer.


In [25]:
R_s = F_s - E_s
R_0 = epsilon_0*F_0 - 2*E_0
R_1 = epsilon_1*F_1 - 2*E_1
R_2 = epsilon_2*F_2 - 2*E_2

R_s, R_0, R_1, R_2


Out[25]:
$$\left ( T_{0}^{4} \epsilon_{0} \sigma + T_{1}^{4} \epsilon_{1} \sigma \left(- \epsilon_{0} + 1\right) + T_{2}^{4} \epsilon_{2} \sigma \left(- \epsilon_{0} + 1\right) \left(- \epsilon_{1} + 1\right) - T_{s}^{4} \sigma, \quad - 2 T_{0}^{4} \epsilon_{0} \sigma + \epsilon_{0} \left(T_{1}^{4} \epsilon_{1} \sigma + T_{2}^{4} \epsilon_{2} \sigma \left(- \epsilon_{1} + 1\right) + T_{s}^{4} \sigma\right), \quad - 2 T_{1}^{4} \epsilon_{1} \sigma + \epsilon_{1} \left(T_{0}^{4} \epsilon_{0} \sigma + T_{2}^{4} \epsilon_{2} \sigma + T_{s}^{4} \sigma \left(- \epsilon_{0} + 1\right)\right), \quad - 2 T_{2}^{4} \epsilon_{2} \sigma + \epsilon_{2} \left(T_{0}^{4} \epsilon_{0} \sigma \left(- \epsilon_{1} + 1\right) + T_{1}^{4} \epsilon_{1} \sigma + T_{s}^{4} \sigma \left(- \epsilon_{0} + 1\right) \left(- \epsilon_{1} + 1\right)\right)\right )$$

Solve for radiative equilibrium

Use sympy.solve to automatically solve the algebraic system.

We will solve for the radiative equilibrium temperatures in two steps:

  • First solve for $T_i^4$, which is a purely linear problem.
  • Then take the fourth roots to solve for the temperatures.

In [26]:
out1 = sympy.solve([R_s + sigma*T_e**4, R_0, R_1, R_2],
            [T_s**4, T_0**4, T_1**4, T_2**4])
out1


Out[26]:
$$\left \{ T_{0}^{4} : \frac{T_{e}^{4} \left(- 2 \epsilon_{0} \epsilon_{1} \epsilon_{2} + 2 \epsilon_{0} \epsilon_{1} + 2 \epsilon_{0} \epsilon_{2} + 3 \epsilon_{1} \epsilon_{2} - 2 \epsilon_{1} - 2 \epsilon_{2} - 4\right)}{\epsilon_{0} \epsilon_{1} \epsilon_{2} - 2 \epsilon_{0} \epsilon_{1} - 2 \epsilon_{0} \epsilon_{2} + 4 \epsilon_{0} - 2 \epsilon_{1} \epsilon_{2} + 4 \epsilon_{1} + 4 \epsilon_{2} - 8}, \quad T_{1}^{4} : \frac{T_{e}^{4} \left(- \epsilon_{1} \epsilon_{2} + \epsilon_{2} + 2\right)}{\epsilon_{1} \epsilon_{2} - 2 \epsilon_{1} - 2 \epsilon_{2} + 4}, \quad T_{2}^{4} : - \frac{T_{e}^{4}}{\epsilon_{2} - 2}, \quad T_{s}^{4} : \frac{2 T_{e}^{4} \left(- \epsilon_{0} \epsilon_{1} \epsilon_{2} + \epsilon_{0} \epsilon_{1} + \epsilon_{0} \epsilon_{2} + \epsilon_{1} \epsilon_{2} - 4\right)}{\epsilon_{0} \epsilon_{1} \epsilon_{2} - 2 \epsilon_{0} \epsilon_{1} - 2 \epsilon_{0} \epsilon_{2} + 4 \epsilon_{0} - 2 \epsilon_{1} \epsilon_{2} + 4 \epsilon_{1} + 4 \epsilon_{2} - 8}\right \}$$

In [27]:
quarter = sympy.Rational(1,4)
out2 = {}
for var4, formula in out1.iteritems():
    var = (var4)**quarter
    out2[var] = sympy.simplify(formula**quarter)
out2


Out[27]:
$$\left \{ T_{0} : T_{e} \sqrt[4]{\frac{- 2 \epsilon_{0} \epsilon_{1} \epsilon_{2} + 2 \epsilon_{0} \epsilon_{1} + 2 \epsilon_{0} \epsilon_{2} + 3 \epsilon_{1} \epsilon_{2} - 2 \epsilon_{1} - 2 \epsilon_{2} - 4}{\epsilon_{0} \epsilon_{1} \epsilon_{2} - 2 \epsilon_{0} \epsilon_{1} - 2 \epsilon_{0} \epsilon_{2} + 4 \epsilon_{0} - 2 \epsilon_{1} \epsilon_{2} + 4 \epsilon_{1} + 4 \epsilon_{2} - 8}}, \quad T_{1} : T_{e} \sqrt[4]{\frac{- \epsilon_{1} \epsilon_{2} + \epsilon_{2} + 2}{\epsilon_{1} \epsilon_{2} - 2 \epsilon_{1} - 2 \epsilon_{2} + 4}}, \quad T_{2} : T_{e} \sqrt[4]{- \frac{1}{\epsilon_{2} - 2}}, \quad T_{s} : \sqrt[4]{2} T_{e} \sqrt[4]{\frac{- \epsilon_{0} \epsilon_{1} \epsilon_{2} + \epsilon_{0} \epsilon_{1} + \epsilon_{0} \epsilon_{2} + \epsilon_{1} \epsilon_{2} - 4}{\epsilon_{0} \epsilon_{1} \epsilon_{2} - 2 \epsilon_{0} \epsilon_{1} - 2 \epsilon_{0} \epsilon_{2} + 4 \epsilon_{0} - 2 \epsilon_{1} \epsilon_{2} + 4 \epsilon_{1} + 4 \epsilon_{2} - 8}}\right \}$$

Now wrap these analytical radiative equilibrium solutions in callable functions:


In [28]:
def Ts(Te, e0, e1, e2):
    return out2[T_s].subs([(T_e, Te), (epsilon_0, e0), (epsilon_1, e1), (epsilon_2, e2)])
def T0(Te, e0, e1, e2):
    return out2[T_0].subs([(T_e, Te), (epsilon_0, e0), (epsilon_1, e1), (epsilon_2, e2)])
def T1(Te, e0, e1, e2):
    return out2[T_1].subs([(T_e, Te), (epsilon_0, e0), (epsilon_1, e1), (epsilon_2, e2)])
def T2(Te, e0, e1, e2):
    return out2[T_2].subs([(T_e, Te), (epsilon_0, e0), (epsilon_1, e1), (epsilon_2, e2)])

Compare numerical and analytical solutions for radiative equilibrium

Define a function that takes a climlab.column.GreyRadiationModel object (which should be first integrated out to equilibrium), and compares the numerical solution to our analytical solution.


In [29]:
tol = 0.01

def test_3level(col):
    (e0, e1, e2)= col.subprocess['LW'].absorptivity
    ASR = (1-col.param['albedo_sfc'])*col.param['Q']
    Te = (ASR/const.sigma)**0.25
    print 'Surface:'
    num = col.Ts
    anal = Ts(Te,e0,e1,e2)
    print '  Numerical: %.2f   Analytical: %.2f    Same:' %(num, anal) , abs(num - anal)<tol
    print 'Level 0'
    num = col.Tatm[0]
    anal = T0(Te,e0,e1,e2)
    print '  Numerical: %.2f   Analytical: %.2f    Same:' %(num, anal) , abs(num - anal)<tol
    print 'Level 1'
    num = col.Tatm[1]
    anal = T1(Te,e0,e1,e2)
    print '  Numerical: %.2f   Analytical: %.2f    Same:' %(num, anal) , abs(num - anal)<tol
    print 'Level 2'
    num = col.Tatm[2]
    anal = T2(Te,e0,e1,e2)
    print '  Numerical: %.2f   Analytical: %.2f    Same:' %(num, anal) , abs(num - anal)<tol

In [30]:
col = column.GreyRadiationModel( num_lev=3 )
col.integrate_years(10.)
test_3level(col)


Integrating for 3652 steps, 3652.422 days, or 10.0 years.
Total elapsed time is 9.99884460229 years.
Surface:
  Numerical: 287.85   Analytical: 287.85    Same: [ True]
Level 0
  Numerical: 224.73   Analytical: 256.29    Same: False
Level 1
  Numerical: 242.05   Analytical: 242.05    Same: True
Level 2
  Numerical: 256.29   Analytical: 224.73    Same: False

In [31]:
e0 = 0.3
e1 = 0.6
e2 = 0.2
col.subprocess['LW'].absorptivity = np.array([e0,e1,e2])
col.integrate_years(10.)
test_3level(col)


Integrating for 3652 steps, 3652.422 days, or 10.0 years.
Total elapsed time is 19.9976892046 years.
Surface:
  Numerical: 291.71   Analytical: 291.71    Same: [ True]
Level 0
  Numerical: 223.20   Analytical: 262.65    Same: False
Level 1
  Numerical: 247.60   Analytical: 242.93    Same: False
Level 2
  Numerical: 264.53   Analytical: 220.04    Same: False

Conclusion: The three-level model works


In [32]:
# The 3-layer solution reduces to two layer solution if we set e_2 = 0
out3 = {}
for var, formula in out2.iteritems():
    if var is not T_2:
        out3[var] = sympy.simplify(formula.subs(epsilon_2,0))
out3


Out[32]:
$$\left \{ T_{0} : T_{e} \sqrt[4]{\frac{- \epsilon_{0} \epsilon_{1} + \epsilon_{1} + 2}{\epsilon_{0} \epsilon_{1} - 2 \epsilon_{0} - 2 \epsilon_{1} + 4}}, \quad T_{1} : T_{e} \sqrt[4]{- \frac{1}{\epsilon_{1} - 2}}, \quad T_{s} : T_{e} \sqrt[4]{\frac{- \epsilon_{0} \epsilon_{1} + 4}{\epsilon_{0} \epsilon_{1} - 2 \epsilon_{0} - 2 \epsilon_{1} + 4}}\right \}$$

Generalizing to N layers

\begin{align} E_i &= \epsilon_i \sigma T_i^4 \\ F_i &= ... + \tau_{i+2} \tau_{i+1} E_{i+3} + \tau_{i+1} E_{i+2} + E_{i+1} + E_{i-1} + \tau_{i-1} E_{i-2} + \tau_{i-1} \tau_{i-2} E_{i-3} + ... \\ F_i &= \sum_{n=2}^{N-i} \bigg( \prod_{j=1}^{n-1} \tau_{i+j} \bigg) E_{i+n} + E_{i+1} + E_{i-1} + \sum_{n=2}^{i-1} \bigg( \prod_{j=1}^{n-1} \tau_{i-j} \bigg) E_{i-n} \\ &= \sum_{n=1}^{N-i} \bigg( \prod_{j=0}^{n-1} \tau_{i+j} \bigg) E_{i+n} / \tau_{i} + \sum_{n=1}^{i-1} \bigg( \prod_{j=0}^{n-1} \tau_{i-j} \bigg) E_{i-n} / \tau_i \\ &= \frac{1}{\tau_i} \left\{ \sum_{n=1}^{N-i} \bigg( \prod_{j=0}^{n-1} \tau_{i+j} \bigg) E_{i+n} + \sum_{n=1}^{i-1} \bigg( \prod_{j=0}^{n-1} \tau_{i-j} \bigg) E_{i-n} \right\} \end{align}

Now substitute $n \rightarrow -n$ in the second summation.

\begin{align} F_i &= \frac{1}{\tau_i} \left\{ \sum_{n=1}^{N-i} \bigg( \prod_{j=0}^{n-1} \tau_{i+j} \bigg) E_{i+n} + \sum_{-n=1}^{i-1} \bigg( \prod_{j=0}^{-n-1} \tau_{i-j} \bigg) E_{i+n} \right\} \\ &= \frac{1}{\tau_i} \left\{ \sum_{n=1}^{N-i} \bigg( \prod_{j=0}^{n-1} \tau_{i+j} \bigg) E_{i+n} + \sum_{n=1-i}^{-1} \bigg( \prod_{j=0}^{-n-1} \tau_{i-j} \bigg) E_{i+n} \right\} \end{align}

And substitude $j \rightarrow -j$ in the second product.

\begin{align} F_i &= \frac{1}{\tau_i} \left\{ \sum_{n=1}^{N-i} \bigg( \prod_{j=0}^{n-1} \tau_{i+j} \bigg) E_{i+n} + \sum_{n=1-i}^{-1} \bigg( \prod_{-j=0}^{-n-1} \tau_{i+j} \bigg) E_{i+n} \right\} \\ &= \frac{1}{\tau_i} \left\{ \sum_{n=1}^{N-i} \bigg( \prod_{j=0}^{n-1} \tau_{i+j} \bigg) E_{i+n} + \sum_{n=1-i}^{-1} \bigg( \prod_{j=0}^{1+n} \tau_{i+j} \bigg) E_{i+n} \right\} \\ &= \frac{1}{\tau_i} \left\{ \sum_{n=1}^{N-i} \bigg( \prod_{j=0}^{|n|-1} \tau_{i+j} \bigg) E_{i+n} + \sum_{n=1-i}^{-1} \bigg( \prod_{j=0}^{1- |n|} \tau_{i+j} \bigg) E_{i+n} \right\} \\ &= \frac{1}{\tau_i} \left\{ \sum_{n=1}^{N-i} \bigg( \prod_{j=0}^{|n|-1} \tau_{i+j} \bigg) E_{i+n} + \sum_{n=1-i}^{-1} \bigg( \prod_{j=0}^{sign(n) (|n|-1)} \tau_{i+j} \bigg) E_{i+n} \right\} \end{align}

Now combine both sums together

\begin{align} F_i &= \frac{1}{\tau_i} \left\{ \sum_{n=1}^{N-i} \bigg( \prod_{j=0}^{|n|-1} \tau_{i+j} \bigg) E_{i+n} + E_i + \sum_{n=1-i}^{-1} \bigg( \prod_{j=0}^{sign(n) (|n|-1)} \tau_{i+j} \bigg) E_{i+n} \right\} - \frac{E_i}{\tau_i} \\ &= \frac{1}{\tau_i} \left\{ \sum_{n=1-i}^{N-i} \bigg( \prod_{j=0}^{sign(n) \big(|n|-1\big)} \tau_{i+j} \bigg) E_{i+n} - E_i \right\} \end{align}

with the convention that $\prod_{j=0}^{-1} = 1$.

Alternatively if we set $\prod_{j=0}^{-1} = 0$ then

$$ F_i = \frac{1}{\tau_i} \left\{ \sum_{n=1-i}^{N-i} \bigg( \prod_{j=0}^{sign(n) \big(|n|-1\big)} \tau_{i+j} \bigg) E_{i+n} \right\} $$

Let's use this notation since it simplifies our expressions.

Net radiation

(absorptivity) * incident - emission

\begin{align} R_i &= \epsilon_i F_i - 2 E_i \\ &= \frac{\epsilon_i}{\tau_i} \left\{ \sum_{n=1-i}^{N-i} \bigg( \prod_{j=0}^{sign(n) \big(|n|-1\big)} \tau_{i+j} \bigg) E_{i+n} \right\} - 2 E_i \end{align}

Alternative...

We will define the transmissivity between layer i and layer i+n

$$ T_{in} = \left\{ \begin{array}{cc} \prod_{j=1}^{n-1} \tau_{i+j} & n > 1 \\ 1 & n = 1 \\ 0 & n = 0 \\ 1 & n = -1 \\ \prod_{j=1}^{-n-1} \tau_{i-j} & n < -1 \end{array} \right\} $$

Then the incident flux follows directly

$$ F_i = \sum_{n=1-i}^{N-i} T_{in} E_{i+n} $$

and the net radiation is

\begin{align} R_i &= \epsilon_i F_i - 2 E_i \\ &= \epsilon_i \sum_{n=1-i}^{N-i} T_{in} E_{i+n} - 2 E_i \end{align}

Now make the substitution $i+n \rightarrow m$

\begin{align} F_i &= \sum_{m=1}^{N} T_{im} E_{m} \\ T_{im} &= \left\{ \begin{array}{cc} \prod_{j=1}^{m-i-1} \tau_{i+j} & m > 1+i \\ 1 & m = i+1 \\ 0 & m = i \\ 1 & m = i-1 \\ \prod_{j=1}^{i-m-1} \tau_{i-j} & m < i-1 \end{array} \right\} \end{align}

Or using the Einstein summation notation, since $m$ is a repeated index, we can just write

$$ F_i = T_{im} E_{m} $$

and the net radiation is

$$ R_i = \epsilon_i F_i - 2 E_i $$

Reformulating in terms of flux between layers

Let the upwelling flux be a vector ${\bf{U}} = [U_0, U_1, ..., U_{N-1}, U_N]$.

If there are $N$ levels then $\bf{U}$ has $N+1$ elements. We will number the layers starting from 0 following numpy index conventions.

  • $U_0$ is the upwelling flux from surface to layer 0.
  • $U_1$ is the upwelling flux layer 0 to layer 1, etc.
  • $U_N$ is the upwelling flux from layer N-1 (the top level) to space.

Same for the downwelling flux ${\bf{D}} = [D_0, D_1, ..., D_N]$. So $D_N$ is the flux down from space and $D_0$ is the backradiation to the surface.

The absorptivity vector is ${\bf{\epsilon}} = [\epsilon_0, \epsilon_1, ..., \epsilon_{N-1}]$ ($N$ elements)


In [33]:
epsilon, epsilon_i, N = sympy.symbols('epsilon, epsilon_i, N', nonnegative=True )

Will do the 3 layer version first


In [34]:
# vector of emissions
E = sympy.Matrix([E_0, E_1, E_2])
E


Out[34]:
$$\left[\begin{matrix}T_{0}^{4} \epsilon_{0} \sigma\\T_{1}^{4} \epsilon_{1} \sigma\\T_{2}^{4} \epsilon_{2} \sigma\end{matrix}\right]$$

In [35]:
# upwelling flux
fromsurface = E_s
U = sympy.Matrix([fromsurface, tau_0*fromsurface + E_0, tau_1*(tau_0*fromsurface + E_0) + E_1, 
                  tau_2*(tau_1*(tau_0*fromsurface + E_0) + E_1) + E_2])
U


Out[35]:
$$\left[\begin{matrix}T_{s}^{4} \sigma\\T_{0}^{4} \epsilon_{0} \sigma + T_{s}^{4} \sigma \left(- \epsilon_{0} + 1\right)\\T_{1}^{4} \epsilon_{1} \sigma + \left(- \epsilon_{1} + 1\right) \left(T_{0}^{4} \epsilon_{0} \sigma + T_{s}^{4} \sigma \left(- \epsilon_{0} + 1\right)\right)\\T_{2}^{4} \epsilon_{2} \sigma + \left(- \epsilon_{2} + 1\right) \left(T_{1}^{4} \epsilon_{1} \sigma + \left(- \epsilon_{1} + 1\right) \left(T_{0}^{4} \epsilon_{0} \sigma + T_{s}^{4} \sigma \left(- \epsilon_{0} + 1\right)\right)\right)\end{matrix}\right]$$

In [36]:
# downwelling flux...
fromspace = 0
D = sympy.Matrix([ tau_0*(tau_1*(tau_2*fromspace + E_2) + E_1) + E_0, 
                 tau_1*(tau_2*fromspace + E_2) + E_1, tau_2*fromspace + E_2, fromspace])
D


Out[36]:
$$\left[\begin{matrix}T_{0}^{4} \epsilon_{0} \sigma + \left(- \epsilon_{0} + 1\right) \left(T_{1}^{4} \epsilon_{1} \sigma + T_{2}^{4} \epsilon_{2} \sigma \left(- \epsilon_{1} + 1\right)\right)\\T_{1}^{4} \epsilon_{1} \sigma + T_{2}^{4} \epsilon_{2} \sigma \left(- \epsilon_{1} + 1\right)\\T_{2}^{4} \epsilon_{2} \sigma\\0\end{matrix}\right]$$

In [37]:
# Net flux, positive up
F = U - D
F


Out[37]:
$$\left[\begin{matrix}- T_{0}^{4} \epsilon_{0} \sigma + T_{s}^{4} \sigma - \left(- \epsilon_{0} + 1\right) \left(T_{1}^{4} \epsilon_{1} \sigma + T_{2}^{4} \epsilon_{2} \sigma \left(- \epsilon_{1} + 1\right)\right)\\T_{0}^{4} \epsilon_{0} \sigma - T_{1}^{4} \epsilon_{1} \sigma - T_{2}^{4} \epsilon_{2} \sigma \left(- \epsilon_{1} + 1\right) + T_{s}^{4} \sigma \left(- \epsilon_{0} + 1\right)\\T_{1}^{4} \epsilon_{1} \sigma - T_{2}^{4} \epsilon_{2} \sigma + \left(- \epsilon_{1} + 1\right) \left(T_{0}^{4} \epsilon_{0} \sigma + T_{s}^{4} \sigma \left(- \epsilon_{0} + 1\right)\right)\\T_{2}^{4} \epsilon_{2} \sigma + \left(- \epsilon_{2} + 1\right) \left(T_{1}^{4} \epsilon_{1} \sigma + \left(- \epsilon_{1} + 1\right) \left(T_{0}^{4} \epsilon_{0} \sigma + T_{s}^{4} \sigma \left(- \epsilon_{0} + 1\right)\right)\right)\end{matrix}\right]$$

In [38]:
# The absorption is then simply the flux convergence in each layer

# define a vector of absorbed radiation -- same size as emissions
A = E.copy()

# Get the convergence
for n in range(3):
    A[n] = -(F[n+1]-F[n])

A


Out[38]:
$$\left[\begin{matrix}- 2 T_{0}^{4} \epsilon_{0} \sigma + T_{1}^{4} \epsilon_{1} \sigma + T_{2}^{4} \epsilon_{2} \sigma \left(- \epsilon_{1} + 1\right) - T_{s}^{4} \sigma \left(- \epsilon_{0} + 1\right) + T_{s}^{4} \sigma - \left(- \epsilon_{0} + 1\right) \left(T_{1}^{4} \epsilon_{1} \sigma + T_{2}^{4} \epsilon_{2} \sigma \left(- \epsilon_{1} + 1\right)\right)\\T_{0}^{4} \epsilon_{0} \sigma - 2 T_{1}^{4} \epsilon_{1} \sigma - T_{2}^{4} \epsilon_{2} \sigma \left(- \epsilon_{1} + 1\right) + T_{2}^{4} \epsilon_{2} \sigma + T_{s}^{4} \sigma \left(- \epsilon_{0} + 1\right) - \left(- \epsilon_{1} + 1\right) \left(T_{0}^{4} \epsilon_{0} \sigma + T_{s}^{4} \sigma \left(- \epsilon_{0} + 1\right)\right)\\T_{1}^{4} \epsilon_{1} \sigma - 2 T_{2}^{4} \epsilon_{2} \sigma + \left(- \epsilon_{1} + 1\right) \left(T_{0}^{4} \epsilon_{0} \sigma + T_{s}^{4} \sigma \left(- \epsilon_{0} + 1\right)\right) - \left(- \epsilon_{2} + 1\right) \left(T_{1}^{4} \epsilon_{1} \sigma + \left(- \epsilon_{1} + 1\right) \left(T_{0}^{4} \epsilon_{0} \sigma + T_{s}^{4} \sigma \left(- \epsilon_{0} + 1\right)\right)\right)\end{matrix}\right]$$

In [39]:
# this should reduce to zero if I did it right
sympy.simplify(A - sympy.Matrix([R_0, R_1, R_2]))


Out[39]:
$$\left[\begin{matrix}0\\0\\0\end{matrix}\right]$$

So that works. I can formulate the tests against numerical code this way


In [ ]: