Hi Johanna,
This problem is similar to the one you mentioned earlier about the calculation of flux and is caused by the use of continuous grounding line migration in the elements, which causes some elements to be partially floating and partially grounded at the same time. So when you do min(ocean_levelset(index),[],2)<0
for example, it does not know that some elements are only partially floating and considers them entirely floating, leading to some discrepancies. The difference is much smaller for the volume above floatation for example because the ice is floating or close to floating at the grounding line and therefore the error is small compared to the volume above floatation over the thicker part of the ice sheet.
Here is some code I used to compare with a similar calculation in Matlab and it gives the same results. I used test 504 adding grounding line migration with md.transient.isgroundingline=1;
and md.groundingline.migration='SubelementMigration';
and refining the mesh.
for step = 1:numel(md.results.TransientSolution)
%Get Grounding line segments for this step
gl1= isoline(md, md.results.TransientSolution(step).MaskOceanLevelset, 'value', 0, 'output','matrix');
%Get speed and thickness for this step at the segment extremities
vx = InterpFromMesh2d(md.mesh.elements,md.mesh.x,md.mesh.y,md.results.TransientSolution(step).Vx,gl1(:,1),gl1(:,2));
vy = InterpFromMesh2d(md.mesh.elements,md.mesh.x,md.mesh.y,md.results.TransientSolution(step).Vy,gl1(:,1),gl1(:,2));
h = InterpFromMesh2d(md.mesh.elements,md.mesh.x,md.mesh.y,md.results.TransientSolution(step).Thickness,gl1(:,1),gl1(:,2));
areas = GetAreas(md.mesh.elements,md.mesh.x,md.mesh.y);
ocean_mask=md.results.TransientSolution(step).MaskOceanLevelset;
floating_area = 0;
grounded_area = 0;
for elem=1:md.mesh.numberofelements,
mask = ocean_mask(md.mesh.elements(elem,: ));
if mask(1)>0 & mask(2)>0 & mask(3)>0, % all grounded
grounded_area = grounded_area + areas(elem);
elseif mask(1)<0 & mask(2)<0 & mask(3)<0, % all floating
floating_area = floating_area + areas(elem);
else
if mask(1)*mask(2)*mask(3)>0,
mainlyfloating = 0;
else
mainlyfloating = 1;
end
if mask(1)*mask(2)>0, % nodes 1 and 2 are similar
s1 = mask(3)/(mask(3)-mask(2));
s2 = mask(3)/(mask(3)-mask(1));
elseif mask(2)*mask(3)>0, % nodes 2 and 3 are similar
s1 = mask(1)/(mask(1)-mask(2));
s2 = mask(1)/(mask(1)-mask(3));
elseif mask(1)*mask(3)>0, % nodes 1 and 3 are similar
s1 = mask(2)/(mask(2)-mask(1));
s2 = mask(2)/(mask(2)-mask(3));
else
error('not possible');
end
if mainlyfloating==1,
phi = (1-s1*s2);
else
phi = s1*s2;
end
grounded_area = grounded_area + phi*areas(elem);
floating_area = floating_area + (1-phi)*areas(elem);
end
end
%Convert form volume to mass in Gt/yr
disp(['Manually computed floating area: ' num2str(floating_area) ]);
disp(['ISSM computed floating area: ' num2str(md.results.TransientSolution(step).FloatingArea) ]);
disp(['Manually computed grounded area: ' num2str(grounded_area) ]);
disp(['ISSM computed grounded area: ' num2str(md.results.TransientSolution(step).GroundedArea) ]);
end
which gives these results:
Manually computed floating area: 1282394022.4674
ISSM computed floating area: 1282394022.4675
Manually computed grounded area: 18538348745.8331
ISSM computed grounded area: 18538348745.833
Manually computed floating area: 1512344737.3878
ISSM computed floating area: 1512344737.3879
Manually computed grounded area: 18308398030.9127
ISSM computed grounded area: 18308398030.9127
Manually computed floating area: 1618345039.4635
ISSM computed floating area: 1618345039.4635
Manually computed grounded area: 18202397728.837
ISSM computed grounded area: 18202397728.837
Manually computed floating area: 1683422013.4702
ISSM computed floating area: 1683422013.4701
Manually computed grounded area: 18137320754.8303
ISSM computed grounded area: 18137320754.8304
Regarding the scaled version of the scalars, this allows to include the area distortion that we introduced by projecting the curved earth surface onto a flat surface and therefore be closer to the exact area. It is available for a number of scalar outputs.
Best
Helene