ladjust_bury_coeff computation should take restoring into account

[klindsay28]

I have a 300 year long run with ladjust_bury_coeff=.true., and global inventories of SiO3 and PO4+DOP+DOPr are drifting. While the drift is small, it is not decaying, which is the expected behavior with ladjust_bury_coeff=.true.. I'm pretty sure that the problem is that the ladjust_bury_coeff computation is not taking restoring, which changes tracer inventory, into account. This suspicion is supported by the observation that the drift in global inventory matches the globally averages of the restoring tendencies of SiO3 and PO4 respectively. So I've got to fix this for the spinup.

[klindsay28]

Thinking about this some more, my proposed solution to the restoring induced drift, including restoring in the ladjust_bury_coeff computation, has a drawback.

In our standard CESM configuration, we are restoring only in the marginal seas. We do this to deal with the imbalance between riverine inputs and loss to sediments in marginal seas. If we do what I suggest, then we are letting imbalances in the marginal seas impact the open ocean. I think this is undesirable.

An alternative is to adjust bury coefficients based on averages that are computed using only open ocean points. I think this would decouple the marginal seas completely from the open ocean. A simple implementation is to have POP compute all global MARBL requested averages. POP doesn't really know details of the fields being averaged, so all or nothing seems to be the simplest approach. However, I'm not entirely satisfied with this simple approach, because POP is then not really doing what MARBL is requesting it do, which is to compute global averages. POP is sorta saying, I know what's best for MARBL.

I'm not really satisfied with either solution. Thoughts? Suggestions?

[matt-long]

@klindsay28 I think your proposed solution to adjust bury coeff's based on open ocean points is the right answer from a scientific point of view. When you say

A simple implementation is to have POP compute all global MARBL requested averages.

I infer that you mean POP will compute "global" averages inclusive of only open-ocean points and that's simple to implement.

I think I understand your discomfort with this approach from the perspective of the interface semantics. I think the response to this is to clearly document our rationale. The GCM is computing a "global" average across interconnected basins because this is where we need to satisfy the integral mass balance constraint (given restoring in the marginal seas).

[klindsay28]

Thanks @matt-long, I was leaning towards my 2nd proposed solution, which you are preferring.

The change is literally a 1-line mod to code (plus more for comments).
I'll work with Mike this afternoon to get that in.
The change is to POP, so I'm changing this issue to reflect that.

This will only change answers in ladjust_bury_coeff=.true..

[matt-long]

Sounds good!

[klindsay28]

I've been thinking more about this (gotta stop doing that)...

I think that doing bury_coeff adjustment in the open ocean is fine for nutrients.

However, for carbon, I think that is (slightly) problematic. If there is an imbalance of carbon in the marginal seas between riverine inputs and loss to sediments, this will result in a net surface CO2 gas flux at equilibrium. One of the primary reasons for adjusting the bury coefficients is to avoid this. (The other is to eliminate drift in nutrient inventories.)

I just checked DIC_RIV_FLUX+DOC_RIV_FLUX+DOCr_RIV_FLUX-pocToSed-calcToSed integrated over the marginal seas is a long cplhist run, and see that it is ~+0.0046 PgC/y. So we would be outgassing that. The riverine input is ~0.0116 and the burial is ~0.0070. A caveat with these numbers is that they are from a run where the bury coefficients were computed using global ocean averages, not open ocean. But hopefully they indicate well the order of magnitude of the imbalance that I'm bringing up.

Some options:

1. do nothing, live with the nutrient drift
2. as proposed, compute bury coeffs using only open ocean integrals, and live with outgassing (it's kinda small)
3. like 2), but also reduce riverine C input to marginal seas by 40%, to reduce the imbalance
4. use global ocean for all bury coeff adjustment calculations, and include nutrient restoring in the computation
5. use global ocean for C bury coeff adjustment, and open ocean for P and Si bury coeff adjustment

From a science point of view, I think 5) seems best. But it is also the most invasive, and I think it would require adding something to the MARBL-GCM interface so that MARBL could inform the GCM over which domain the requested global averages are to be computed.

Personally, I'm leaning towards 1) or 2), primarily because of practicality, but I think I could be swayed to another option.

Thoughts?

[mnlevy1981]

5. use global ocean for C bury coeff adjustment, and open ocean for P and Si bury coeff adjustment

From a science point of view, I think 5) seems best. But it is also the most invasive, and I think it would require adding something to the MARBL-GCM interface so that MARBL could inform the GCM over which domain the requested global averages are to be computed.

I agree with this being the most appealing solution from a "getting the science right" perspective, and also with the concerns about the interface complexity. It seems like we would need to introduce a global averaging datatype to replace

     ! public data - global averages
     real (r8)                                 , public, allocatable  :: glo_avg_fields_interior(:)   ! output (nfields)
     real (r8)                                 , public, allocatable  :: glo_avg_averages_interior(:) ! input (nfields)
     real (r8)                                 , public, allocatable  :: glo_avg_fields_surface(:,:)  ! output (num_elements,nfields)
     real (r8)                                 , public, allocatable  :: glo_avg_averages_surface(:)

This new datatype would also communicate to the GCM whether the average should include marginal seas, and then it's a matter of determining how to decide whether or not the average should include marginal seas (maybe as easy as "if C bury coeff include marginal seas otherwise do not", maybe as complicated as letting the user control it somehow or basing it off of what tracers are being restored)

Maybe we should leave things as they are for CESM 2, but try to implement the above for the CMIP runs? (Is that practical, or are we going to use the CESM 2 spin-up to determine the bury coefficients for CMIP?)

[matt-long]

Would another solution be to turn burial off in the marginal seas, zero-out riverine C inputs, and do all balancing based on open ocean?

This would require extending the interface to enable a burial switch.

[klindsay28]

That is another option, but I don't see a straightforward way to turn burial off in the marginal seas. Burial occurs in the marbl_set_interior_forcing calling tree of MARBL, and that part of the code currently doesn't know if a column is in the open ocean or in a marginal sea.

[klindsay28]

didn't read completely.

I'd rather not have the GCM know about burial.

[matt-long]

Fair point. However, I could see extending the marbl_domain_type to include something like a basin field; burial could then be turned off optionally within MARBL depending on the value of this field.

[klindsay28]

I'd rather not go that route. Changing the nature of riv fluxes and burial based on basin seems more kludgy to me than specifying which basin averages are to be computed over. So I prefer 5) over your 6). This is admittedly subjective.

[matt-long]

I don't really see how option (6) is inherently more kludgy, but I am fine with either (5) or (6).