Hi,
I am trying to run two transient models of Thwaites Glacier with two different bed topography datasets (Bedmap2 and BedMachine). I can get both transient runs working, but the results of the Bedmap2 run show almost no change of thickness, velocity or grounding line position with time, whereas the results of the BedMachine run show large changes in these fields through time. The only difference in the code for the two models is where I bring in the bed datasets to build the geometry, and both models initialize correctly, so I'm not sure where to look for the cause of this error. The models are based on the PIG example tutorial. Any suggestions of where I could look to fix this issue? Thanks!
Code to load in Bedmap2 and build geometry:
%Binary Loading
disp(' Loading Bedmap2 data from Binary');
fid = fopen('../Data/bedmap2_bin/bedmap2_bed.flt','r','l');
topg = fread(fid,[6667,6667],'float32');
fclose(fid);
fid = fopen('../Data/bedmap2_bin/bedmap2_surface.flt','r','l');
usrf = fread(fid,[6667,6667],'float32');
fclose(fid);
fid = fopen('../Data/bedmap2_bin/bedmap2_thickness.flt','r','l');
thick = fread(fid,[6667,6667],'float32');
fclose(fid);
%read coordinates
ncols =6667;
nrows =6667;
xll =-3333500;
yll =-3333500;
yllmax =3333500;
gridsize =1000;
x_m=xll+(0:1:ncols)'*gridsize;
y_m=(yllmax-nrows*gridsize)+(0:1:nrows)'*gridsize;
%Geometry
disp(' Interpolating surface and ice base');
md.geometry.bed = InterpFromGridToMesh(x_m,y_m,rot90(topg),md.mesh.x,md.mesh.y,0);
md.geometry.surface = InterpFromGridToMesh(x_m,y_m,rot90(usrf),md.mesh.x,md.mesh.y,0);
clear usrf, topg;
%solve for thickness
md.geometry.thickness=md.geometry.surface-md.geometry.bed;
%solve for floating ice thickness
HB=md.geometry.surface*(md.materials.rho_water/(md.materials.rho_water-md.materials.rho_ice));
md.geometry.thickness(md.geometry.thickness>HB)=HB(md.geometry.thickness>HB);
md.geometry.base=md.geometry.surface-md.geometry.thickness;
%ensure hydrostatic equilibrium
disp(' Ensuring hydrostatic equilibrium on ice shelf')
md.geometry.hydrostatic_ratio=ones(md.mesh.numberofvertices,1);
% Set min thickness to 1 meter % {{{
disp(' Min thickness is set to 1 meter')
pos = find(md.geometry.thickness <= 1);
md.geometry.thickness(pos) = 1;
md.geometry.surface(pos)=md.geometry.surface(pos)+1;
md.geometry.base=md.geometry.surface-md.geometry.thickness;
md.mask.groundedice_levelset(HB==md.geometry.thickness)=-1;
% md.mask.groundedice_levelset(md.geometry.bed< md.geometry.base+0.01)=-1;
md.inversion.thickness_obs = md.geometry.thickness;`
Code to load in Bedmachine and build geometry:
`% Get bedmachine data to build up the surface topography grid
surface= double(ncread('../Data/BedMachineAntarctica-2019-05-14-002.nc','surface'));
bed= double(ncread('../Data/BedMachineAntarctica-2019-05-14-002.nc','bed'));
ysurf= double(ncread('../Data/BedMachineAntarctica-2019-05-14-002.nc','y'));
xsurf= double(ncread('../Data/BedMachineAntarctica-2019-05-14-002.nc','x'));
y_surf=flipud(ysurf);
x_surf=(xsurf);
surfflip=flipud(surface');
bedflip=flipud(bed');
%interpolate bed data onto our mesh vertices
md.geometry.surface=InterpFromGridToMesh(x_surf,y_surf,surfflip,md.mesh.x,md.mesh.y,0);
md.geometry.bed=InterpFromGridToMesh(x_surf,y_surf,bedflip,md.mesh.x,md.mesh.y,0);
%solve for thickness
md.geometry.thickness=md.geometry.surface-md.geometry.bed;
%solve for floating ice thickness
%HB is flotation criterion
HB=md.geometry.surface*(md.materials.rho_water/(md.materials.rho_water-md.materials.rho_ice));
md.geometry.thickness(md.geometry.thickness>HB)=HB(md.geometry.thickness>HB);
md.geometry.base=md.geometry.surface-md.geometry.thickness;
%ensure hydrostatic equilibrium
disp(' Ensuring hydrostatic equilibrium on ice shelf')
md.geometry.hydrostatic_ratio=ones(md.mesh.numberofvertices,1);
% Set min thickness to 1 meter % {{{
disp(' Min thickness is set to 1 meter')
pos = find(md.geometry.thickness <= 1);
md.geometry.thickness(pos) = 1;
%adjust surface according to adjusted thicknesses
md.geometry.surface(pos)=md.geometry.surface(pos)+1;
md.geometry.base=md.geometry.surface-md.geometry.thickness;
md.mask.groundedice_levelset(HB==md.geometry.thickness)=-1;
% md.mask.groundedice_levelset(md.geometry.bed< md.geometry.base+0.01)=-1;
%Double Check geometry for consistency
pos=find(md.mask.groundedice_levelset>=0);
md.geometry.base(pos)=md.geometry.bed(pos);
md.geometry.thickness(pos)=md.geometry.surface(pos)-md.geometry.base(pos);
md.inversion.thickness_obs = md.geometry.thickness;