% 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
end
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
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
end
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))