WDE17_iceberg_model.m

% Iceberg Model -----------------------------------------------------------

% NOTE: To run this script you'll need the 5 .mat files that are

% found in this folder:

% https://www.dropbox.com/sh/9qc7l3yacyqgn3c/AAC69yLJUH39VsjbEwAdRRMOa?dl=0

% This model accompanies the paper "An Analytical Model of Iceberg Drift"

% Wagner, Dell, Eisenman, JPO (2017), henceforth WDE17

% Other articles making use of this model:

% Gone with the wind: How model biases skew iceberg meltwater distributions

% Wagner & Eisenman, GRL (2017).

% "On the representation of capsizing in iceberg models"

% Wagner, Stern, Dell, Eisenman, Ocean Model (2017)

% "Wave inhibition by sea ice enables trans-Atlantic ice rafting of debris

% during Heinrich Events"

% Wagner, Dell, Eisenman, Keeling, Severinghaus, EPSL, (submitted)

% The model computes the drift and decay of icebergs, by solving their

% motion analytically for given velocity and SST conditions. Those

% conditions are loaded first (see "input fields").

% Then model parameters and analytical expressions are set.

% Then run parameters (time, domain, space domain etc) are specified.

% Then the release locations are loaded (this is a .mat-file, constructed

% in "construct_seeding.m".

% Finally the different bergsizes are loaded and specified.

% Then we compute the Lagrangian trajectories, looping over every iceberg

% size class. These loops consist of 2 components for every timestep.

% 1) drifting - computes iceberg velocity and translation

% 2) melting - computes change in iceberg dimensions and volume

% Trajectories for every iceberg size class are saved as individual files

% in the output trajectory.

% Till Wagner & Ian Eisenman, May 2017,

% tjwagner@ucsd.edu, eisenman@ucsd.edu

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

% pick type of input ------------------------------------------------------

model = 'ECCO2'; %you can either run this with ECCO2 input fields or

% model = 'CCSM4'; %CCSM4 input fields

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

% load input fields -------------------------------------------------------

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

% ECCO2 data can be downloaded from

% ftp://ecco2.jpl.nasa.gov/data1/cube/cube92/lat_lon/quart_90S_90N/

% with JRA-25 surface wind reanalysis data

% (https://rda.ucar.edu/datasets/ds625.0/)

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

% CCSM4 data can be downloaded from

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

% the demo data set here is for the year 1992, and for the Atlantic only

% here we load surface ocean currents, winds, SST, and landmask from 1 file

load(strcat(model,'_1992_Atlantic.mat'))

fprintf('model data loaded \n')

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

% Iceberg Parameters and analytical expressions of

% gamma, Lambda, alpha, beta (eqns 7,8,9)

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

R = 6378*1e3; % earth radius in m

rhow = 1027; % density of water (kg/m^3)

rhoa = 1.2; % density of air (kg/m^3)

rhoi = 850; % density of shelf ice (kg/m^3), Silva et al (2006)

drho = rhow-rhoi;

Cw = 0.9; % bulk coefficient water (Bigg et al 1997)

Ca = 1.3; % bulk coefficient air (Bigg et al 1997)

Om = 7.2921*1e-5; % rotation rate of earth (rad/s)

ff = @(lat) 2*Om*sin(abs(lat)*pi/180); % latitude in degrees

ga = sqrt(rhoa*drho/rhow/rhoi*Ca/Cw); % gamma = sqrt(ca/cw), Eq. 7

S = @(l,w) l.*w./(l+w); % harmonic mean length

La = @(u,lat,S) Cw*ga/ff(lat).*u/S/pi; % Lambda, Eq. 9

% alpha and beta (eq 8 in WDE17)-------------------------------------------

a = @(La) 1./(2*La.^3).*(sqrt(1+4*La.^4)-1);

b = @(La) 1./(sqrt(2)*La.^3).*sqrt((1+La.^4).*sqrt(1+4*La.^4)-3*La.^4-1);

% Melt parameters (WDE17 Appendix)-----------------------------------------

Ti = -4;

a1 = 8.7e-6; a2 = 5.8e-7;

b1 = 8.8e-8; b2 = 1.5e-8;

c = 6.7e-6;

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

% specify the space domain ------------------------------------------------

LAT = double(input.latw); LON = double(input.lonw);

if strcmp(model,'CCSM4')

LAT2 = LAT(:); LON2 = LON(:);

vec = reshape(1:length(LON2),151,192);

end

minLAT = min(LAT(:)); maxLAT = max(LAT(:));

minLON = min(LON(:)); maxLON = max(LON(:));

msk = input.landmask;

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

% set run parameters ------------------------------------------------------

if strcmp(model,'ECCO2') % Input fields time step in days

DT = 3;

else

DT = 1;

end

trajnum = 25; % total number of iceberg trajectories to compute

final_t = 366/DT; % time span of input field

startrange = round(linspace(1,final_t/2,trajnum)); %evenly space seeding

dt = 24*3600; % model timestep in seconds (dt = 1 day)

dtR = dt/R*180/pi; % need this ratio for distances

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

t = 1:final_t; %input time

nt= length(t)*DT; %number of model timesteps

tt = linspace(1,length(t),nt); %model time

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

% Load Seeding fields -----------------------------------------------------

load(strcat(model,'_Laurentide_Seed'))

seed_X = repmat(Seed_X(:),[100,1]); %cycle through each location 100x

seed_Y = repmat(Seed_Y(:),[100,1]); %i.e. this can run 3600 icebergs

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

% these are the circulation fields-----------------------------------------

uwF = input.uw(:,:,t); vwF = input.vw(:,:,t); %water vels input

uaF = input.ua(:,:,t); vaF = input.va(:,:,t); %air vels input

sst = input.sst(:,:,t); %sst vels input

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

% loop over individual initial iceberg size classes -----------------------

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

bvec = 1:10; %vector of which size classes to compute - has to be [1,10]

bergdims=[100 67 67 % [L W H]

200 133 133

300 200 200

400 267 267

500 333 300

600 400 300

750 500 300

900 600 300

1200 800 300

1500 1000 300];

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

for bb = bvec

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

bergsize = bb; % current berg size class

fprintf('run bergsize B%d \n',bergsize)

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

XIL = nan(trajnum,nt); YIL = XIL; VOL = XIL; %set output arrays

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

L = bergdims(bergsize,1); %initialize iceberg dimensions

W = bergdims(bergsize,2);

H = bergdims(bergsize,3);

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

% run drift and melt---------------------------------------------------

mm=0; ss=0; ob=0; %set counters for melted, survived, and escaped bergs

for j = 1:trajnum

if mod(j,10)==0; fprintf('%d trajectories computed \n',j); end

ts = startrange(j); %pick seeding time (of input field)

tts= ts*DT; %trajectory start time (of model)

lt = nt-tts; %trajectory run length

xil = nan(1,lt); yil = xil; v = xil; %initialize output vectors

yig = seed_Y(j); xig = seed_X(j); % cycle through start locations

if strcmp(model,'ECCO2')

xil(1) = LON(xig); yil(1) = LAT(yig); %initial lon and lat

else

xil(1) = LON(yig,xig); yil(1) = LAT(yig,xig);

end

l = L*ones(1,lt); w = l*W/L; h = l*H/L; %initial berg dimensions

v(1) = L*W*H; %initial volume and dvol

% integrate as long as the iceberg is in the domain and not

% melted and over the time period specified above

i = 0; outofbound = 0; melted = 0;

while outofbound == 0 && melted == 0 && i<lt-1

i = i+1;

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

% Computes iceberg drift component, WDE17, Section 3

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

% find nearest neighbor of surface condition field

if strcmp(model,'ECCO2') % this only works on a regular grid

YI = dsearchn(LAT,yil(i));

XI = dsearchn(LON,xil(i));

else %CCSM4 is on an irregular grid

indic = dsearchn([LAT2,LON2],[yil(i),xil(i)]);

[XI,YI]=find(indic==vec);

end

% now interpolate fields linearly between timesteps------------

timestep = tt(tts+i);

t1 = floor(timestep); t2 = t1+1;

dt1 = timestep-t1; dt2 = t2-timestep;

ua = uaF(XI,YI,t1)*dt1+uaF(XI,YI,t2)*dt2;

va = vaF(XI,YI,t1)*dt1+vaF(XI,YI,t2)*dt2;

uw = uwF(XI,YI,t1)*dt1+uwF(XI,YI,t2)*dt2;

vw = vwF(XI,YI,t1)*dt1+vwF(XI,YI,t2)*dt2;

SST= sst(XI,YI,t1)*dt1+sst(XI,YI,t2)*dt2;

% compute wind speed and Lambda at location (for given iceberg

% size)--------------------------------------------------------

Ua = sqrt(ua^2+va^2);

LA = La(Ua,yil(i),S(l(i),w(i)));

% now compute iceberg velocity Eq (6)--------------------------

ui = uw + ga*( a(LA)*va + b(LA)*ua);

vi = vw + ga*(-a(LA)*ua + b(LA)*va);

% iceberg translation (convert from m to deg lat/lon)----------

dlon = ui*dtR;

dlat = vi*dtR;

yil(i+1) = yil(i) + dlat;

xil(i+1) = xil(i) + dlon/cos((yil(i+1)+yil(i))/2*pi/180);

% check you haven't gone out of bounds-------------------------

if xil(i+1)>maxLON || xil(i+1)<minLON || ...

yil(i+1)>maxLAT || yil(i+1)<minLAT

outofbound = 1;

ob = ob+1;

fprintf('iceberg %d left domain at timestep %d \n',j,i);

else % now check iceberg not on land---------------------------

if strcmp(model,'ECCO2')

yi2(1) = find(LAT<=yil(i+1),1,'last');

yi2(2) = find(LAT>yil(i+1),1,'first');

xi2(1) = find(LON<=xil(i+1),1,'last');

xi2(2) = find(LON>xil(i+1),1,'first');

% when new position within one grid box of land:

if any(find(msk(xi2,yi2)==0))

yil(i+1) = yil(i); %iceberg reset to time step i

xil(i+1) = xil(i);

end

else % similar for CCSM4

indic = dsearchn([LAT2,LON2],[yil(i+1),xil(i+1)]);

[XI,YI]=find(indic==vec);

if input.landmask(XI,YI)==1

yil(i+1) = yil(i);

xil(i+1) = xil(i);

end

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

% Here we compute the iceberg melting component, WDE15 Appendix

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

% Compute Melt Terms Eq A1 in WDE17----------------------------

Me = a1*Ua^0.5 + a2*Ua;

Mv = b1*SST + b2*SST^2;

Mb = c*sqrt((ui - uw).^2+(vi-vw).^2)^0.8*(SST-Ti)*l(i)^(-0.2);

% Apply Melt Rates --------------------------------------------

dldt = - Mv - Me; dhdt = - Mb;

l(i+1) = l(i)+dldt*dt;

w(i+1) = w(i)+dldt*dt;

h(i+1) = h(i)+dhdt*dt;

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

% Make sure the berg is not negative size

if l(i+1)<0 || w(i+1)<0 || h(i+1)<0

l(i+1)=0; w(i+1)=0; h(i+1)=0;

melted = 1;

mm = mm+1;

end

% Rollovers (Eq A2 in WDE17)

if w(i+1)/h(i+1) < sqrt(6*rhoi/rhow*(1-rhoi/rhow))

hn = w(i+1); w(i+1) = h(i+1); h(i+1) = hn;

end

% Make sure length>width:

if w(i+1)>l(i+1)

wn = l(i+1); l(i+1)=w(i+1); w(i+1) = wn;

end

% Compute new volume

v(i+1) = l(i+1)*w(i+1)*h(i+1);

% Check whether iceberg survived-------------------------------

if i == lt-1 && v(i+1) > 0

ss = ss+1;

end

ind = 1:i+1;

XIL(j,ind)=xil(ind); YIL(j,ind)=yil(ind); %store trajectory

VOL(j,ind)=v(ind); %store volume

end

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

fprintf('%d icebergs died, %d lived, %d left the domain \n',mm,ss,ob)

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

save(sprintf('%s_Size%d',model,bb),'XIL','YIL','VOL'); %save output

end

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

% Plot Iceberg Trajectories------------------------------------------------

% This loads in output files from "iceberg_model.m" for different iceberg

% sizes and plots them (colorcoded).

% Requires the matlab mapping toolbox. I made the atl.mat file to quickly

% plot a low-resolution landmask of the north atlantic.

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

figure(1); clf

% make North Atlantic Map plot (requires mapping toolbox)------------------

load atl

latlim = [33 70];

lonlim = [-80 0];

worldmap(latlim,lonlim)

gri = gridm('off');

geo = geoshow(atl(1:10), 'FaceColor', [.9 .9 .9]);

setm(gca,'FFaceColor',[1 1 1])

hold on

hm = mlabel('on');

set(hm(:),'fontsize',13);

set(hm(:),'verticalalignment','baseline');

pm = plabel('on');

set(pm(:),'fontsize',13);

set(pm(:),'verticalalignment','baseline');

%uncomment the following if you've precomputed output files and -----------

%want to call just the plotting cell --------------------------------------

% model = 'ECCO2';

N = 25; %number of trajectories to be plotted

bvec = [9 7 5 3 1]; %specify different iceberg sizes you're interested in

p = cell(1,length(bvec)); pl = p; leg = p;

ii = 0;

for i = bvec

load(sprintf('%s_Size%d',model,i))

bb = i; bs = bb;

XL = XIL;

YL = YIL;

ind = 1:N;

ii = ii+1;

for tind = 1:length(ind)

vend = find(VOL(ind(tind),:)<.1*VOL(ind(tind),1),1,'first');

%plot until iceberg is 90% decayed

if isempty(vend)

vend = size(VOL,2);

end

if tind == 1

p{ii} = plotm(YL(ind(1),1:vend)',XL(ind(1),1:vend)'...

,'-','col',1-[bb/15 1-bb/15 1-bb/15]);

else

plotm(YL(ind(tind),1:vend)',XL(ind(tind),1:vend)'...

,'-','col',1-[bb/15 1-bb/15 1-bb/15]);

end

for i = 1:ii; pp = p{i}; pl{i}=pp(1); leg{i} = sprintf('%d',bvec(i)); end

ll = legend([pl{:}],leg(:),'location','northeastoutside',...

'interpreter','latex');

legendtitle(ll,'Size Class','fontweight','normal');

title(sprintf('%s Forcing',model))