ISSM-SIMIP-HOM（高阶和全斯托克斯冰盖模型）的示例程序如网站issm.jpl.nasa.gov。
本研究选用的是以Matlab代码的形式驱动。我将以官网上的ismipF.par参数文件为例进行演示。
在无气候输入情况下，驱动ISMIP-HOM的方法依次是——
（1）建立三维网格。输入网格的东西向长度、横向节点数量；南北向长度，纵向节点数量；选择模型层数。
（2）地形初始化。输入模型表面（DEM）、冰川厚度分布、冰川基底（DEM-冰川厚度）。需要注意的是，在模型表面的任意一处，冰川厚度都不能定义为0，这会导致模型不连续，所以我的建议是，冰川厚度为0的区域，都取值为0.0000001。
（3）输入冰流系数，选择模型演化时间。
（4）在runme.m 与ismipF.par文件夹下面建立一个Models文件夹（否则程序会报错）。
（5）准备好数据DEM、冰川厚度、冰川表面速度长（请重采样到一致的行列数）。如果你在网格生成那一步选择了88x99网格，那么我建议你的DEM、冰川厚度、冰川表面速度也重采样到88x99。
（6）运行runme.m。
下面这个是依据冰川基底剪应力，结合冰川表面地形和DEM反演冰川厚度的代码，运行前请准备好SLOPE.tif、2013_V_Final.tif、DEM.tif——
%% SHORT SCRIPT TO NUMERICALLY SOLVE ICE FLOW/SLIDING EQUATIONS AND INVERT FOR ICE THICKNESS
%% inputs
slope = imread("SLOPE.tif");
velocity = imread("2013_V_Final.tif");
DEM = imread("DEM.tif");
n = 3; %Glen's flow constant
A_c = 2.5e-24; %Arrhenius creep constant
f = 0.9; %lateral drag correction ('shape factor')
p_i = 917; %density of ice
g = 9.79; %gravity
spery = 60*60*24*365; %seconds per year
C1_1 = 0.0012/spery; %0.0012 Sliding constant
C2_1 = 0; %Basal water fraction
% threshold_1 = 6000; %if threshold added
%
% C1_2 = 0.0005/spery; %0.0012
%
% C2_2 = 0.1;
%
% threshold_2 = 5500;
%
% C1_3 = 0.0010/spery; %0.0012
%
% C2_3 = 0.25;
% %Smooth slope along 150m coupling length
% slope = movmean(slope,3,2);
% slope = movmean(slope,3,1);
% %Make slope clip same size as velocity matrix if necessary
% Slope_clip = interp2(Slope_clip, linspace(1, size(Slope_clip,2), size(Ice_velocity_resamples,2)).', linspace(1, size(Slope_clip,1), size(Ice_velocity_resamples,1)));
% Reduce size to speed up calculations if necessary
% slope = interp2(slope, linspace(1, size(slope,2), size(slope,2)/5).', linspace(1, size(slope,1), size(slope,1)/5));
% velocity = interp2(velocity, linspace(1, size(velocity,2), size(velocity,2)/5).', linspace(1, size(velocity,1), size(velocity,1)/5))
% save_volume = NaN(5,5); %If iterating
% save_summit = NaN(5,5);
% % iterate C1 and C2
% C1list = [0.0001,0.0005,0.001,0.002];
% for C1loop = 1:length(C1list)
% C2list = [0,0.1,0.25,0.5,0.95];
% for C2loop = 1:length(C2list)
%% Solve
slope_long=deg2rad(reshape(slope,[size(slope,1)*size(slope,2),1]));
velocity_long=reshape(velocity,[size(velocity,1)*size(velocity,2),1])/spery;
DEM_long=reshape(DEM,[size(DEM,1)*size(DEM,2),1]);
storage_mat = NaN(size(velocity_long,1),1);
% curr_pix = zeros(size(velocity_long,1)*4,1);
tic
for loop = 1:size(velocity_long,1)
if ~isnan(velocity_long(loop,1))
% if DEM_long(loop,1) > threshold_1 %If threshold added
syms H
S = vpasolve((H^(n+1)*2*A_c*(f*p_i*g*sin(slope_long(loop,1)))^n/(n+1))...
+ (H^2*sin(slope_long(loop,1))^2*((C1_1*(f*p_i*g)^2)/((1-C2_1)*(H*p_i*g))))...
== velocity_long(loop,1), H);
storage_mat(loop,1) = S(2,1);
% elseif DEM_long(loop,1) <= threshold_1 && DEM_long(loop,1) > threshold_2
% syms H
%
% S = vpasolve((H^(n+1)*2*A_c*(f*p_i*g*sin(slope_long(loop,1)))^n/(n+1))...
% + (H^2*sin(slope_long(loop,1))^2*((C1_2*(f*p_i*g)^2)/((1-C2_2)*(H*p_i*g))))...
% == velocity_long(loop,1), H);
%
% storage_mat(loop,1) = S(2,1);
%
% elseif DEM_long(loop,1) <= threshold_2
% syms H
%
% S = vpasolve((H^(n+1)*2*A_c*(f*p_i*g*sin(slope_long(loop,1)))^n/(n+1))...
% + (H^2*sin(slope_long(loop,1))^2*((C1_3*(f*p_i*g)^2)/((1-C2_3)*(H*p_i*g))))...
% == velocity_long(loop,1), H);
%
% storage_mat(loop,1) = S(2,1);
% end
%
end
end
toc
inverted_thickness=reshape(storage_mat,size(slope));
save("Thinkness.mat","inverted_thickness")
----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
下面这是主程序(runme.m)的代码
%step 7 is specific to ISMIPA
%step 8 is specific to ISMIPF
%我选择了ISMPI的A模型，所以选择Steps = [1:7]
%steps=[1:7]; %ISMIPA
steps=[1:6,8]; %ISMIPB
%选择参数文件，IsmipA.par or IsmipF.par
%ParamFile='IsmipA.par';
ParamFile='IsmipF.par'
%运行步骤
%网格生成 #1
if any(steps==1)
%将 MD 初始化为新模型
md=model();
%生成方形网格，62.19km，82个网格点;37.5km，41个网格节点
if(ParamFile=='IsmipA.par')
md=squaremesh(md,62190,37500,82,41);
elseif(ParamFile=='IsmipF.par')
md=squaremesh(md,62190,37500,82,41);
end
% 绘制方形网格
plotmodel(md,'data','mesh');
pause(1);
% 保存给定模型
save ./Models/ISMIP.Mesh_generation md;
end
%Masks #2
if any(steps==2)
% 加载上一步的模型
md = loadmodel('./Models/ISMIP.Mesh_generation');
% 设置掩膜
% all MISMIP nodes are grounded
md=setmask(md,'','');
% 绘制掩膜 #md.mask to locate the field
plotmodel(md,'data',md.mask.ocean_levelset);
pause(1);
% 保存给定模型
save ./Models/ISMIP.SetMask md;
end
%参数化 #3
if any(steps==3)
% 加载前一步
md = loadmodel('./Models/ISMIP.SetMask');
% 模型参数化，
% 由参数文件完成参数化
md=parameterize(md,ParamFile);
% 保存模型
pause(1);
save ./Models/ISMIP.Parameterization md;
end
%挤压 #4
if any(steps==4)
%加载前一步的模型
md = loadmodel('./Models/ISMIP.Parameterization');
% 垂直拉伸前面的网格 #help extrude
% only 5 layers exponent 1
md=extrude(md,5,1);
% 绘制3D几何图形 #plotdoc
plotmodel(md,'data',md.geometry.base)
pause(1);
%保存模型
save ./Models/ISMIP.Extrusion md;
end
%设置流量计算方法 #5
if any(steps==5)
%加载前一步的模型
md = loadmodel('./Models/ISMIP.Extrusion');
% 设置流量计算的近似值
% 使用高阶模型(HO)
md=setflowequation(md,'HO','all');
% 保存模型
pause(1);
save ./Models/ISMIP.SetFlow md;
end
%设置边界条件 #6
if any(steps==6)
% 加载前一步的模型
md = loadmodel('./Models/ISMIP.SetFlow');
% 狄利克雷边界条件称为SPCs
% 冰冻结在底部，没有速度
% SPC 初始化为 NaN 每个顶点一个值
md.stressbalance.spcvx=NaN*ones(md.mesh.numberofvertices,1);
md.stressbalance.spcvy=NaN*ones(md.mesh.numberofvertices,1);
md.stressbalance.spcvz=NaN*ones(md.mesh.numberofvertices,1);
% 提取基数
basalnodes=find(md.mesh.vertexonbase);
% 将（冰）床上的滑动设置为零
md.stressbalance.spcvx(basalnodes)=0.0;
%->
md.stressbalance.spcvy(basalnodes)=0.0;
% 周期性边界必须固定在两侧
% 求域两侧的索引，先求 x，再求 y
% for x
% 创建 maxX，x 等于 x 的最大值的索引列表（使用>>帮助查找）
maxX=find(md.mesh.x==max(md.mesh.x));
% create minX，x 等于 x 的最小值的索引列表
minX=find(md.mesh.x==min(md.mesh.x));
% for y
% 创建 maxY，y 等于 y 的最大值的索引列表
% 但不是 x 等于 x 的最大值或最小值
%（即 maxX 和 minX 中的索引应从 maxY 和 minY 中排除）
maxY=find(md.mesh.y==max(md.mesh.y) & md.mesh.x~=max(md.mesh.x) & md.mesh.x~=min(md.mesh.x));
% 创建 minY，y 等于 y 的最大值的索引列表
% 但不是 x 等于 x 的最大值或最小值
minY=find(md.mesh.y==min(md.mesh.y) & md.mesh.x~=max(md.mesh.x) & md.mesh.x~=min(md.mesh.x));
% 设置应该配对在一起的节点，minX 与 maxX，minY 与 maxY
md.stressbalance.vertex_pairing=[minX,maxX;minY,maxY];
if (ParamFile=='IsmipF.par')
% 如果我们正在处理 IsmipF，则解决方案在masstransport
md.masstransport.vertex_pairing=md.stressbalance.vertex_pairing;
end
% 保存模型
pause(1);
save ./Models/ISMIP.BoundaryCondition md;
end
%Solving #7
if any(steps==7)
% 加载模型
md = loadmodel('./Models/ISMIP.BoundaryCondition');
% 设置集群 #md.cluster
% 通用参数 #help generic
% 仅设置进程的名称和数量
md.cluster=generic('name',oshostname(),'np',2);
% 设置要查看的信息 #help verbose
md.verbose=verbose('convergence',true);
% Solve, we are solving a StressBalanc
md=solve(md,'Stressbalance');
% 保存给定模型
save ./Models/ISMIP.StressBalance md;
% 绘制表面速度场
plotmodel(md,'data',md.results.StressbalanceSolution.Vel)
pause(1);
end
%Solving #8
if any(steps==8)
% 加载前一步的模型
md = loadmodel('./Models/ISMIP.BoundaryCondition');
% 设置集群 #md.cluster
% 通用参数 #help generic
% 仅设置进程的名称和数量
md.cluster=generic('name',oshostname(),'np',2);
% 设置要查看的信息 #help verbose
md.verbose=verbose('convergence',true);
% 将瞬态模型设置为忽略热模型
md.transient.isthermal=0;
% 定义时间步长方案
% 这里的一切都应该以年为单位提供。
% 设置时间步长
md.timestepping.time_step=1;
% 设置结束时间 (20*4 years time_steps)
md.timestepping.final_time=1*33;
% Solve #help solve
% we are solving a TransientSolution
md=solve(md,'Transient');
% 保存模型
save ./Models/ISMIP.Transient md;
% 绘制表面速度场
plotmodel(md,'data',md.results.TransientSolution(33).Vel);
pause(1);
end
----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
下面这是主程序的参数文件(ismipF.par)的代码——
%ISMIP F实验的参数化
%设置模拟通用名称#md.miscellaneous
disp('Constructing Geometry');
% 定义上表面的几何模型的形状
filename = 'Dem_90m_New_reverse.tif';
img = imread(filename);
img_Array = double(img);
img_0_1 = img_Array(:);
n = length(img_0_1);
B = zeros(size(img_0_1));
md.geometry.surface=img_0_1;
% L是这个方形的大小
L=max(md.mesh.x)-min(md.mesh.x);
filename_1 = 'Dem_Mask_1_0_New_reverse.tif';
img_1 = imread(filename_1);
img_Array_1 = double(img_1);
img_1_1 = img_Array_1(:);
%下表面与厚度的几何模型的形状
img_1_1 = 60*img_1_1 + 0.0000001
md.geometry.thickness=img_1_1;
md.geometry.base=md.geometry.surface-md.geometry.thickness;
%绘制几何图形进行查验
plotmodel(md,'data',md.geometry.thickness);
pause(1);
disp('Defining friction parameters');
%这些参数不会被使用，但需要修复 #md.friction
%每个节点一个弗里西顿系数 (md.mesh.numberofvertices,1)
%将年转换为秒 #md.constants.yts
md.friction.coefficient=sqrt(md.constants.yts/(1000*2.140373*10^-7))*ones(md.mesh.numberofvertices,1);
%每个元素一个弗里西顿指数（p，q）
md.friction.p=ones(md.mesh.numberofelements,1);
md.friction.q=zeros(md.mesh.numberofelements,1);
disp('Construct ice rheological properties');
%流变参数位于材料部分 #md.materials
%B 每个顶点有一个值
md.materials.rheology_B=(1/(2.140373*10^-7/md.constants.yts))*ones(md.mesh.numberofvertices,1);
%n 每个元素有一个值
md.materials.rheology_n=1*ones(md.mesh.numberofelements,1);
disp('Set boundary conditions');
%设置冰盖的默认边界条件
md=SetIceSheetBC(md);
disp('Initializing velocity and pressure');
%初始化速度场和压力场 #md.initialization
md.initialization.vx=zeros(md.mesh.numberofvertices,1);
md.initialization.vy=zeros(md.mesh.numberofvertices,1);
md.initialization.vz=zeros(md.mesh.numberofvertices,1);
md.initialization.pressure=zeros(md.mesh.numberofvertices,1);
2397
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
