Antarctica ice shelves are being treated as wetlands rather than glaciers

[billsacks]

Brief summary of bug

In recent surface datasets, including the ones being used for CMIP6 runs, Antarctica's ice shelves are being treated as wetlands rather than glaciers.

General bug information

CTSM version you are using: As far as I can tell, affects all versions for the last year or two.

Does this bug cause significantly incorrect results in the model's science? No

Although the landunit type is incorrect over the ice shelves, at least in the one run I checked (b.e21.BHIST.f09_g17.CMIP6-historical.001), the affected regions are still covered in deep snow year-round, so the impact on land-atmosphere fluxes should be minimal.

Configurations affected: All (I think)

Details of bug

General description

Jan Lenaerts noticed that there are no glaciers over Antarctica's ice shelves in the CMIP6 runs (for example, 80 S, 180 E/W (Ross Ice Shelf)). This is very apparent if you look at the 4th element (glc_mec) of PCT_LANDUNIT in a history file from one of these runs.

I think I have traced the problem back to a change in the LANDMASK field on the new vegetation raw data files that we've been using for the last year or two: This field is 0 over the ice shelves in /glade/p/cesm/cseg/inputdata/lnd/clm2/rawdata/pftcftdynharv.0.25x0.25.LUH2.histsimyr1850-2015.c170629/mksrf_landuse_histclm50_LUH2_1850.c170629.nc, whereas it is 1 over ice shelves in the files we used before a couple of years ago (/glade/p/cesmdata/cseg/inputdata/lnd/clm2/rawdata/pftlandusedyn.0.5x0.5.simyr1850-2005.c090630/mksrf_landuse_rc1850_c090630.nc).

Probable cause

For better or worse (probably worse), mksurfdata_map uses the landmask on the PFT rawdata file to determine where to set landcover to 100% wetland:

https://github.com/ESCOMP/ctsm/blob/a34404419aaa81f5021602963222c50daa10c0f6/tools/mksurfdata_map/src/mksurfdat.F90#L707-L724

My understanding of the code is that the LANDMASK field isn't used directly in the mksurfdata_map code, but (if I remember correctly) I think it's used indirectly: I think the LANDMASK field on that file is used to set the mask field on the mapping files, and the mask field on the mapping files determines the mask in mksurfdata_map.

What this does and doesn't impact in CMIP6 runs

Fortunately, since we're using virtual glacier columns over Antarctica in the non-WACCM CMIP6 runs, I think this bug doesn't cause any loss of information from the land ice perspective – but some diagnostics may require more post-processing than they otherwise would have.

In principle, using the wrong landunit here could lead to scientifically incorrect results. However, at least in the one run I checked (b.e21.BHIST.f09_g17.CMIP6-historical.001), the affected regions are still covered in deep snow year-round, so the impact on land-atmosphere fluxes should be minimal.

Suggested fix(es)

I can see two ways to fix this issue:

(1) Fix the LANDMASK field on the raw vegetation files. This would require putting valid data on that file for the Antarctic ice shelves (this could be called bare ground everywhere in that region).

(2) Remove the code referenced above which uses the PFT data's LANDMASK to 0 out everything. I really don't like that code, and would like us to remove it, but some analysis is required to determine if its removal would have unforeseen consequences, and whether it needs to be replaced by other, similar code.

@whlipscomb @gunterl @Katetc @dlawrenncar @lawrencepj1 @ekluzek @olyson

[billsacks]

I'm adding the "impacts science" label because, in discussions with @whlipscomb , I realized that this could still impact science in scenario runs where the snowpack melts.

[billsacks]

Based on discussion with @dlawrenncar and @lawrencepj1 - in the short-term, we will fix this problem by expanding the land mask on the pft raw data files so that it includes areas classified as glacier in the glacier raw data set, for Antarctica. Something like:

if (latitude(i) <= -60 .and. pct_glacier(i) > 0) then
   landmask(i) = 1
end if

@lawrencepj1 - here's the glacier raw data file you can use for this: /glade/p/cesmdata/cseg/inputdata/lnd/clm2/rawdata/mksrf_glacier_3x3min_simyr2000.c120926.nc

In the longer term we should do something like #553 , which moves the above logic (and more) into mksurfdata_map, but that carries a bit more risk.

[whlipscomb]

@billsacks, This looks good to me for the short term. To make sure we've done our due diligence: Has someone checked that the wetland conversion wasn't erroneously applied to glaciers elsewhere, such as Greenland or Ellesmere?

[dlawrenncar]

Related to Bill Lipscomb's point, why are we including the restriction to only make this correction where latitude <=-60 . Is that just to be safe so as not to change anything anywhere else? But, then maybe if there is a misclassified land point (wetland where it should be glacier) in the northern high latitudes, we should change that also?

…

On Tue, Oct 30, 2018 at 12:02 PM William Lipscomb ***@***.***> wrote: @billsacks <https://github.com/billsacks>, This looks good to me for the short term. To make sure we've done our due diligence: Has someone checked that the wetland conversion wasn't erroneously applied to glaciers elsewhere, such as Greenland or Ellesmere? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#545 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AUAcVO9jY72DZWOq0LsVG-8w5tXdlpaJks5uqJRJgaJpZM4X6iLy> .

[billsacks]

@whlipscomb and I talked with @lawrencepj1 over lunch. We decided that we'll at least try doing this correction globally. My concern was exactly what @dlawrenncar raised, and particularly that it might change land-surface fluxes if it's applied in a place that is not snow-covered in the historical run. But I'd say: let's apply it globally and then look carefully at where there are differences.

[billsacks]

See this comment: #553 (comment)

[billsacks]

From discussion with @lawrencepj1 - he's going to try a different, more conservative approach than the one laid out in #553 (comment) - within the block that typically sets wetlands everywhere, he'll check if pct_glacier > 0, and if so, set the grid cell coverage to glaciers rather than wetlands. This doesn't come as close to the long-term solution he wants, but this way he won't have to deal with what to do about pct_pft, soil color, etc., and we agreed that it's better to take this safer approach for the cmip6 runs.

[billsacks]

The necessary mksurfdata_map changes are implemented in #564 . We also will need updated surface datasets generated with the new mksurfdata_map. As noted in #564 (comment), we want these both on master and on the clm5 release branch. @ekluzek has agreed to do that, so I'm assigning this issue to him for this remaining piece.

[ekluzek]

Fix in release-clm5.0.16