WE15_EBM_simple.m

% Reference: "How Model Complexity Influences Sea Ice Stability",

% T.J.W. Wagner & I. Eisenman, J Clim (2015)

% WE15_EBM_simple.m:

% This code describes the EBM as discussed in Sec. 2b of the article above,

% hereafter WE15. Here we use central difference spatial integration and

% time stepping with MATLAB's ode45.

% The code WE15_EBM_fast.m, on the other hand, uses a faster, but more

% complicated formulation of the diffusion operator and Implicit Euler time

% stepping.

% Parameters are as described in WE15, table 1. Note that we do not include

% ocean heat flux convergence or a seasonal cylce in the forcing

% (equivalent to S_1 = F_b = 0 in WE15). This code uses an ice albedo when

% T<0 (WE15 instead uses the condition E<0, which is appropriate for the

% inclusion of a seasonal cycle in ice thickness). In this code, we define

% T = Ts - Tm, where Ts is the surface temperature and Tm the melting point

% (WE15, by contrast, defines T = Ts).

% Till Wagner & Ian Eisenman, Mar 15

% tjwagner@ucsd.edu or eisenman@ucsd.edu

function [x_out,T_out] = WE15_EBM_simple

%%Model parameters (WE15, Table 1 and Section 2d) -------------------------

global D B n cw aw ai S A F x dx

D = 0.6; % diffusivity for heat transport (W m^-2 K^-1)

A = 193; % OLR when T = T_m (W m^-2)

B = 2.1; % OLR temperature dependence (W m^-2 K^-1)

cw = 9.8; % ocean mixed layer heat capacity (W yr m^-2 K^-1)

S0 = 420; % insolation at equator (W m^-2)

S2 = 240; % insolation spatial dependence (W m^-2)

a0 = 0.7; % ice-free co-albedo at equator

a2 = 0.1; % ice=free co-albedo spatial dependence

ai = 0.4; % co-albedo where there is sea ice

F = 0; % radiative forcing (W m^-2)

% -------------------------------------------------------------------------

n = 50; % grid resolution (number of points between equator and pole)

x = linspace(0,1,n);

dx = 1/(n-1);

S = S0-S2*x.^2; % insolation [eq. (3) in WE15, with S_1 = 0]

aw = a0-a2*x.^2; % open water albedo

T0 = 10*ones(1,n); % initial condition (constant temp. 10C everywhere)

tspan = [0 30]; % timespan in years [t_start t_end]

% time integration---------------------------------------------------------

[t,T] = ode45(@rhsWE15,tspan,T0);

if nargout==0 % plot results-----------------------------------------------

figure(1), clf

subplot(1,2,1)

plot(t,T)

xlabel('$t$ (years)');

ylabel('$T$ ($^\circ$C, for different latitudes)')

xlim([0 t(end)])

box on

subplot(1,2,2)

hold all

plot(x,T(end,:), 's-')

xlabel('$x = sin \theta$');

ylabel('$T$ ($^\circ$C)')

box on

else % save results--------------------------------------------------------

x_out=x;

T_out=T;

end

% -------------------------------------------------------------------------

% ODE with spatial finite differencing for ode45 --------------------------

function dTdt = rhsWE15(~,T)

global D B n cw aw ai S A F x dx

Tdot = zeros(n,1);

% forcing

alpha = aw.*(T>0)'+ai*(T<0)'; % WE15, eq. (4)

C = alpha.*S-A+F;

% solve c_wdT/dt = D(1-x^2)d^2T/dx^2 - 2xDdT/dx + C - BT [cf. WE15, eq.(2)]

% use central difference

Tdot(2:n-1) = (D/dx^2)*(1-x(2:n-1).^2)'.*(T(3:n)-2*T(2:n-1)+T(1:n-2)) ...

-(D*x(2:n-1)/dx)'.*(T(3:n)-T(1:n-2));

% boundary condition: let dT/dt(x=0) = 0

Tdot(1) = D*(2*(T(2)-T(1))/dx^2); % WE15, eq.(7)

% use asymmetric (backward) spatial difference for the pole

Tdot(n) = - D*(2*x(n)*(T(n)-T(n-1))/dx); %O(dt) accuracy at x=1

% temperature tendency

dTdt = (Tdot+C' - B*T)/cw;