Modeling the Greenland Ice Sheet with IceBridge data
Goals
Go to IntroductionTutorial steps to be taken:
MeshWe modify the experiment from the Greenland SeaRISE tutorial, and improve from there. Run the first step in >> plotmodel(md,'data','mesh');
Now, we want to refine the mesh in JI area. An outline of this area >> exptool('Jak_outline.exp');
Next, we modify the Use MATLAB's zoom tool in the figure to make a close-up of the JI domain. ParameterizationWe want to include high-resolution bedrock and surface elevation data acquired in the OIB mission. The data is accessible at \href{http://data.cresis.ku.edu/data/grids/Jakobshavn_2008_2011_Composite_XYZGrid.txt}. Save the file in the To do this, the bedrock data is read, transformed into a usable grid, and interpolated to the mesh in the parameter file %Reading IceBridge data for Jakobshavn
disp(' reading IceBridge Jakobshavn bedrock');
fid = fopen('../Data/Jakobshavn_2008_2011_Composite/grids/Jakobshavn_2008_2011_Composite_XYZGrid.txt');
titles = fgets(fid);
data = fscanf(fid,'%g,%g,%g,%g,%g',[5 266400])';
fclose(fid);
[xi,yi]= ll2xy(md.mesh.lat,md.mesh.long,+1,45,70);
bed = flipud(reshape(data(:,5),[360 740])); bed(find(bed==-9999))=NaN;
bedy = flipud(reshape(data(:,1),[360 740]));
bedx = flipud(reshape(data(:,2),[360 740]));
%Insert Icebridge bed and recalculate thickness
bed_jks=InterpFromGridToMesh(bedx(1,:)',bedy(:,1),bed,xi,yi,NaN);
in=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,...
'Jak_grounded.exp','node',1);
bed_jks(~in)=NaN;
pos=find(~isnan(bed_jks));
md.geometry.base(pos)=bed_jks(pos);
md.geometry.thickness=md.geometry.surface-md.geometry.base;
Modify the Solution: %Reading IceBridge data for Jakobshavn
disp(' reading IceBridge Jakobshavn bedrock');
fid = fopen('../Data/Jakobshavn_2008_2011_Composite_XYZGrid.txt');
titles = fgets(fid);
data = fscanf(fid,'%g,%g,%g,%g,%g',[5 266400])';
fclose(fid);
[xi,yi]= ll2xy(md.mesh.lat,md.mesh.long,+1,45,70);
bed = flipud(reshape(data(:,5),[360 740])); bed(find(bed==-9999))=NaN;
surf = flipud(reshape(data(:,4),[360 740])); surf(find(surf==-9999))=NaN;
bedy = flipud(reshape(data(:,1),[360 740]));
bedx = flipud(reshape(data(:,2),[360 740]));
%Insert Icebridge bed and recalculate thickness
bed_jks=InterpFromGridToMesh(bedx(1,:)',bedy(:,1),bed,xi,yi,NaN);
surf_jks=InterpFromGridToMesh(bedx(1,:)',bedy(:,1),surf,xi,yi,NaN);
in=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,...
'Jak_grounded.exp','node',1);
bed_jks(~in)=NaN;
surf_jks(~in)=NaN;
pos=find(~isnan(bed_jks));
md.geometry.base(pos)=bed_jks(pos);
md.geometry.surface(pos)=surf_jks(pos);
md.geometry.thickness=md.geometry.surface-md.geometry.base;
Next, let's plot the surface elevation, the ice thickness, and base:
To plot the difference in the ice base topography between SeaRISE and OIB datasets do (1) modify the parameterization step in your >> md2=loadmodel('Models/Greenland.Parameterization2')
>> md=loadmodel('Models/Greenland.Parameterization')
>> plotmodel(md,'data',md2.geometry.base-md.geometry.base)
Zoom to the JI basin for better visibility. Stress BalanceWe now use inverse control methods to solve for Greenland friction coefficient. The velocity map below contains some gaps. Exclude the gaps from the inversion by creating a new >> exptool('data_gaps.exp')
Exclude these data gaps in the inversion by giving them zero weight during the inversion process: in=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y, 'data_gaps.exp','node',1);
md.inversion.cost_functions_coefficients(find(in),1)=0.0;
md.inversion.cost_functions_coefficients(find(in),2)=0.0;
Launch the stressbalance simulation, and plot velocity and basal friction coefficient. A logarithmic plot scale reveals more highlights of the velocity field structure: >> plotmodel(md,'data',md.results.StressbalanceSolution.Vel,'log',10,'caxis',[0.5 5000]);
>> plotmodel(md,'data',md.results.StressbalanceSolution.FrictionCoefficient);
They should look like this: Even at this coarse resolution we can identify the high friction values inland and lower values towards the coast, which may be related to the basal thermal regime of the ice sheet. TransientFinally, do a transient run (step 4) for 20 years, and decrease the surface mass balance linearly by 1 m w.e./yr over the last 10 years ( %Set surface mass balance
x1 = ncread(ncdata,'x1');
y1 = ncread(ncdata,'y1');
smb = ncread(ncdata,'smb');
smb = InterpFromGridToMesh(x1,y1,smb',md.mesh.x,md.mesh.y,0)*1000/md.materials.rho_ice;
smb = [smb smb smb-1.0];
md.smb.mass_balance = [smb;1 10 20];
Your results will be located in You can plot time series of surface mass balance, mean velocity and ice volume: ResultsWell done! Here are some suggestions on what to explore further:
|