Track free water from each lithology in the Tian 2019 Reaction Model

[danieldouglas92]

A hacakthon project. I haven't tested if this works yet, but I wanted some feedback on what I was implementing. This might be a simplified version of what I was initially thinking, but I think at least for my purposes what I'm doing here will work.

Before this PR, using the Tian 2019 reactions for the Reactive Fluid Transport material model only determined the total bound water and the total free water at each point. This PR adds the option for the user to track the volume fraction contribution from each of the 4 lithologies to the total volume fraction of free water in the model domain. This requires the definition of 4 more compositional fields: sediment_porosity, MORB_porosity, gabbro_porosity, peridotite_porosity, which will purely act as tracers that are advected with the total fluid velocity (which is determined based on the porosity compositional field). The simplifications that I'm making that I want an outside perspective on are:

1. The change in these 4 sub-porosity fields are being determined as simply the change in the total porosity multiplied by the fraction of each composition that exists at that point. For example, a 10% (total) porosity change at a point which contains 75% peridotite and 25% sediment, would result in a peridotite_porosity of 7.5%, a sediment_porosity of 2.5%, and 0% for MORB_porosity and gabbro_porosity. I think this is fine, because the total porosity change is already being determined using the average composition at each given point.

2. The problem with the above approach is that the change in the total porosity, which will always be equal to or greater than the values of the subfields, could result in negative values of each sub-porosity field. What I did to get around this is to just prevent the reaction rates for these subfields from becoming negative, so that the values for these sub-fields can only ever increase or remain constant. I think because these fields are meant to act as tracers for where the water is moving, this is ok. An example for why I think this is useeful is imagine you have the 4-subfields advecting with some non-zero total porosity, and then the total porosity goes to 0 because it gets reabsorbed into the solid somewhere in the model, you would be left with the value of the 4 subfields at this point, which could be useful for determining where the free water from each lithology is ending up in the model as a post-processing step.

I'll work on a test case for this to check that it is actually doing what I want it to do over the next week.

[gassmoeller]

Sorry it took a while to get to this PR. I have not reviewed this in depth yet, let me know when you are ready for that.
However to your questions:

1. I can see this work for models that use discrete composition distributions (where most cells are completely filled with only a single composition), but for mixed cases like the one you mentioned it seems physically very unlikely. The fluid is released by the breakdown of certain minerals, which may be purely existing in one of the compositions in a cell (say clays in sediments, or serpentine in oceanic crust). So if your cell is filled with a mixture of sediments and dry peridotite, and you see an increase in free water, it is likely purely sediment derived fluid, not derived from the peridotite. Do we have information somewhere in the code about which of the compositions release the fluid, or are we limited to the aggregate fluid release?
2. I do not understand why the reduction in global fluid content could cause negative individual fluid contents. Lets say you have a global fluid content of 10% and individual contents of 7% (sediment derived fluids) and 3% (MORB derived fluids), and the global content is reduced by 5%. Then you wouldnt subtract 5% from both of the individual fluid fields, you would distribute the reduction onto the two fields according to their relative fraction (7% -> 5%/10% * 7%; and 3% -> 3% * 5%/10%). This should never give you negative values (except for possible oscillations). Could you elaborate how you reduce the individual fields? I am pretty sure never reducing the values will give you wrong values in the long run (what happens if you have repeated influx of fluids into a region where part of them get absorbed? Will the individual fractions of the free fluid still be correct? One way around this problem while still keeping track of absorbed fluid is to introduce another set of individual fields that tracks the bound portion of the individual fluid fields for each composition. It gets a bit unwieldy with all these fields, but maybe that is the more correct approach?