Rework glacier regions to prevent conversion of negative ice runoff to liquid and sensible heat flux

[billsacks]

Is it a problem that conversion of snow capping to liquid also operates on negative snow capping fluxes?

@lvankampenhout raised this question by email; I am moving this discussion to here to get broader feedback. For now I'm bringing @dlawrenncar @olyson @whlipscomb into the loop; we may want to bring others into the loop, too.

It's important to note that Leo's analysis was done on a run where he had accidentally set up the glacier regions so that all of Greenland had glacier_region_ice_runoff_behavior = 'melted'. In out-of-the-box cases, Greenland itself has glacier_region_ice_runoff_behavior = 'remains_ice', and I think the problem he's seeing only would occur in the region of the world that is part of the CISM domain but outside Greenland itself (this includes the Canadian archipelago and some other small land areas).

Here's Leo's email:

I just took a second look at FSH_RUNOFF_ICE_TO_LIQ because I found it remarkable that NCView indicated that there were positive values. Yet from the description of the variable (and the comments in the code) you’d expect only negative values. The positive values occur over large extents from April to September, with peaks over 150 W/m2 in summer, which are quite substantial. Here is July: https://www.dropbox.com/s/8hxwcvx9ewez74m/fig1.png?dl=0

First I figured this could arise when there is irrigation over the tundra or something , like I saw in one of my earlier test simulations with modified soil parameters, leading to a negative runoff. However, this explanation is at least not the full story, since the positive SH values also occur over 100% glacier grid cells, as seen here (white line) https://www.dropbox.com/s/vr2mwmkgai87rnk/fig2.png?dl=0

Looking at the code I believe that this bit may be responsible for the behaviour that I’m seeing (from GlacierSurfaceMassBalanceMod.F90):

      ! Yet another way to think about this is: When ice melted, we let the liquid run
      ! off, and replaced it with new ice from below. But that new ice needed to come
      ! from somewhere to keep the system in water balance. We "request" the new ice from
      ! the ocean by generating a negataive ice runoff equivalent to the amount we have
      ! melted (i.e., equivalent to the amount of new ice that we created from below).

      ! As above: Melt is only generated for glacier columns. But it doesn't hurt to do
      ! this for all columns in the do_smb filter: this melt term will be 0 for other
      ! columns.

      qflx_ice_runoff_snwcp(c) = qflx_ice_runoff_snwcp(c) - (1.0_r8 - glc_dyn_runoff_routing(g)) * qflx_glcice_melt(c)

The snow capping term that is modified here serves as the basis for the solid runoff term

lnd2atmMod.F90

         qflx_ice_runoff_col(c) = waterflux_inst%qflx_ice_runoff_snwcp_col(c) + &
              waterflux_inst%qflx_ice_runoff_xs_col(c)

So the line of events would be that by taking water from the ocean to replace melting ice, we are generating a negative runoff term that then triggers the code that generates a sensible heat flux. If this is the case, then there is an artificial positive feedback loop where melting of ice generates sensible heat, thereby speeding up the melt . Are you aware of this issue, i.e. is this the trade-off for the Labrador Sea problem, or is this a bug?

I agree with Leo's sense that this is probably due to the negative ice runoff term that can be generated in some cases. This negative ice runoff term can appear over glacier columns that have glacier_region_melt_behavior = 'replaced_by_ice', and the SH flux in question can appear over glacier columns that have glacier_region_ice_runoff_behavior = 'melted'. Out-of-the-box, the only area that has both of these behaviors is the relatively small land area that is within the CISM domain but outside of Greenland itself (including the Canadian archipelago and some other small areas). So in a global sense, this may not be hugely important out-of-the-box – although the areas affected may be important for sea ice, and also note that users are free to set up these glacier region behaviors however they want (so some users' runs may have larger affected regions).

It's hard to know what the "right" thing is to do here, since this arises from the interaction of two non-physical terms. On an annual-average basis, for a system near steady state, I can convince myself that the current behavior (where negative ice runoff is converted to negative liquid runoff and a positive sensible heat flux) is "correct" in some sense: Just as we need the negative ice runoff – in order to compensate for too-large positive ice runoff at other times and places, and in order to conserve energy, as noted in https://github.com/ESCOMP/ctsm/blob/fd7a88a1a052014c5c9efafe7aa12ead3c2dfa5e/src/biogeophys/GlacierSurfaceMassBalanceMod.F90#L397-L425 – so, too, this positive SH flux can be seen as compensating for too large of a negative SH flux at other times and places. However, it may be that this positive SH flux causes problems in practice, such as the positive feedback that Leo suggests.

From a conservation standpoint, it would be just as fine to keep this negative ice runoff as ice, thus adding heat to the ocean rather than to the atmosphere. There would be some asymmetry here: positive ice runoff (from snow capping) extracts heat from the atmosphere, whereas negative ice runoff (generated during glacial melt) adds heat to the ocean. But, given that these are all somewhat non-physical fluxes anyway, I'm not sure that that asymmetry is inherently wrong. Implementing this asymmetry – if desired – should be fairly easy: we can either put the negative ice runoff in a term other than qflx_ice_runoff_snwcp in https://github.com/ESCOMP/ctsm/blob/fd7a88a1a052014c5c9efafe7aa12ead3c2dfa5e/src/biogeophys/GlacierSurfaceMassBalanceMod.F90#L431 or we can add a conditional in handle_ice_runoff so that whether the conversion is done depends on the sign of qflx_ice_runoff_col.

@lvankampenhout @dlawrenncar @olyson @whlipscomb - do any of you have feelings on what is the "right" thing to do here? Also, should we bring anyone else into this discussion at this point?

[billsacks]

I talked briefly with @olyson about this. He didn't have any immediate thoughts, but suggested:

1. We might want to do some simulations looking at the impact of treating this in different possible ways, as an aid to making this decision.

2. We might want to consider what would happen under 2 or 4xCO2 - where glacial melt would be greater.

3. We might want to look more closely at how big of an impact this is having in the current CESM2 runs (or maybe that was just my thought).

(There may have been one other suggestion, too – @olyson , please chime in if you remember any of your other thoughts.)

[billsacks]

I have talked with @dlawrenncar @olyson @whlipscomb and others about this over the last few days. Here is a summary of these discussions:

Recall that, with snow capping, in regions where ice runoff is converted to liquid, we have:

snow capping -> ice runoff -> liquid runoff & -SH

It turns out that, in regions where this behavior is combined with the behavior of having glacial melt runoff, we also have:

ice melt -> +liquid runoff & -ice runoff Then: -ice runoff -> -liquid runoff & +SH

On net, then, ice melt -> +SH

Dave Lawrence points out that the final, net result is that we're just not letting the glacial column melt, and are taking the energy that is trying to melt it and sending it back to the atmosphere as a sensible heat flux. This is similar to treating glacier as soil for the sake of these energy fluxes, except that all of the energy goes back as SH, rather than LW radiation, etc.

Keith Oleson showed plots of July SH showing that this is mainly a problem in very warm runs, e.g., 4xCO2. In the 4xCO2 run, the SH from this conversion flux can be a few times as large as the background SH flux.

Dave, Keith and I discussed the possibility of doing something else with this negative ice runoff; for example:

• We could keep it as negative ice, but Dave and/or Keith worried that might have bad implications for ocean circulation
• We could put it in longwave, but that might similarly heat the atmosphere
• We could put it in reflected solar, but that could lead to > 1 albedo

We decided that the best solution is to rework the glacier regions so there's no instance of this combination of behaviors (glacier_region_melt_behavior = replaced_by_ice, glacier_region_ice_runoff_behavior = melted). Specifically, we will remove the small region that is currently defined as being within the CISM domain but outside Greenland itself, treating that region as we do the rest of the world outside of Greenland and Antarctica. This means that we won't be computing SMB over quite the whole CISM domain – instead, we'll always send 0 SMB from Ellesmere island and other places that are within the CISM domain but outside of Greenland.

Members of the LIWG are happy with this solution.

The current Greenland GLACIER_REGION on the surface datasets seems as good as anything for this purpose: it was created to be consistent with CISM's default 4-km input file, treating everything that CISM considers land or ice-covered to be within this region.

I might do this in multiple stages:

1. Keep the current GLACIER_REGION, but change the behavior of the in-CISM-outside-Greenland region to be virtual, remains_in_place, melted. The reason for keeping this virtual for now is (a) to avoid needing to interpolate initial conditions, and (b) to avoid even more answer changes. This would also avoid the need for recreating all surface datasets now. However, I may want to put in place all of the necessary code mods while they're fresh in my head, and do some testing with (2).

2. Change this region's glacier_region_behavior to "single" – so it's now single, remains_in_place, melted, just like most of the world.

3. Remake surface datasets with 3 rather than 4 glacier regions. This should be bit-for-bit with (2).

Note that, once we get to stage (2), there would be a small part of the CISM domain that no longer has virtual columns. We'll need to be careful to ensure that we don't try to change glacier cover within that region - i.e., we'll need to make sure that we only update glacier cover in places that are both (a) within CISM's "icemask", and (b) have "virtual" behavior. (Previously we did this by ensuring that all grid cells within the icemask had virtual behavior; this will no longer be possible so we'll need to update the logic to handle this.)

Similarly, in glc2lndMod, we'll need to change update_glc2lnd_dyn_runoff_routing to not just set glc_dyn_runoff_routing_grc = icemask_coupled_fluxes_grc – but instead to take into account glacier behavior: wherever we're not computing SMB, we should set glc_dyn_runoff_routing_grc = 0.

[whlipscomb]

@billsacks, Thanks for the summary. Your multiple-stages plan makes sense to me. Please let me know if any CISM changes would be helpful, in either the short or longer term.