Initial implementation of past spreading history in mass conserving slab thermal structure

[danieldouglas92]

Began an implementation into the oceanic plate and subducting plate thermal plugins, in response to #529. The idea is that over time spreading ridges are not stationary, and as they move closer/further from the trench the age of the plate entering the trench decreases/increases. This impacts what the state of the thermal of the slab is. I have an initial implementation that is not done, but has promising results. I've attached an image that shows an example where the ridge has moved farther from the trench, increasing the age of the plate as it subducts for 20 Myr. The first shows an example where the ridge does not move relative to the trench, and the image shows isotherms for cold slab temperatures. The second image shows the case where the ridge has moved closer to the trench, resulting in a younger, hotter slab at the trench, and the third image shows the reverse case where the ridge has moved farther from the trench resulting in a colder, older slab at the trench. I've also attached an example .wb file.

This needs to be extended to account for the temperature in the unsubducted oceanic lithosphere as well, because as is evident in the pictures there is a clear discontinuity at the trench, but I think this is close to being completed, maybe 60-70% done.

ridge_drift.txt

[danieldouglas92]

@lhy11009 @mibillen it would be great to get your feedback on this whenever you're able!

[danieldouglas92]

Did some more testing and realized that without a schematic diagram showing what I'm trying to do this code is probably pretty hard to dive into.

The idea is to move a Mid-ocean ridge based off of the difference between how fast the plate is being subducted at the trench, and how fast the plate is being generated. If the plate is being generated faster than it is being subducted, this would mean that the ridge is being pushed away from the trench. If the plate is being subducted faster then the plate is being generated, this means that the ridge is being pulled toward the trench. The mass_conserving temperature plugin would then require the user to input not just the spreading rate, but also the subducting rate. If these are equal, then this scenario resolves to what worldbuilder currently assumes, which is that the ridge is stationary through time and the subducting rate equals the spreading rate. When the ridge is moving, the user specifies a time over which the ridge has been moving. This is used with the velocity of the ridge to calculate the distance that the ridge has moved over that time. The distance is then used to calculate a shift in the age of the trench, which is proportional the spreading rate, and also a shift in the effective age of the slab, which is proportional to the subducting rate.

This produces a shift in the temperature of the slab that is proportional to the duration of ridge movement, the magnitude of the ridge velocity, and also the direction of ridge movement. Here is a case where the ridge moves towards the trench, we see that the isotherms near the trench are shallower, indicating the the slab entering the trench is now younger. Temperature contours are from 300 - 1000 K in intervals of 100 K.

And here is a case where the ridge moves away from the trench, this time we see that the isotherms are deeper, indicating that the age of the slab entering the trench is now older.

I still need to apply the shift in the age_at_trench to the unsubducted portion of the plate, and also create tests and possibly adjust the parsing function get_values_at_array(), but I think this is looking much better.

[lhy11009]

Why the subducting velocity is stored in the .first in some cases and in the .second in other cases?
Is it possible to make them always in .first or always in .second

[danieldouglas92]

Yes I agree that this is confusing and I plan on making this more consistent in the next commit!

[lhy11009]

My main comment is: is there a way to walk around the effective_age_shift ? this variable is not very intuitive and it takes me several times to learn its meaning.
At the same time, I am confused about why we need the time_of_ridge_migartion in the first place.
Maybe it's possible to use the subducting time and:
maybe we can define a parameter of ridge distance as:
ridge_distance = ridge_velocity * time_of_subduction
then the age_at_trench could be modified as
"age_at_trench = ridge_parameters[1] / spreading_velocity + ridge_distance / ridge_velocity

[danieldouglas92]

The names of the variables could definitely use some improvement. But what I intend effective_age_shift to be is a value that increases/decreases the plate age after it has been subducted, white trench_age_shift is a value that increases/decreases what the age of the plate is at the trench, these are different variables because trench_age_shift depends on the spreading velocity, while effective_age_shift depends on the subduction velocity.

The idea for why I made time_of_ridge_migration a variable is so that the user can define how long the ridge has been moving, but I'm certainly open to ideas! It's not clear to me what you mean by time_of_subduction, though. The time that subduction started?

[lhy11009]

Isn't that just the plate_age_sec + (distance_along_plane / subducting_velocity)? So the difference to the plate_age_sec is that the material is subducted so an age is added?

[danieldouglas92]

I renamed plate_velocity to subducting_velocity so it was clearer to me which velocity was which while writing this code, but effective_age_shift is an additive age yes. The difference between distance_along_plane / subducting_velocity and effective_age_shift is that effective_age_shift only applies over a certain length of the slab

[lhy11009]

So a further confusion I notice is this variable of "plate_age_sec".
The "effective_age_shift" is added to the expression here. But the "plate_age_sec" is from the "age_at_trench * seconds_in_year." The "age_at_trench" is from "ridge_parameters[1] / spreading_velocity + trench_age_shift". So in the end, both the "trench_age_shift" and the "effective_age_shift" are added in the expression of effective_plate_age.
This I don't think make sense.

[danieldouglas92]

Ahh yes you're definitely right. I need to be applying that to the unsubducted oceanic plate so that the ages match up, not the subducting plate!

[lhy11009]

@danieldouglas92 Thanks for the explanation, but I am still confused why the ages are modified in case distance_along_plane < (d1 - d2), what happens if distance_along_plane > (d1 - d2)?

[danieldouglas92]

@danieldouglas92 Thanks for the explanation, but I am still confused why the ages are modified in case distance_along_plane < (d1 - d2), what happens if distance_along_plane > (d1 - d2)?

Nothing happens if distance_along_plane > (d1 - d2). To get something to happen to the whole slab the time_of_ridge_migration would have to increase and/or the ridge_velocity would have to increase, it all depends on how much the ridge has moved relative to how long the slab is

[lhy11009]

Is that because if distance_along_plane > (d1 - d2), the point wasn't returned to the surface when looking at t = t2 in the diagram? I guess I am slow to get the idea, but can you write a little more of this criteria? Thanks

[lhy11009]

If it fits into your schedule, maybe we can meet and discuss some time tomorrow or Tuesday? I think Magali will be very busy, but I'll try to touch base with her on Tuesday afternoon.

[danieldouglas92]

I guess the way that I've been thinking about it while working on this implementation is that in that schematic diagram that I drew, we have t1 < t2, where t1 is the present. For times t1 <= t <= t2, the ridge is moving, but for t > t2, the ridge is not moving, so for all times before a user defined time, the ridge was stationary.

If we could meet on Wednesday that would be better for me, My schedule is completely open for Wednesday (Eastern Time)!

[codecov]

Codecov Report

Attention: Patch coverage is 96.00000% with 7 lines in your changes are missing coverage. Please review.

Project coverage is 98.23%. Comparing base (c67addd) to head (3a5fa35).
Report is 13 commits behind head on main.

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #654      +/-   ##
==========================================
- Coverage   98.30%   98.23%   -0.07%     
==========================================
  Files         107      107              
  Lines        7357     7464     +107     
==========================================
+ Hits         7232     7332     +100     
- Misses        125      132       +7     

Files	Coverage Δ
...eanic_plate_models/temperature/half_space_model.cc	98.52% <100.00%> (+0.09%)	⬆️
...es/oceanic_plate_models/temperature/plate_model.cc	100.00% <100.00%> (ø)
...ucting_plate_models/temperature/mass_conserving.cc	99.59% <100.00%> (+0.02%)	⬆️
source/world_builder/utilities.cc	97.22% <98.33%> (-0.14%)	⬇️
source/world_builder/parameters.cc	97.40% <87.50%> (-0.58%)	⬇️

... and 6 files with indirect coverage changes

---

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 c67addd...3a5fa35. Read the comment docs.

[github-actions]

Benchmark	Main	Feature	Difference (99.9% CI)
Slab interpolation simple none	1.012 ± 0.004 (s=470)	1.020 ± 0.008 (s=418)	+0.7% .. +1.0%
Slab interpolation curved simple none	1.023 ± 0.004 (s=466)	1.025 ± 0.006 (s=415)	+0.0% .. +0.3%
Spherical slab interpolation simple none	1.166 ± 0.007 (s=397)	1.171 ± 0.006 (s=376)	+0.3% .. +0.5%
Slab interpolation simple curved CMS	1.052 ± 0.010 (s=418)	1.054 ± 0.009 (s=439)	-0.0% .. +0.4%
Spherical slab interpolation simple CMS	1.546 ± 0.009 (s=293)	1.548 ± 0.009 (s=291)	-0.0% .. +0.3%
Spherical fault interpolation simple none	1.173 ± 0.008 (s=367)	1.173 ± 0.006 (s=402)	-0.1% .. +0.2%
Cartesian min max surface	2.304 ± 0.017 (s=199)	2.318 ± 0.021 (s=193)	+0.3% .. +0.9%
Spherical min max surface	7.219 ± 0.039 (s=60)	7.222 ± 0.029 (s=67)	-0.2% .. +0.3%

[MFraters]

Nice work! And thank you @lhy11009 for your reviewing of the code :)

I have a few initial comments, which mostly stem from valgrind tester. Let me know if you want to discuss them further.

[MFraters]

valgrind is tripping over this line, can you check if this is actually filled?

valgrind output for gwb-dat tests/gwb-dat/cartesian_multiridge.wb tests/gwb-dat/cartesian_multiridge.dat --limit-debug-consistency-checks

...
39: 1500e3 500e3 0 10e3 1538.93  0
39: ==4763== Invalid read of size 8
39: ==4763==    at 0x3C8AFC: WorldBuilder::Utilities::calculate_ridge_distance_and_spreading(std::vector<std::vector<WorldBuilder::Point<2>, std::allocator<WorldBuilder::Point<2> > >, std::allocator<std::vector<WorldBuilder::Point<2>, std::allocator<WorldBuilder::Point<2> > > > >, std::vector<std::vector<double, std::allocator<double> >, std::allocator<std::vector<double, std::allocator<double> > > >, std::unique_ptr<WorldBuilder::CoordinateSystems::Interface, std::default_delete<WorldBuilder::CoordinateSystems::Interface> > const&, WorldBuilder::Objects::NaturalCoordinate const&, std::vector<std::vector<double, std::allocator<double> >, std::allocator<std::vector<double, std::allocator<double> > > >, std::vector<double, std::allocator<double> >) (utilities.cc:1337)
39: ==4763==    by 0x37232B: WorldBuilder::Features::OceanicPlateModels::Temperature::PlateModel::get_temperature(WorldBuilder::Point<3> const&, WorldBuilder::Objects::NaturalCoordinate const&, double, double, double, double, double) const (plate_model.cc:168)
39: ==4763==    by 0x35DF6D: WorldBuilder::Features::OceanicPlate::properties(WorldBuilder::Point<3> const&, WorldBuilder::Objects::NaturalCoordinate const&, double, std::vector<std::array<unsigned int, 3ul>, std::allocator<std::array<unsigned int, 3ul> > > const&, double, std::vector<unsigned long, std::allocator<unsigned long> > const&, std::vector<double, std::allocator<double> >&) const (oceanic_plate.cc:176)
39: ==4763==    by 0x3CE692: WorldBuilder::World::properties(std::array<double, 3ul> const&, double, std::vector<std::array<unsigned int, 3ul>, std::allocator<std::array<unsigned int, 3ul> > > const&) const (world.cc:387)
39: ==4763==    by 0x31450E: main (main.cc:305)
39: ==4763==  Address 0x5f82118 is 0 bytes after a block of size 8 alloc'd
39: ==4763==    at 0x4849013: operator new(unsigned long) (in /usr/libexec/valgrind/vgpreload_memcheck-amd64-linux.so)
39: ==4763==    by 0x43341D: __gnu_cxx::new_allocator<double>::allocate(unsigned long, void const*) (new_allocator.h:127)
39: ==4763==    by 0x4265D3: std::allocator_traits<std::allocator<double> >::allocate(std::allocator<double>&, unsigned long) (alloc_traits.h:464)
39: ==4763==    by 0x412877: std::_Vector_base<double, std::allocator<double> >::_M_allocate(unsigned long) (stl_vector.h:346)
39: ==4763==    by 0x3FFBD8: void std::vector<double, std::allocator<double> >::_M_range_initialize<double const*>(double const*, double const*, std::forward_iterator_tag) (stl_vector.h:1582)
39: ==4763==    by 0x3EBB6B: std::vector<double, std::allocator<double> >::vector(std::initializer_list<double>, std::allocator<double> const&) (stl_vector.h:629)
39: ==4763==    by 0x3721F9: WorldBuilder::Features::OceanicPlateModels::Temperature::PlateModel::get_temperature(WorldBuilder::Point<3> const&, WorldBuilder::Objects::NaturalCoordinate const&, double, double, double, double, double) const (plate_model.cc:168)
39: ==4763==    by 0x35DF6D: WorldBuilder::Features::OceanicPlate::properties(WorldBuilder::Point<3> const&, WorldBuilder::Objects::NaturalCoordinate const&, double, std::vector<std::array<unsigned int, 3ul>, std::allocator<std::array<unsigned int, 3ul> > > const&, double, std::vector<unsigned long, std::allocator<unsigned long> > const&, std::vector<double, std::allocator<double> >&) const (oceanic_plate.cc:176)
39: ==4763==    by 0x3CE692: WorldBuilder::World::properties(std::array<double, 3ul> const&, double, std::vector<std::array<unsigned int, 3ul>, std::allocator<std::array<unsigned int, 3ul> > > const&) const (world.cc:387)
39: ==4763==    by 0x31450E: main (main.cc:305)
39: ==4763== 
39: 1500e3 1000e3 0 10e3 500.274  0
...

[lhy11009]

@danieldouglas92 It's nice that you have these implemented and had a discussion with me to improve it.

Previous discussion:

1. Agreement on the bulk part of the implementation
2. Question about the time shifts: it's not quantitatively tested
3. Issue that the time shifts are added twice for the effective_plate_age

my WIP:

1. New function for calculating the effective trench and plate ages
2. unit tests for this function and the function of calculate_ridge_distance_and_spreading.

When I am done, I'll create a PR for your branch.

I have some additional questions about the calculate_ridge_distance_and_spreading:

1. It seems that in the last for loop over the segments of the relevant ridge, the implementation determines the nearest one, but the spreading velocities are not updated accordingly. I think if you follow the variable "spreading_velocity_at_ridge", you'll get the value at the last segment from the for loop.
2. the ridge_migration_time: what is the start point and the end point of this time duration? Does it end at the current state?

It's your call if you want to reply here or arrange another 30 min meeting.

Just to save further emails:
My available time:
Friday: 11:30 am - 13:00 pm in Eastern time
Next Monday: 11:00 am - 12:30 pm
Zoom: 7989138001

[mibillen]

Hi all, I've read through the back and forth conversation here but I don't have time to review code right now, and I'd to wait until Haoyuan and Daniel are agreed that things are working. From the schematic and the demo figures it looks very promising. Key thing is to think of tests that can address "edge" cases. For example, will this work for any slab length (is there ever an issue where the slab it too long or too short?). Can you put in an assert that catches if the age at time of subduction becomes negative? This should never happen, and since this is not just an input variable now, someone could do something tricky with ridge location and spreading rate that ends up a-physical. Thank you both for working on this... its very cool feature to add.

[lhy11009]

@mibillen Just to keep you updated. I separated the calculation of slab ages into a separate function here:
danieldouglas92#1 (comment)
Some tests are under way for this.

[danieldouglas92]

I addressed most of the comments and added some asserts dealing with the input parameters (making sure that they are consistent with the dimensions specified in ridge coordinates, but yes I think that adding an assert for the age to make sure it is never negative would be a smart thing to do

[lhy11009]

Once you add the assert function, I'll create a unit test to check the throw.

[danieldouglas92]

Done! @lhy11009

[danieldouglas92]

@lhy11009 is there anything that I can be doing to help finish this up? I think it's just the unit_test for the age function you implemented that is left (unless I'm forgetting something else), and then I can clean up the PR and try to get this merged.

[lhy11009]

Hi Daniel Have you seen the new PR opened last night? It's mainly a test of the ridge distance function in spherical. Can you address the question I have there? Cheers

…

On Wed, Feb 28, 2024 at 8:02 AM Daniel Douglas ***@***.***> wrote: @lhy11009 <https://github.com/lhy11009> is there anything that I can be doing to help finish this up? I think it's just the unit_test for the age function you implemented that is left and then I can clean up the PR and try to get this merged. — Reply to this email directly, view it on GitHub <#654 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AG2WMFOSWIT2NZB5WLYVFMDYV5IJJAVCNFSM6AAAAABDNBOHMKVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTSNRZGMYDGOJSHA> . You are receiving this because you were mentioned.Message ID: ***@***.***>

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

[danieldouglas92]

I did not see that PR thank you for directing me to it!!

[lhy11009]

Hi Daniel

As you said in the previous conversation, this is ready for a merge. It's nice we can get this done together.

I'll have a discussion with Magali later about the calculate_effective_trench_and_plate_ages function and hopefully summarize it into another PR later.

@MFraters What do you think?

[danieldouglas92]

@lhy11009

I changed the argument in the function get_temperature_analytic from plate_velocity to subducting_velocity, and I just want to make sure that this is the correct velocity to be passing to this function, other than this I just need to squash the commits and address any comments from Magali and Menno

[lhy11009]

Yeah, this should be correct. At this point, as long as all tests are passed, it's ready to be merged.

[lhy11009]

Just another small comment.

Can you add a comment on the variable of "ridge_migration_time"? Like the definition of this variable.

The reason I ask is that we don't assign any time to ridge spreading, then why we need to assign a time to ridge migration?

[MFraters]

Cool, thanks for your work @danieldouglas92 and @lhy11009! It has been great to see you work together on this pull request :)

I will try to have a look at it today. In the mean time:

1. can you fix the typos found by the typos test (yay, the typos tester works as intended!). They show up as yellow warning boxes in the changed files tab.
2. Improve the coverage of the changes/additions in both source/world_builder/parameters.cc and source/world_builder/utilities.cc. This doesn't have to be a 100%, but there seems to be still some low hanging fruit for better coverage.

[lhy11009]

Okay, I see, I'll deal with the test coverage in utilities.cc. @danieldouglas92 Can you deal with the other file and the typos?

[danieldouglas92]

@lhy11009 I'm about to submit a change that addresses everything Menno said so I think it's all good! You don't have to worry about it but thanks

[MFraters]

here are already a few comments

[MFraters]

todo?

[MFraters]

Why did you change the default from 0.01 to 0.05?

[danieldouglas92]

I just thought that 1 cm/yr seemed to slow and 5 cm/yr made more sense!

[MFraters]

Could you add a comment somewhere what distance1 and distance2 represent, or make the name more clear? Secondly, it doesn't seem to be used in any of the tests, so could you add a comment on in what cases this is used, and/or add a testcase where it is used?

[lhy11009]

@danieldouglas92 nice job. Magali and I just touched base on what the effective age and the age_at_trench mean in the mass conserving model.

The effective age could be calculated by the total distance of the query point on the slab to the current ridge, while the age at the trench was the trench age when this query point is at the trench (purple Age_tr_1 in this case) and it could be derived by subtracting the time of subduction. (d_plane / Vsubduction) if the subduction rate doesn't change.

The idea of following the query point for it's current effective age and past trench age is to conserve energy while subducting back in time.

The figure shows a case with no ridge migration but with trench retreat motion (Vsubduction > Vspreading). So it's straightforward to compute these two.

Moreover, I think to move forward with ridge migration, the approach is to propagate that effect into the subducting velocity and create a function to calculate t_subduction (time of subduction of the query point).

The current implementation is good in that we match the previous test. But I'll soon make an update to address this in the following PR.

Daniel, can you go to the function of calculate_effective_trench_and_plate_ages and remove the two time shifts at this point? Such that we are dealing with only a trivial case at this point?

Let me know what you think. Should we touch base on this as well? I can meet Fri between 11:00 am - 14:00 pm, Eastern time

[danieldouglas92]

Yes we can call tomorrow at 11:00am EST! I think it would be good to discuss what you mean because as far as I can tell what you say is what is currently implemented? The effective age is being set based on the distance that a point has moved (d_plane in your diagram) relative to the ridge but is scaled based on the subducting velocity. I think the trench_age_shift should be removed yes but the effective_age_shift seems correct to me.

Also, trench_age_shift is currently multiplied by 0 anyway so it's effectively pointless having it in the PR anyway.

[danieldouglas92]

effective_age_shift = (1 - (distance_along_plane / ridge_drift_distance)) * sign_v * distance_along_plane / subducting_velocity

Just to clarify my previous comment, this is how we calculate the effective age shift, which is equal to what you call t_subducting in your diagram, and this is added to effective_age which is again what you are saying in your diagram, I think whether you call it trench migration or ridge migration, the effect is the same in that the age of the slab is changing relative to how the trench/ridge are moving relative to each OTHER, not relative to WHICH is moving

[lhy11009]

Hi Daniel, see you tomorrow at 11:00 EST !

we are yet to come to the same page with ridge migration and this sign_v multiplier.

I agree the distance of ridge migration would affect the total distance of subduction. But wouldn't that be handled in the subducting velocity? Can't we just use that value?

Also, are you preparing to return different values of ridge drift for different query points on the slab? Like, if the ridge starts to move at 0 Ma, and subduction also initiates, then at 10 Ma, a deeper point would have a larger ridge migration time than the shallower point.

See this formulation

effective_age_shift = (1 - (distance_along_plane / ridge_drift_distance)) * sign_v * distance_along_plane / subducting_velocity

A naive way to think is the ridge_drift distance must be different cause it divides the distance_along_plane so they should both vary with the position of the query point.

No worry, it wouldn't require much effort to address these questions and move forward.

And before, I haven't asked these questions. And I think these are not what we need to address anyway in this PR.

[danieldouglas92]

For your example of ridge migration and subduction both beginning at 0 ma, the current formulation does not increase the age more down dip, it does the opposite. In fact, the largest shift happens at the trench, and when distance_along_plane > ridge_drift_distance there is no age shift at all.

Anyway, it will be good to talk about this to tomorrow because i'm sure the problem is actually easy to figure out if we discuss it for a few minutes.

[danieldouglas92]

@lhy11009 I actually sat down and thought about what the implementation is actually doing and you're right that the implementation is definitely not doing what I thought it was doing! Looking forward to hearing your plan tomorrow

[MFraters]

Based on the discussion this looks good to me. @lhy11009, can you confirm that this what you and @danieldouglas92 discussed, and if so approve the pull request?

[MFraters]

please update this

[lhy11009]

looks good to me

[lhy11009]

nice work

[MFraters]

Thanks again for your work on this!