Reducing circular modeling and clarifying production flow for new parameterizations in ClimaLand.Snow

[a-charbon]

Is your feature request related to a problem? Please describe.
Writeup of a circular modeling phenomenon, requested by @kmdeck :
When creating and implementing new parameterizations for snow depth or density to be used by the Snow model, the obeyed equation is

$$\rho_{water} * SWE = \rho_{snow} * z $$

where z is the height of the snowpack, and SWE is the snow-water-equivalent (mass of water in the snow column, or the height if the snow column was melted into liquid water), and $\rho_{snow}$ is the snow bulk density.

SWE is always a prognostic variable in the model, so modeling $\rho_{snow}$ and $z$ is a matter of choosing/defining a parameterization for one, and deriving the other. Different phenomena dictating the evolution of SWE and energy (albedo, runoff, interactions with other model elements etc.) will depend on $\rho_{snow}$ and $z$ directly, so it is important to be able to reference/call both when necessary.

The original model used a constant density model (a constant value for $\rho_{snow}$), meaning whenever depth was needed, one merely needed to return a scaled value of SWE using a snow_depth() function (which looks at SWE and scales it by the constant $\rho_{water} / \rho_{snow}$, which is > 1 ). This provides the physical consistency $Z \to 0$ when $SWE \to 0$ and guarantees Z > SWE in a manner that is numerically stable even for SWE near float precision and as $SWE \to 0$. This is not the most accurate parameterization, however, and more dynamic models were necessary (preserving this stability and physical consistency would also be useful when introducing such dynamics).

To do this, one can either treat either density or depth as a diagnostic/auxiliary (or prognostic) variable, prescribe the update rule for that variable, and derive the other from the derived value and SWE. Implementing the back-end control-flow of the Snow model to handle all possibilities of user choice (whether they prescribe depth, or density, or both) in a clean and consistent way has proven convoluted. Currently, density is always a diagnostic variable (p.snow.ρ_snow always exists) and SWE is always a prognostic variable (Y.snow.S), we let the user choose whether depth is a prognostic variable or not (specify whether Y.snow.Z exists or not), and have four functions:

• the scalar function snow_bulk_density(z, SWE) which returns SWE/z * ρ_water for floats SWE, z (i.e. it is mapped over a field)
• snow_depth(...) which has default return ρ_water .* Y.snow.S ./ p.snow.ρ_snow but takes the model type as an argument (the arguments "...") so the user can dispatch off it and extend it to just return Y.snow.Z if they use a model that includes it
• update_density!(...) within update_aux!(), which has default behavior
  p.snow.ρ_snow .= snow_bulk_density.(Y.snow.S, snow_depth(...)) but takes the model type as an argument (the arguments "...") so the user can dispatch off it and extend it to their created parameterization
• update_density_prog!(...) within compute_exp_tendency!() (created to handle potential inclusion of prognostic variables associated with density/depth parameterizations, like dY.snow.Z), which has the default return of nothing but takes the model type as an argument so the user can dispatch off it and extend it to update Y.snow.Z with their intended parameterization, should they choose a case where it exists.

This way a user can run a custom snow model and try out their own parameterizations without ever having to interact with (or redefine) the Snow model back-end, though this leaves them a convoluted flowchart to follow depending on how they want to define their model and pick/design their parameterizations (do they extend snow_depth() or update_density!(), or even both? do they need to extend update_density_prog()? See the NeuralSnow extension starting lines 106 to see an implementation of a snow depth parameterization that only needs to redefine snow_depth(...) = Y.snow.Z and not update_density!() but must also redefine update_density_prog!(), or the default ConstantDensityModel only has to extend update_density!() and nothing else)

If we do not do this dispatch option, then internal functions that define effects requiring depth/density face potential redundancy or instability issues. For instance, the snow runoff parameterization (which dictates SWE tendencies) requires the snow depth, so for all parameterizations like this one, we could:

• A) define each function per model type, and the back-end determines the right version by dispatch (one for models that include depth and one for those which do not include depth)...but this approach would require dispatching every depth-dependent snow parameterization, which will be exhaustive, bulky, and not-scalable as more complex parameterizations are integrated or users try to create their own
• B) only ever define these functions to take SWE and $\rho_{snow}$ (since both always exist)...then $z$ is derived within the function as SWE/z * ρ_water, but this is a redundant computation and wasted resources for models where Y.snow.Z already exists as a field that could just be called/referenced. This choice can also be unstable as SWE, Z $\to$ 0 - for model choices where Y.snow.Z and Y.snow.S exist and become small, then updates follow p.snow.ρ_snow = ρ_water .* Y.snow.S ./ Y.snow.Z (numerically unstable) and immediately afterwards, the runoff calculates depth from only SWE and $\rho_{snow}$ as z = ρ_water .* Y.snow.S ./ p.snow.ρ_snow, meaning
  = ρ_water .* Y.snow.S ./ ( ρ_water .* Y.snow.S ./ Y.snow.Z)
  so two unstable calculations were carried out (when one could alternatively just return Y.snow.Z). A similar issue could happen when using models where Y.snow.Z does not exist and a non-constant density parameterization is used (calculating z from SWE and ρ, then using this z to get ρ from SWE and z).

Dispatching the depth and density themselves (so such functions receive the dispatched input snow_depth(...), etc.) removes the redundancies in A) and B) but creates the convoluted development path described above for users trying to define new parameterizations.

Are there ways to mitigate these tradeoffs and create a coherent development flow that is concise and stable?

Describe alternatives you've considered
We could require that snow depth is always a prognostic variable, though this would create some extraneous storage/computation for simple models like constant-density models that do not require tracking of depth. Any inclusion of snow depth of a prognostic variable modeled separately also requires additional handling to make sure we always have z > SWE and $z\to0$ when $SWE \to 0$ and avoid any numerical instabilities. This could ultimately limit the types of snow models a user could have the chance to create without redefining the Snow model back-end, as well.