Calculation of mantle mass by Aragog for PROTEUS

[nichollsh]

PROTEUS expects that the interior modules are able to return values for total mantle mass, as well as the liquid and solid components.

This is handled by SPIDER internally, and read from the JSON files. This is currently not handled by Aragog internally, and is done in aragog.py here: https://github.com/FormingWorlds/PROTEUS/blob/b0d2ae0058dc799299687235dd953ba0ba750256/src/proteus/interior/aragog.py#L214

Note that this calculation reflects the Dummy module, but I do not think it is correct for a resolved case. The global melt-fraction (Phi_global) is calculated using the volume-averaged melt fraction in each cell (this is fine). However, using this value as in the line above is only true for constant density, which can be shown analytically.

When density is not constant the mass of melt should be calculated as the sum of the melt fractions and masses of each layer in the Aragog mesh. The mass of solid can still be calculated as the remainder of this and the total mantle mass.

@djbower @planetmariana @lsoucasse Please can you tell me if this logic is correct? If so, we should think about updating the calculation M_mantle, M_mantle_liquid, and M_mantle_solid used by PROTEUS/Aragog.

[lsoucasse]

You say that Spider handles liquid and solid mantle mass calculation internally but is it really different from what we are doing in Aragog now? We could compute liquid/solid mantle mass from liquid/solid density fields but then it will not be consistenty with the total mantle mass and pressure field so I am not sure if it is better than relying on the melt fraction.

[djbower]

The design philosophy behind SPIDER/Aragog was the following:

• The mass of the planet is constant and set during initialisation.
• Setting the mass and an initial P-T relationship creates a fixed mapping between radius-P which is useful because thermodynamic quantities are often expressed in terms of P, but to compute the energy fluxes through the mesh requires radius.
• The liquid/solid properties are used to compute material properties and determine energy terms (fluxes and interior contributions), but these are not deemed significant enough to affect the planet's overall volume, which would change the cell boundaries and thus the energy flux through the spherical shell.

The choices above necessarily introduce several inconsistencies that are absorbed elsewhere in the code, and which are fine for Earth-sized planet but will invariably break down as the planet size increases. For example, the initial assumed Adams-Williamson EOS will not produce a density that exactly agrees with the solid, liquid, or mixed phase. Equivalently, the mass of each cell is effectively different from the perspective of the fixed mesh versus the actual amount of melt and solid that is accommodated. It can sometimes be thought of as a compromise between preferencing the calculation to conserve mass or conserve volume. Again, as the planet size increases this choice will be more difficult to make, and it may necessitate re-deriving the mass-volume-pressure-radius relationships from the static and dynamic perspective to deduce what the long-term best compromise is (short of solving a static-structure problem every time-step, which I think we've agreed is not feasible at the moment).

I will add some notes on mass coordinates in the Aragog GitHub issues in the next few days, which also relates to this problem.