Include option for "plate model" in the "mass conserving temperature"

[lhy11009]

The option we have for the "mass conserving" temperature is to use the hpc (half-space cooling) model in the temperature of the slab. This is a good choice, but if the user chooses to use the "plate" model for their oceanic plate, there will be a abrupt change in temperature around the trench between the slab and the subducting plate. To be more specific, this issue is not so apparent when the subducting plate has a young age, as the two models agree well with a plate age to at least 80 Ma. But with an older subducting plate, this will certainly be an issue. That's why I want to have the "plate" model as an option to use in the "mass conserving" temperature.
Below is an example of using the "plate model" for the oceanic plate and the "hpc model" inside the "mass conserving" temperature model for the slab. Note the line-shape discontinuity in the temperature field below the trench.

Essentially there are two things to implement: one is a temperature formulation to give to the temperature below the slab. And another one is an energy formulation to give to the calculation of the temperature above the slab.

Formulation

The equation for the plate model temperature is:
$T - T_p = \left(T_m - T_p\right)\left(\left(1-\frac{y}{L}\right)-\sum_{i=0}^{\infty} \frac{2}{n\pi}sin\left(\frac{n\pi y}{L}\right)exp\left[\left(\frac{u L}{2 \kappa} - \sqrt{\frac{u^2 L^2}{4 \kappa^2} + n^2\pi^2}\right)\frac{u t}{L}\right ]\right)$
Where $T_m$ is the minimum temperature and the $T_p$ is the background temperature (i.e. the potential temperature) . L is the max depth for the plate model. Compared to the classic expression for the plate model, the surface temperature is replaced by the minimum temperature while the mantle temperature is replaced by the background temperature.
An integration to get the heat contents:
$H\left(u, t, T_m, T_p \right) = \rho c_p \int_{0}^{L} (T - T_p)dy = \rho c_p \left(T_m - T_p\right)\left(\frac{L}{2} - \sum_{k=0}^{\infty} \frac{4L}{(2k+1)^2 {\pi}^2}exp\left[\left(\frac{u L}{2 \kappa} - \sqrt{\frac{u^2 L^2}{4 \kappa^2} + (2k+1)^2\pi^2}\right)\frac{u t}{L}\right ]\right)) $
With these equations. I added an option to use the plate model for the temperature model of "mass conserving". The temperature on top accounts for the differences in the heat contents from the plate cooling, while temperature below the slab is assigned with the plate cooling model directly.
Specifically, the initial heat content is given by
$H_{ini} = H\left(u, t_{sb} , T_s, T_p \right)$
by substituting the velocity and age with the spreading velocity and the age of the subducting plate at the trench. Moreover, the temperature at the surface of the model domain is applied rather than the minimum temperature on the slab surface.
Subsequently, a heat content at the bottom of the slab is given by:
$H_{bot} = H\left(u, t_{eff}, T_m, T_p \right)$
Where the time of an effective plate age is adapted.

Analysis: Differences in the effects on the mantle temperature

The take-home message is the temperature field in the mantle wedge shouldn't vary much when the reference model is switched from the "hpc" to the "plate" model

Here, I defined the "heating thickness" of the thermal model as:
$L_h\left(t\right) = \frac{H\left(u, t_{sb} , T_s, T_p \right)}{\rho c_p \left(T_s - T_p \right)}$
with the aforementioned definition of H. These new value would have the virtue of being independent of the value of $T_s$ and $T_p$ and has a clear meaning. Physically, it's that the heat from this thickness of the mantle is removed through the surface in a period of time.
Of course, the equation for the plate model and the HPC model is different. In the next plot, I should the difference between these two models in the figure on top:
Here the values for the HPC model will keep increasing, which is expected as the HPC model will keep cooling the mantle. On the other hand, the plate model will first increase, and then reach a maximum as the plate keeps aging.
In the 2nd subfigure down below, I include the effects of a warming slab surface (i.e. The minimum temperature is warming than the surface temperature of the model). In this example, the bottom half of the slab is warming up by the mantle wedge, so that the "heating thickness" is negative. The curve would start at value 0 with the age at the trench (40 Ma). Then it ages and warms up till the age at the tip (44 Ma, by a convergence rate of 5 cm/yr, the slab is 200 km long in this case).

The prm file

a new option of "reference model name" is added with the section of the "mass conserving" temperature model. "plate model" (default is "half space model"). Another option of importance is the "max distance slab top", as this value is used in the plate model formulation in this section as the thickness of the plate. In order to render a continuous transition from the plate onto the slab, it should be the same value used in the "plate model" temperature model in the "oceanic plate" section of the sp plate.

Here is an example:

  "temperature models": [
    {
      "model": "mass conserving",
      "density": 3300,
      "thermal conductivity": 3.3,
      "adiabatic heating": true,
      "plate velocity": 0.05,
      "ridge coordinates": [
        [
          [
            0,
            -5.0
          ],
          [
            0,
            5.0
          ]
        ]
      ],
      "coupling depth": 50000.0,
      "shallow dip": 60.0,
      "taper distance": 100000.0,
      "min distance slab top": -100000.0,
      "max distance slab top": 150000.0,
      "reference model name": "plate model"
    }

Results

Results for an old subducting plate of 150 Ma (overriding plate is 40 Ma in age)

In the Cartesian geometry, the results from using a hpc model (left) is compared to a plate model (right). The one with the hpc model has a incontinuous transition at the trench while the one with the plate model ensure a continuous transition at the trench:

In the spherical geometry, the differences between the hpc (left) and the plate (right) models still hold, but both of these models have induce a singularity around the tapering at the tip:

I am going to show more of this issue with a younger plate (40 Ma).

Results for a young subducting plate of 40 Ma (overriding plate is 20 Ma in age)

For a 40 Ma slab, the transferring from the subducting plate to the slab is continuous as expected. Left is a HPC model and right is a plate model:

Same results but in the chunk (spherical) geometry. The issue with the tapering singularity is not as proclaimed as the case of the old subducting slab, but the singularity is still there:

Tests

I have added two tests. One is the "mass_conserving_slab_temperature_plate_model_young" with a subducting plate of 40 Ma at the trench (left sub-figure) and the other one is the "mass_conserving_slab_temperature_plate_model_old" with a subducting plate of 150 Ma at the trench (right sub-figure). These tests are visualized down below:

[codecov]

Codecov Report

Attention: 1 lines in your changes are missing coverage. Please review.

Comparison is base (40655ff) 93.46% compared to head (cfb5408) 93.49%.

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #471      +/-   ##
==========================================
+ Coverage   93.46%   93.49%   +0.03%     
==========================================
  Files          92       92              
  Lines        6301     6346      +45     
==========================================
+ Hits         5889     5933      +44     
- Misses        412      413       +1     

Files	Coverage Δ
...ucting_plate_models/temperature/mass_conserving.cc	97.47% <98.03%> (+0.06%)	⬆️

---

Continue to review full report in Codecov by Sentry.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 40655ff...cfb5408. Read the comment docs.

[github-actions]

Benchmark	Main	Feature	Difference (99.9% CI)
Slab interpolation simple none	1.023 ± 0.005 (s=447)	1.023 ± 0.007 (s=436)	-0.1% .. +0.2%
Slab interpolation curved simple none	1.026 ± 0.003 (s=411)	1.026 ± 0.003 (s=469)	-0.1% .. +0.1%
Spherical slab interpolation simple none	1.257 ± 0.009 (s=377)	1.256 ± 0.005 (s=342)	-0.2% .. +0.1%
Slab interpolation simple curved CMS	1.153 ± 0.010 (s=381)	1.153 ± 0.008 (s=402)	-0.2% .. +0.2%
Spherical slab interpolation simple CMS	1.863 ± 0.005 (s=230)	1.860 ± 0.007 (s=256)	-0.2% .. -0.0%
Spherical fault interpolation simple none	1.273 ± 0.004 (s=356)	1.272 ± 0.005 (s=354)	-0.2% .. -0.1%
Cartesian min max surface	2.376 ± 0.017 (s=175)	2.363 ± 0.021 (s=207)	-0.8% .. -0.3%
Spherical min max surface	7.379 ± 0.020 (s=62)	7.366 ± 0.020 (s=63)	-0.3% .. -0.0%

[lhy11009]

There are the wb and prm file I used to generate the test.
case.wb.txt
case_f.prm.txt

[lhy11009]

Solution for the mentioned issue

Discussed in the issue here:
#364

Change in the parameter file: change the "depth method" from "begin segment" to "begin at end segment":

  "coordinate system": {
    "model": "spherical",
    "depth method": "begin at end segment"
  }

A result of an old slab (150 Ma at the trench) after modification is shown below. The temperature field is continuous at the change of the segment.

[MFraters]

Now that #470 is merged, could you rebase and let me know when it is ready for review or if you need more input?

[lhy11009]

After rebasing, the temperature in my young slab looks nice. But there is a thin thermal boundary in the other part of the model. As I didn't assign any overriding plate there, I am a little confused about it being there. As in the previous visualization of the same test, there was no trace of a thermal boundary in the other part of the model.

@mibillen Is this from some updates on the model or it's something else in the WORLDBUILDER?

[MFraters]

hmm, strange. Which worldbuilder file did you use for this figure?

[lhy11009]

It's the tests/gwb-dat/mass_conserving_slab_temperature_plate_model_young.wb file I included for a young plate. This plot is generated with this grid:
build_mass_conserving_plate_model_test.grid.txt
It's basically the same test, but ran with the version before and after the rebase onto the new update of WB.

[lhy11009]

Maybe you can take a quick look at the files first and see if there is some quick issue you can notice. Then we can together look at the visualization, with these different branches before and after the rebase on my computer.

[MFraters]

Sorry it took me so long to look into the issue. I think the surface temperature is forced by these entries in the input file you mentioned: "surface temperature":273, "force surface temperature":true,.

Can you rebase as well?

[lhy11009]

Hi @MFraters
It took me a while to get back on top of the PR. I have fixed the tests and rebased it, can you take a look again? Thanks.

[MFraters]

This looks good to me. Only a two small comments. @mibillen, can you also take a look at this given that you are more familiar with the mass conversing model?

[mibillen]

I've reviewed the whole code and besides the mispelling of summation I do not see any other issues. One comment however, is that there is a typo in the heat flow equation given at the top of this pull request: the "4" before the exponential term should be "4L" - based on units and me checking the integral. It is correct in the code. Any way to fix that above? I recommend that this can be merged once theses last few points are corrected.

[lhy11009]

I have addressed the comments. And @mibillen is right about the missing "L" in the equations above. I have modified the text.
@MFraters, if you have time, can you take another look at the "enum" variable from a previous comment? I changed the variable "reference_model_name" to enum type at the end of the file "include/world_builder/features/subducting_plate_models/temperature/mass_conserving.h". It's working but I wonder if it's what you want.

[MFraters]

Yes, that was what I meant. Can you fix the tests?

[lhy11009]

@MFraters @mibillen, I have fixed the tests and indentation. Let's give it another try.

[MFraters]

There are some comments marked as resolved, but don't seem to be. Can you take another look at the comments to make sure that they are resolved?

[lhy11009]

Yep, I might mess up using the local folder because I also have a local folder on my laptop. I'll take a look.

…

On Fri, Jan 12, 2024 at 6:51 AM Menno Fraters ***@***.***> wrote: ***@***.**** commented on this pull request. There are some comments marked as resolved, but don't seem to be. Can you take another look at the comments to make sure that they are resolved? ------------------------------ In include/world_builder/features/subducting_plate_models/temperature/mass_conserving.h <#471 (comment)> : > @@ -110,7 +110,7 @@ namespace WorldBuilder bool adiabatic_heating; std::vector<std::vector<Point<2>>> mid_oceanic_ridges; Operations operation; - + std::string reference_model_name; This doesn't seem resolved yet to me. ------------------------------ In source/world_builder/features/subducting_plate_models/temperature/mass_conserving.cc <#471 (comment)> : > @@ -226,6 +231,8 @@ namespace WorldBuilder { const double distance_from_plane = distance_from_planes.distance_from_plane; + const int plate_model_sommation_number = 100; // for the plate model same here — Reply to this email directly, view it on GitHub <#471 (review)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AG2WMFKCFFFDBJTYCCVZ533YOFEYHAVCNFSM6AAAAAAUPPLOW2VHI2DSMVQWIX3LMV43YUDVNRWFEZLROVSXG5CSMV3GSZLXHMYTQMJYGQZTOMBSGA> . You are receiving this because you authored the thread.Message ID: ***@***.*** com>

-- Ph.D student of Geophysics Earth & Planetary Sciences Dept., UC Davis Davis, CA 95616

[lhy11009]

@MFraters Indeed, I have addressed these comments in the local branch on my other laptop. Now that I have fix the branch and tests, it should work.

[lhy11009]

@MFraters, Hi Menno, do we have extra work to do here?

[MFraters]

I asked on slack a bit ago, but you may have missed it between the other messages. Currently your changelog is at the 0.5.0 entry, which is already released. Could you add it to the unreleased section? You probably need to rebase before. Other than that it is good to go.

[lhy11009]

Sorry that I missed that. I vaguely remember you sent something before but I didn't find that on github. I have migrated the change log to the unreleased part.

[MFraters]

Great, thanks for contributing this!

[lhy11009]

No problem! Great job for you to review this too.

…

On Fri, Feb 2, 2024 at 3:38 PM Menno Fraters ***@***.***> wrote: Great, thanks for contributing this! — Reply to this email directly, view it on GitHub <#471 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AG2WMFPGK4AV6GR6MMLHL2TYRV2HPAVCNFSM6AAAAAAUPPLOW2VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTSMRUHEYDOMZTGQ> . You are receiving this because you authored the thread.Message ID: ***@***.***>

-- Ph.D student of Geophysics Earth & Planetary Sciences Dept., UC Davis Davis, CA 95616