Consider one or more changes to reduce DYNBAL fluxes (dynamic landunit conservation correction fluxes)

[billsacks]

Introduction

When subgrid areas change due to dynamic landunits / columns, fluxes are generated to compensate for changes in gridcell-average water and energy (QFLX_LIQ_DYNBAL, QFLX_ICE_DYNBAL, EFLX_DYNBAL). These fluxes can (and should) be reduced via careful setting of reset_dynbal_baselines, but that still doesn't reduce these fluxes to zero.

The current implementation seems good enough in situations where there is a gradual change in a given grid cell's subgrid areas over many years. Unfortunately, though, that is not always the situation when coupling to CISM, as has been seen in recent CESM3 runs: In grid cells around the ice sheet margin, large fractions of a CTSM grid cell can blink between mostly-vegetated/bare and mostly-glacier on a 2-year cycle, and this oscillatory behavior can go on for many years, possibly indefinitely. This can lead to large dynbal fluxes in each year, even when dribbled throughout the year (as is currently done).

It is problematic for the ocean to get negative dynbal fluxes, but it is also problematic to remove these negative dynbal fluxes by spreading them globally, as is currently done, because then regional fluxes that should average to 0 spatially and temporally (e.g., into the Arctic ocean) end up being biased positive (as explained here). Therefore, I have been thinking about possible strategies to further reduce these DYNBAL fluxes. I'm laying out a few possibilities here, and we may want to implement one or more of these.

(1) Add pools that slowly release the dynbal fluxes, with multi-year residence times

Currently, the dynbal fluxes - which are created on year boundaries - are dribbled out over the next year using infrastructure that ensures conservation on an annual basis. This is satisfying in the sense that water and energy are not allowed to get out of balance from one year to the next. But this can lead to larger-than-desirable fluxes, especially in the case of oscillatory behavior as described above.

To further reduce these fluxes – especially in the oscillatory situation where the year y+1 fluxes are approximately equal and opposite of the year y fluxes in a given grid cell – we could instead introduce holding pools for these DYNBAL fluxes, where the liquid water, ice, and energy are released slowly, with a multi-year residence time. The residence time should be chosen to be short enough to not have a significant impact on system-wide water and energy balance, but long enough that the fluxes are small. Note that (I think...), for oscillatory fluxes, the outfluxes will be biased at first (positive if the first-year flux is positive, negative if the first-year flux is negative), but will eventually settle to a multi-year average of 0 (in the asymptote, anyway); longer residence times will mean it takes longer to settle to this near-zero multi-year average. I would somewhat arbitrarily suggest a 10-year residence time; for oscillatory fluxes, this would reduce fluxes by a factor of 10 initially, and a factor of 20 in the long-term behavior after a few decades of spinup (at which point out-fluxes from the holding pool would alternate between roughly 0.05 and -0.05 times the magnitude of the annual oscillatory flux... I think). (Update 2026-08-28: I think I'd actually be inclined to a slightly longer residence time, like 20 years, to reduce these fluxes even further.)

I think this would be relatively easy to implement. I think it would also be relatively easy to test, at least for a basic sanity test that gives a fair amount of confidence: After the first change in glacier areas (i.e., at a year boundary), compare the fluxes from these pools with the fluxes from a baseline run. In the first few days, the fluxes from these pools should be almost exactly equal to the baseline fluxes divided by the residence time (in years), but these new fluxes should gradually decline in magnitude compared with those from the baseline.

(2) Include H2OSNO in the dynbal baselines

Currently, the dynbal baselines only account for below-ground water states. This is physically justifiable because these below-ground states are virtual states for glacier (the glacial ice, which has an assumed depth), and glacier columns don't have an actual soil representation. The remaining dynbal fluxes are then largely driven by differences in H2OSNO between vegetated and glacier columns.

I think we could include H2OSNO (and possibly other, small water/ice states) in the dynbal baselines. Then, as long as reset_dynbal_baselines is done after an initial spin-up, I think this would reduce the dynbal ice flux.

I think this would be pretty easy to implement. But I'd want to give this more thought - and probably go back through notes from when this was first developed - to feel comfortable that this is an okay thing to do. One possible concern is what would happen if the snow state has evolved significantly since the dynbal baselines were reset: could this actually lead to larger dynbal fluxes than before? (I'm not sure that it would, but I'd want to think this through.) In addition, this inclusion of H2OSNO in the baselines feels harder to justify and explain, even if it arguable still conserves mass and energy (because I think we're free to set the baselines however we want for conservation purposes, as long as we keep them fixed once the "real" coupled model simulation begins).

(3) Completely zero out dynbal fluxes in the early phases of coupled model spinup

The dynbal fluxes are needed for conservation in a coupled simulation, but it seems arguable whether conservation is actually important in the early stages of coupled model spinup. In particular, note that we do some initial coupled model spinup followed by additional offline spinup of the land and ocean components, which breaks conservation, before continuing the coupled model spinup. In that first phase of coupled model spinup, we could arguably zero out the dynbal fluxes. This feels equivalent to acknowledging that the initial ice sheet state is wrong / out of balance, and we let it come into initial balance without worrying about conservation.

Some issues with this, though, are:

• This would probably need to be done manually - setting this zeroing flag to true and false at the correct times - and this is error prone.
• Whenever this flag is changed, model behavior will change, potentially shocking the system. So, for example, it could be problematic to do this right before switching from the spinup to the transient run, because then the behavior in the transient run will differ from the spinup.
• While this could help with problems of the initial system being out of balance (giving it some time to come into balance), it wouldn't help with the oscillatory situation described above, where no amount of spinup will be enough to stop these fluxes.

[billsacks]

My feeling is that (1) is probably the most worthwhile in terms of a combination of ease of implementation, effectiveness, and justifiability.

[adamrher]

Thanks for explaining all these nuances arising from dynamic landunits and an evolving (expanding/retreating) ice sheet. I am struck by your description of this numerical instability:

Around the margins of the ice sheet, there is an oscillatory behavior where the ice sheet blinks in and out of existence on a 2-year cycle. I think this arises because: In years where there is no ice sheet in a given land grid cell, there is no capacity for the land model to generate negative surface mass balance. The ice sheet can then freely advance into that grid cell without being melted. Once it advances, however, in the next year the land model does generate a net-negative surface mass balance. This negative surface mass balance then melts the ice in that grid cell to a greater extent than the inflow of ice, resulting in zero ice in the following year.

I feel like there's enough information in a non-glaciated grid cell to provide a first guess of the surface mass balance of an advancing glacier front into that grid cell. For example, what if we had a default "advancing ice reference profile" used to compute a surface energy balance and SMB at non-glaciated grid points, at each EC, which would get activated once the area fractions of any of the ECs exceed 0? Apologies if this is ridiculous, Im just riffing.

[billsacks]

I feel like there's enough information in a non-glaciated grid cell to provide a first guess of the surface mass balance of an advancing glacier front into that grid cell. For example, what if we had a default "advancing ice reference profile" used to compute a surface energy balance and SMB at non-glaciated grid points, at each EC, which would get activated once the area fractions of any of the ECs exceed 0? Apologies if this is ridiculous, Im just riffing.

@adamrher - thanks a lot for sharing your thoughts here. I was actually thinking the same thing a few days ago, wondering why we didn't do something like what you suggested. (I didn't actually get as far into it as you seem to have gotten with your idea of an advancing ice reference profile... I was just thinking about this more vaguely....) I can't remember for sure, but I think at least part of the issue is how to ensure conservation: With the current coupling, it's safe for CTSM to send any amount of positive SMB to CISM, but I don't think it's safe to send arbitrary negative SMB: any negative SMB has already created a runoff flux and so needs to be satisfied by CISM in order to conserve.

As I thought about that, I realized that we might already have issues with non-conservation in the current implementation, if the negative SMB in a year exceeds the amount of ice in the grid cell. I have a vague recollection of conversations where we recognized this as something we should resolve at some point, but as far as I can remember, I think this is still an issue. But sending a pre-emptive negative SMB everywhere would (I think) cause a much bigger potential conservation issue in this respect.

That said, this is probably a solvable problem... it would just take some thought and work. For example, we could change the coupling so that negative SMB does not actually generate runoff from CTSM - instead it is treated more like positive SMB, in that it's just sent to CISM to deal with. Then CISM would apply as much of the negative SMB as it can, and CISM would be responsible for generating a liquid runoff flux corresponding to the actual amount of melted ice. (CISM already generates a small liquid runoff flux for the basal melt, so this could piggy-back on that existing coupling field.)

I do think that something like this could be worth doing to prevent these oscillations, at least if the CISM work doesn't feel too hard. But this feels like a much bigger task than any of the CTSM-specific ideas I laid out above, and so would not be feasible for CESM3.

@whlipscomb - I'd be interested in your thoughts on this... though maybe we should move this piece to a CISM issue for further discussion.

[billsacks]

Here are some other ideas:

(4) Do some regional spatial canceling / offsetting of negative and positive values

Rather than relying on CMEPS to cancel negative and positive values – which does a global spreading of negative values – we could instead do something similar within CTSM, but separated by region. That is, we'd sum the total DYNBAL fluxes within a given region; if the sum is positive, we'd zero all negative values and downweight the positive values proportionally; and vice versa if the sum is negative. A logical way to do this could be based on the existing glacier regions - so, for example, this would be done regionally over Greenland. We could potentially also have a new glacier region behavior flag saying whether to do this regional cancellation for a given region.

I think this would be harder to implement than the first few options, though would be feasible. It would require a global communication, but we would only need to do this global communication once per year, so the performance cost should be low. To ensure processor count independence of results, we would want to use shr_reprosum_mod for the global sums rather than relying on standard MPI reductions.

(5) Handle negative DYNBAL ice fluxes through adjustments in the snow pack rather than a negative ice runoff

This idea was motivated by thinking about what we'll want to do when it comes time to tracerize this code, and realizing that negative DYNBAL ice and water fluxes are going to present problems for this tracerization (similarly to #4134 ). I'm focused on negative ice fluxes here because I think they are currently the most problematic fluxes for the CESM3 simulations and also are likely to be the most problematic fluxes for tracerization, but we might want to do something similar for DYNBAL liquid fluxes (and maybe even DYNBAL energy fluxes).

Rather than dealing with DYNBAL ice fluxes via a negative runoff flux, we could instead keep the imbalance in a state variable that indicates an unmet demand; note that this would include the unmet demand both for bulk water and for each tracer (based on the computed imbalance that we're trying to correct with these DYNBAL fluxes). We then keep it there until we're able to meet it locally, by decreasing existing states or incoming fluxes. One easy way to meet this demand would be via earlier or later positive DYNBAL fluxes; strategy (1) from my initial comment would facilitate that. But to avoid relying on that entirely, we would also try to gradually remove some existing state variable to chip away at this deficit state.

I think the best way to do this would be to remove small amounts of the lowest snow layer. This is similar to what is already done with snow capping, so there is precedent - and code - for this already. We might want a minimum snow threshold (in terms of snow water equivalent, snow depth, or number of existing layers) for applying this. Or, maybe better, we could just say something like: remove no more than X% of the mass of the bottom snow layer in a given time step.

A benefit of this approach is that we know the tracer masses in this state, so we can more explicitly and confidently account for the deficits in tracer state by removing tracer mass that we know to exist (as opposed to sending a negative tracer flux that we hope will be met). However, we'll need to figure out what to do when the tracer ratio in the source state (e.g., the bottom snow layer) differs from the tracer ratio in the deficit pool, as will typically be the case. I'm not sure off-hand whether it will be better to remove tracers according to their ratio in the source state or according to the ratio in the demand pool.

I also considered doing this by decreasing the incoming snowfall (rather than by decreasing the bottom snow layer). But I think this could cause more problems. We at least wouldn't want to decrease the snowfall entirely to 0, since that would impact the surface albedo. We could decrease the snowfall by up to some specified fraction. But it seems best to take state out of the snow pack, since removing snow from the bottom layer seems least likely to impact the climate. Also, in typical cases, we'll have a large snow state immediately following the dynamic landunit transition anyway, so this seems like a good source to use to account for this imbalance. And this seems sort of equivalent to removing incoming snowfall anyway... just with some extra patience needed for the snowfall to accumulated into a snow pack before then removing it.

I considered doing something similar with positive DYNBAL fluxes. The closest analog I can think of off-hand would be applying these as a snow flux rather than an ice runoff flux, but it seems like this could mess with the system too much: for example, a gradually-dribbled snow flux will maintain an anomalously high albedo. We could try to apply this to the bottom of the snow pack, but I'm not sure we have an easy way to do that in the current code. So my initial feeling is that it's probably best to continue to handle these positive DYNBAL fluxes as they currently are - via an ice runoff term. And this feels reasonably physically justifiable.

In this draft formulation, note that there is an asymmetry: A positive DYNBAL imbalance would always be fluxed out at the specified rate, but a negative imbalance would only be drawn down if there is sufficient state or incoming flux. I would expect there to typically be sufficient or nearly sufficient state (since a relatively large state variable is what leads to this imbalance in the first place), but that's not guaranteed. So, on average, I would expect these DYNBAL pools to equilibrate at slightly positive values, leading to a small non-physical storage of mass in the system. I'm not sure whether this is something to be concerned about.

My current preferences

I'm currently thinking that we should deal with this via a combination of approaches (1) and (5). Approach (1) will reduce the fluxes considerably, especially in the oscillatory case. This feels important and helpful with (5) (to decrease the amount by which we're adjusting state variables) and without (5) (to decrease the typical adjustment fluxes). Adding approach (5) then helps us deal with the negative ice DYNBAL fluxes locally, which will both (a) put less reliance on the global removal of negative runoff in CMEPS (and consequently, I think, lead to better regional fluxes), and (b) make it easier to handle tracer fluxes.

[samsrabin]

Adding "next" to discuss milestone.