Incorporating PICOP in transient model

AHill1010

Hello!

I am attempting to incorporate PICOP, a warm cavity plume model, into my transient scenario. The information on PICOP is in the 2019 User Guide, section 5.4, beginning on page 147. My set up of the transient model is based off of the Pine Island Stability model in the tutorial.

`if any(steps==4) %PICOP Model

%load model from other folder and call it md4
md4 = loadmodel('../ThwaitesFull2/Models/Thwaites_Control_friction');

%Set controls on transient solution
md4.groundingline.migration='SubelementMigration';
md4.inversion.iscontrol=0;
md4.transient.ismasstransport=1;
md4.transient.isstressbalance=1;
md4.transient.isgroundingline=1;
md4.transient.ismovingfront=0;
md4.transient.isthermal=1;

%Define simulation length and timestep
md4.timestepping.time_step=1;
md4.timestepping.final_time=10;
md4.transient.requested_outputs={'default','IceVolume','IceVolumeAboveFloatation'};

%Basal Forcings
md4.basalforcings = basalforcingspico();
md4.basalforcings.basin_id = ones(md4.mesh.numberofelements,1);
md4.basalforcings.num_basins = 1;
md4.basalforcings.maxboxcount = 5;
md4.basalforcings.groundedice_melting_rate = zeros(md4.mesh.numberofvertices,1);
%Here I define the farocean temperature and salinity in a matrix, as
%is done in other steps of the transient modle.
md4.basalforcings.farocean_temperature = 272.*ones(md4.mesh.numberofvertices,1);
md4.basalforcings.farocean_salinity = 34.73.*ones(md4.mesh.numberofvertices,1);
%Here I just assign the temperature and salinity a value. Regardless of
%which option I choose, the same error results.
%     md4.basalforcings.farocean_temperature = 272;
%     md4.basalforcings.farocean_salinity = 34.73;
   
%Turn on/off the PICOP model, if off, =0, then just the PICO model is
%run. I figured I would try to get that running before adding in the
%plume part.
md4.basalforcings.isplume = 0;

%Redefine parameters.
disp('   Loading SeaRISE data from NetCDF');
ncdata='../Data/Antarctica_5km_withshelves_v0.75.nc';
gflux = ncread(ncdata,'bheatflx_fox')';
x1 = ncread(ncdata,'x1');
y1 = ncread(ncdata,'y1');
disp('   Set geothermal heat flux');
md4.basalforcings.geothermalflux=InterpFromGridToMesh(x1,y1,gflux,md4.mesh.x,md4.mesh.y,0);
clear gflux;

%Solve
md4 = solve(md4,'Transient');

%Save
save ./Models/Thwaites_TEST md4;

end`

The resulting errors say

checking model consistency
model not consistent: field 'basalforcings.farocean_temperature' dimension # 1 should be of size 2
model not consistent: field 'basalforcings.farocean_salinity' dimension # 1 should be of size 2
Error using ismodelselfconsistent (line 36)
Model not consistent, see messages above
Error in solve (line 90)
ismodelselfconsistent(md),
Error in ThwaitesSensitivity2 (line 137)
md4 = solve(md4,'Transient');

I define both far ocean temperature and salinity, and I don't know why the first dimension needs to be 2. I tried messing with the dimensions but continued to get the same error. Any suggestions?

Thanks!
Ann!

mathieumorlighem

Hi Ann,

you should have one temperature and one salinity per basin and the second raw is to define the time. In your case, you only have 1 basin, so your far field forcing should look like:

md.basalforcings.farocean_temperature = [271.15 272.15 273.15; 0.5 1 1.5]; %K
md.basalforcings.farocean_salinity    = [31 32 33; 0.5 1 1.5]; %PSU

the second row is the time (in years). A linear interpolation between the different time is applied. I hope this helps!
Mathieu

AHill1010

mathieumorlighem Thank you so much for your response!

I added your two lines of code, both after my two lines where I define the matrices of far ocean temperature and salinity, as well as substituting for those two lines. Unfortunately, I had no success. When replacing the lines of code I got the additional errors shown below, and when adding your lines of code I got the same errors as I was previously getting. Since the error says the first dimension must be 2, I changed both matrices to 2x2, but that didn't work. Since you said I should have one temperature and salinity per basin, I changed both matrices to 1x1. This was also unsuccessful. Is there something that I am missing?

ThwaitesSensitivity2
Loading SeaRISE data from NetCDF
Set geothermal heat flux
checking model consistency
model not consistent: field 'basalforcings.farocean_temperature' dimension # 1 should be of size 2
model not consistent: NaN values found in field 'basalforcings.farocean_temperature'
model not consistent: field 'basalforcings.farocean_salinity' dimension # 1 should be of size 2
model not consistent: NaN values found in field 'basalforcings.farocean_salinity'
Error using ismodelselfconsistent (line 36)
Model not consistent, see messages above
Error in solve (line 90)
ismodelselfconsistent(md),
Error in ThwaitesSensitivity2 (line 139)
md4 = solve(md4,'Transient');

Thank you so much for your help!

Ann

mathieumorlighem

Dear Ann,

it is working on my end so I don't think the code you added is correct 😀 Maybe you forgot to change md to md4? make sure md.basalforcings.farocean_temperature looks right before md=solve(md,'tr');
Mathieu

AHill1010

@mathieumorlighem That was my problem! I simply didn't change md to md4. Thank you so much for all of your help, I am extremely appreciative! If you don't mind, I would like to make sure I understand what you did. You specified the temperature at three different times, allowing the model to fill in the temperature for all times. This is how you make temperature variable throughout the model run. I also tried just giving one temperature and one time, which also worked, so I assume that defines temperature as a constant. Thanks again!
Best,
Ann

mathieumorlighem

Hi Ann, yes you are absolutely right. Each column of the far field properties correspond to the time given in the last row (in years), and if a time step falls in between two of these times, ISSM uses a linear interpolation. Don't hesitate if you have any other questions.
Mathieu

AHill1010

@mathieumorlighem Ok, that sounds great! Thank you again for all of your help, I greatly appreciate it!
Best,
Ann