ABSTRACT (original article): The temperature in the transient climate response is lower than the equilibrium temperature for the same amount of forcing. The degree of disequilibrium is not constant in time and depends on various climate parameters. We derive intuition for this by solving the heat equation with a surface temperature feedback for linearly increasing forcing. The surface temperature initially evolves at a slower rate than the corresponding steady state (SS) temperature and it accelerates until quasi-steady state (QSS), when the SS and QSS temperatures evolve in parallel with a constant offset. The offset depends on the rate of forcing and total heat capacity of the system divided by the square of the climate feedback. The timescale over which the climate system approaches QSS depends also on the effective ocean mixing and is order of thousands of years. Over societally relevant timescales (around 100 years), the top-of-atmosphere energy imbalance increases, and the actual temperature moves farther from the steady-state temperature expected for the same forcing. CITATION (original article): Chtirkova, B. (2026). Analytical insights into the transient climate response. arXiv:2603.01674. https://doi.org/10.48550/arXiv.2603.01674
The problem of climate change can be modelled mathematically as the radiative response to a perturbed energy input. The perturbation of the energy input is referred to as forcing (F
[
-2
Wm
]
) and the system response can be described by a change in the surface temperature (Ts [K]) and corresponding change in the emission to space (λ
T
s
[
-2
Wm
]
). In steady state, the flux perturbation should be balanced by the change in surface temperature needed to reach the new flux equilibrium:​F -
λ
T
s
= 0 flux equilibrium ​
T
s
= F/λ surface temperature required to reach flux equilibrium​The climate feedback λ[Wm-2K-1] is comprised of a change in the Plank emissivity, increased water vapour content, etc. It is assumed constant and the temperature response to forcing is assumed linear for the small ranges of interest. An Earth’s atmosphere representative value is λ = 1
[
-2
Wm
-1
K
]
. A doubling of carbon dioxide concentrations amounts to a forcing of around F = 4
[
-2
Wm
]
.​To gain intuition for the temperature change over time, the effective heat capacity of the system needs to be accounted for. A simple model for this is a 1-dimensional diffusive domain, in which the heat conductivity coefficient is would account for the ocean processes which mix the warmed waters in the upper layers with those below. We assume mixing is proportional to the temperature difference and the heat conductivity parameter K
[
-1
WK
-1
m
]
is assumed constant in depth. An Earth’s oceans representative value is K = 350
-1
WK
-1
m
with a heat capacity of seawater of c =
6
4.10
-3
Jm
-1
K
.​One can set the 1D heat equation with upper Robin and lower Neumann boundary condition which describe the forcing-feedback energy input at the top and ensure energy conservation by no flux leakage at the bottom. We want to solve the problem in which the forcing increases linearly with time F = at. For present-day, a representative value of the forcing rate is a = 0.03
-2
Wm
/year.​​
c
∂T
∂t
K
2
∂
T
∂
2
z
​​Boundary conditions:​
z0:F(t,z0)-K
∂T
∂z
z0
F-λT(t,z0),λ>0​​zL:F(t,zL)-K
∂T
∂z
zL
0
​​And we assume an initial condition of unperturbed steady state:​
t=0:
T(t0,z)0
Clear
ClearAll[T,t,z,a,λ,c,K,L]
Assumptions for parameters, flux is positive downwards convention
$Assumptions={z>=0,a>0,λ>0,c>0,K>0,L>0};
Define the diffusion equation
eqn=c*D[T[t,z],t]-K*D[T[t,z],{z,2}];
Boundary conditions
bc1=(-K*D[T[t,z],z]-a*t+λ*T[t,z])/.z->0;(*Robin*)​​bc2=(-K*D[T[t,z],z])/.z->L;(*Neumann*)
Initial condition
ic=T[t,z]/.t->0;
Create the full system
eqnbc={eqn==0,bc1==0,bc2==0};​​eqnbcic={eqn==0,bc1==0,bc2==0,ic==0};
Since we constantly increase the energy input of the system, we do not expect it to approach a steady state in which dT/dt = 0. However, a solution of quasi-steady state in which dT/dt = const exists. In this solution, the rate of change of the temperature is a/λ and the depth profile T(z) is quadratic - the surface is warmer than the deeper parts of the domain.
In[]:=
Print["Quasi-steady state (QSS) solution, T(t, z)"]​​sol=DSolve[eqnbc,T,{t,z}];​​solbc=T[t,z]/.First[First[sol]]​​​​Print["QSS Surface temperature, T(t, z=0)"]​​tsqss=FullSimplify[solbc/.z->0]​​​​Print["QSS quadratic depth profile, T(z)"]​​Tzqss=solbc/.t->0​​​​(*Plot*)​​a=0.1;λ=1;c=1;K=1;L=10;​​Plot3D[solbc,{t,0,5},{z,0,10},AxesLabel{"t","z","T(t,z)"}]​​Clear[a,λ,c,K,L]
Quasi-steady state (QSS) solution, T(t, z)
Out[]=
-2acKL+2aKtλ-2acLzλ+ac
2
z
λ
2K
2
λ
QSS Surface temperature, T(t, z=0)
Out[]=
a(-cL+tλ)
2
λ
QSS quadratic depth profile, T(z)
Out[]=
-2acKL-2acLzλ+ac
2
z
λ
2K
2
λ
Out[]=
If we are interested in the path the system takes to reach this quasi-steady state, we need to solve for a system with an initial condition in equilibrium. In this case, the true analytical solution is an infinite sum of decaying exponentials for their own spatial chunks, the coefficients of which are the solutions of the transcendental equation tan(Lx) = λ/(Kx).
Solution with relaxation to QSS
Print["Full solution with relaxation to QSS"]​​sol=Assuming[a>0,DSolve[eqnbcic,T,{t,z}]];​​solbcic=T[t,z]/.First[First[sol]]
Full solution with relaxation to QSS
Out[]=
at
λ
+
∞
∑
K[1]=1
-2
2
ac
-
KtK[2,K[1]]
c

-1+
KtK[2,K[1]]
c

λ-λCosL
K[2,K[1]]
+K
K[2,K[1]]
SinL
K[2,K[1]]
Cosz
K[2,K[1]]

K[2,K[1]]
+
λSinz
K[2,K[1]]

K

K
λ
K[2,K[1]]
-
3
λ
2
CosL
K[2,K[1]]

3
K
+
2λK[2,K[1]]
K
+
L
2
λ
K[2,K[1]]
2
K
-
λ
2
CosL
K[2,K[1]]

K[2,K[1]]
K
+L
2
K[2,K[1]]
-2
3
λ
2
CosL
K[2,K[1]]

+Kλ3K+2Lλ-KCos2L
K[2,K[1]]
K[2,K[1]]+2
3
K
L
2
K[2,K[1]]

TanL
K[2,K[1]]

λ
K
K[2,K[1]]
&&K[1]∈&&K[1]≥1&&K[2,K[1]]>0
Indeterminate
True
For Earth’s climate, we can assume that the depth of the ocean if very large, hence the approximation Lλ >> K allows us to find the first root of the transcendental equations:
2
(π/(2L))
, which corresponds to the slowest decaying mode of the system.
Approximation 1
Print["Approximation 1: slowest decay mode. Extract only the first term from the infinite sum across coefficients and replace K[2, K[1]] = pi^2/(4L^2)."]​​​​solbcicsterm1=
at
λ
-2
2
ac
-
KtK[2,K[1]]
c

-1+
KtK[2,K[1]]
c

λ-λCosL
K[2,K[1]]
+K
K[2,K[1]]
SinL
K[2,K[1]]
Cosz
K[2,K[1]]

K[2,K[1]]
+
λSinz
K[2,K[1]]

K

K
λ
K[2,K[1]]
-
3
λ
2
CosL
K[2,K[1]]

3
K
+
2λK[2,K[1]]
K
+
L
2
λ
K[2,K[1]]
2
K
-
λ
2
CosL
K[2,K[1]]

K[2,K[1]]
K
+L
2
K[2,K[1]]
-2
3
λ
2
CosL
K[2,K[1]]

+Kλ3K+2Lλ-KCos2L
K[2,K[1]]
K[2,K[1]]+2
3
K
L
2
K[2,K[1]]
;​​solicbc1=solbcicsterm1/.K[2,K[1]]->(Pi^2)/(4*L^2)
Plot
A second approximation would be that the slowest decay mode has the spatial chunk of the entire domain, i.e. the quadratic profile from the QSS solution.
Approximation 2
Plot
Redefine analytical expressions for the surface temperature and imbalance
Parameters representative for the climate model
Evaluate Ts, N and their derivatives
To an order of magnitude, our analytical expectations for the surface temperature, system imbalance and their tendencies, agree with the climate model simulations which explicitly resolve the large-scale atmosphere-ocean fluid dynamics. For global mean scales, the climate response to forcing can be reasonably well approximated to that of a linear system with a surface feedback and a diffusive domain.

CITE THIS NOTEBOOK

Analytical insights into the transient climate response​
by Boriana Chtirkova​
Wolfram Community, STAFF PICKS, May 15, 2026
https://community.wolfram.com/groups/-/m/t/3716486