Process subcycling and time-averaged diffusion

[sdeastham]

@chinahg noted that contrail behaviour was non-convergent with changing time steps. This plot shows, for the issl_rhi140 simulation (restricted to 2 hours of run time), how the ice mass (upper) and ice crystal number (lower) vary with changing timesteps - either all set to 1 minute (alldt1), all set to 5 minutes (alldt5), or with transport set to 5 minutes and ice growth set to 1 minute:

Increasing the transport timestep causes the initial drop in ice crystal number to disappear, while also resulting in slower growth in the ice mass. Having disparate growth and transport steps results in what appears to be a totally different contrail. Subsequent investigation reveals at least three causes:

1. Diffusion coefficients are estimated for the full transport step based on the value at the time step mid point: this results in inconsistent representation of diffusion at the start of the simulation as the diffusion coefficients decay towards their final value (after 7.5 minutes). This is particularly problematic if non-linear decay is used, or if the transport time step exceeds 7.5 minutes.
2. Reducing the ice growth timestep below the transport timestep resulted in non-physical behaviour: until Timestep choice in input.yaml does not do what you would expect #107 was resolved, making the ice growth time step shorter than the transport time step would result in transport being erroneously overestimated (transport would occur every ice growth time step, but using the transport time step length!).
3. Ice growth and transport are fast and tightly coupled: we allow users to separate ice growth and transport time steps, but this causes endless difficulties because the two processes are actually tightly coupled.

To resolve this, multiple changes are implemented. First, timestep settings were corrected in #122 (resolving #107). This corrected the issue with split timesteps, but changing the time step still means that the ice crystal number still does not fall correctly and the ice mass growth is slowed.

Implementation of a semi-Lagrangian solver (#125) changes the reference (alldt1) case somewhat, but does not address the fundamental issues - it is included here only to show that the additional numerical diffusion due to the Eulerian solver did not drive the base issue, although the move to an S-L scheme does prevent spurious sublimation of small crystals:

I then decided to take the step of changing how we approached the problem. The fundamental issue is that the architecture isn't really designed to allow decoupling of the ice crystal growth and transport. What we really want is for the ice crystal growth to be performed at multiple points during transport, but that in turn requires that the transport also be sub-stepped. As such, I decided to simplify things: rather than specifying a transport time step and a growth time step, you now specify an "outer time step" (the time step on which met data is updated and regridding is performed), and an "inner physics time step" (on which both transport and ice growth are calculated). This means that, no matter what the outer time step is, the actual transport and physics are performed using the inner step. In the following, the outer time step is "tdt" while the inner is "idt". The additional "fine" case uses an inner step of 30 seconds.

This at least shows that the solution is similar whether a substep of 30 or 60 seconds is used, and whether an outer step of 1 or 5 minutes is used; we are no longer seeing the drastic divergences observed earlier. This also means that increasing the outer time step alone has little effect on run time. The good news is that a 2-hour simulation with 1 minute substepping is now 16% faster than the reference case; if the simulation is allowed to extend to completion, the average rate (in terms of simulation hours completed per wall clock minute) has increased by 50%. To really assess convergence though a more thorough study (ideally involving consideration of time resolution, grid resolution, and aerosol bin resolution) is needed.

[sdeastham]

FYI - this is on hold pending documention and a grid resolution study (see #125).

[sdeastham]

@lrobion @speth @Calebsakhtar @chinahg - I've drafted a report summarizing APCEMM in its current state, including convergence testing and proposed simulation configurations for the broader user base. I would be very grateful for your thoughts on this; the report was developed collaboratively with Google Gemini (I'm not convinced this ended up being any faster than just doing it myself after all the edits and reiteration, but it was at least a more interesting process):
APCEMM v1.3 documentation

[Calebsakhtar]

Thanks so much for this @sdeastham! I have scheduled some time on Thursday for this review, so you should expect to hear back from me on or before Friday.

[Calebsakhtar]

Review of the science and wording completed:
APCEMMDocumentation-CAM27082026.pdf

No major issues found. I did not check the equations extremely thoroughly this time, but the ones I am familiar with all looked good.

I think this will be extremely valuable for the wider contrails community, thanks so much @sdeastham!

[speth]

A couple comments on the report: I had a little trouble understanding the connection between the description of the "Strang Splitting" that's used to couple the advection and diffusion operators (5.2) and the description of the "Inner Microphysics Subcycling" loop in 5.3. The latter describes the inner loop as:

• growth / sublimation
• sedimentation and diffusion
• horizontal advection

However, I believe the actual process in the code is:

• Strang split advection + diffusion, where each half advection step is horizontal advection followed by sedimentation
• growth / sublimation

Or, written out more fully:

• 1/2 step horizontal advection
• 1/2 step sedimentation
• full 2D diffusion step
• 1/2 step horizontal advection
• 1/2 step sedimentation
• full growth / sublimation

Also, I remember a discussion we had at some point about it being problematic for particles to potentially fall multiple layers in a single advection step, without any intervening growth step, since the particle should take up water from those intermediate layers. I may have hallucinated that you had implemented something to handle this. Are the cases run with the inner substep of 10 s short enough that this doesn't happen?

It might be interesting to compare the 2D profiles (ice crystal number density or other) at a couple points in the lifetime for your four different "tiered" profiles, to see how they look. I wouldn't anticipate anything too wild, given the relatively good agreement in terms of the integrated optical depth over time.

[sdeastham]

Thanks @speth and @Calebsakhtar for the comments!

@Calebsakhtar: really appreciate the thorough review! I've made minor modifications accordingly where necessary. I would consider anyone who has provided input on the report to be a co-author. I decided not to include some of the suggestions not because I thought they were incorrect, but rather because I think they are priorities for next steps; I'm loath to let v1.3 slide much further as I think some of these updates really need to get out there.

@speth: I took a look through Section 5.2 and 5.3, comparing to the code - I agree with your assessment. I've updated the report accordingly. I'm also intending to add the proposed comparison of 2D cross sections (I'm curious too, if nothing else!).

@lrobion and @chinahg - I'd still very much value any thoughts you might have on the report! One big question I have is whether we need to roll in a solution to #115. I'm very reluctant, because I don't know how much time I'm likely to have to incorporate this in the next few weeks. However I'm open to other perspectives!

[lrobion]

Will have feedback for you by tomorrow EOD (US). Looking into #115 again, I do think it's worth fixing, but the impact is smaller than I thought. Given that APCEMM diffusion coefficients for ice are binary (they are uniform for all other quantities), the issue only occurs at the boundary between the contrail mask and the ambient air. In addition, the diffusion coefficient values inside the mask only differ due to diffusion enhancement which with its current profile reverts to the background value after 780s. I am hoping to get the chance to do a drive by fix for this before next week.

I would also like to assess the impact of #106 before deciding if that's a must fix for v1.3.

[lrobion]

APCEMM-v1.3-Louis.pdf

Thanks for putting this together Seb, here's my annotations on the doc. This is super helpful. I flagged a couple of things I thought were not 100% clear (differences between PDF and code where I don't know if this is something that is slightly unclear in the PDF or an issue in the code).

Otherwise, I second Ray's comments on looking at the 2D slices. I was also interested in seeing the effect of adding even more bins (order of 500-1000) and seeing the effect it has on contrail simulations for your different tiers. I know performance is a big limiting factor for this right now, but I think there are a lot of ways in which we could make APCEMM faster so these high-fidelity cases might become within reach for regular simulations at some point.

EDIT: updated the PDF with more comments on the TVD section

[lrobion]

I know this is still a draft and depends on #125 which has also been updated so I just put a couple of high level comments for now.

[lrobion]

We also need to add the ice size bins options here r_min, r_max, and the ratio (+ remove them from Core/Parameters.hpp so that they is a single source for the values)

[lrobion]

I am against setting defaults in the header files, I think all defaults should live in defaults/input.yaml, then its the YAML parser's job to reconcile inputs together

[lrobion]

I think we can validate at input parse that this value is > 0 and crash APCEMM otherwise.
This seems to fall back to a hardcoded 60s default value which I don't think we should do (we can put that default in the input file if we want instead)