Isostatic initial topography

[lhy11009]

Pull Request Checklist. Please read and check each box with an X. Delete any part not applicable. Ask on the forum if you need help with any step.

This adds an "Isostatic topography" module to the "Inital topography model".

Before your first pull request:

• I have read the guidelines in our CONTRIBUTING.md document.

For all pull requests:

• I have followed the instructions for indenting my code.

For new features/models or changes of existing features:

• I have tested my new feature locally to ensure it is correct.
• I have created a testcase for the new feature/benchmark in the tests/ directory.
• I have added a changelog entry in the doc/modules/changes directory that will inform other users of my change.

[lhy11009]

A model set up with a free surface on top typically runs with initial steps of isostatic adjustment, where the velocity could be very big. This is the discussion we have in #7003.
In the PR, I use the model initial conditions (e.g., temperature, composition) to approximate the isostatic topography as the initial topography.

[lhy11009]

This is an example from the tests I included, with a continent on the left side and an ocean on the right side. The blue curve results from running the model with no initial topography, and plotting the equilibrated profile at t = 2e5 yr. The orange curve is the result from the isostatic initial topography (Units: X (km), topography (m)).

[lhy11009]

Another example is from the subduction-collision model I am running. This model has a Continent-Ocean-Slab-Continent setup. There are low-wavelength features in the calculated isostatic topography, near the slab trench (x = 3000 km). Actually, the previous example also shows an abrupt change in a step-function, a low-pass filter might be needed. After, isostacy works for long-wavelength topography.

[lhy11009]

The implementation of this current has several tricky points, which are marked as the following:

[lhy11009]

This "additional" initialize step is needed, as the "Isostatic topography" needs to query the initial composition, temperature, as well as evaluate the material model after these are all initialized. These information are used to derive column-wise density below every surface point to determine how isostacy could be attained.

[lhy11009]

During this "additional" initiation state, sample points are distributed on the surface where the column-wise isostacy is computed.

[lhy11009]

This "additional" initiation is added here in the core.cc, after other related modules.

[lhy11009]

I am current unsure what is the best way to work through the interfaces. Also the current implementation only works in 2-d box geometry, so following PR would potentially handle the other geometries.

[lhy11009]

@anne-glerum, would you take a look and tell me what you think.

[Djneu]

Hi @lhy11009 and @anne-glerum , we've been meaning to get to this for awhile, and I think it'll be a huge quality of life improvement, so it is great to see you working on it! I haven't had a thorough look through the code yet, but here are just a few questions/ideas I was having on tackling this problem.

1. When you add initial topography does it add mantle material or does it only thicken the material that is already at the surface? (Anne mentioned this may depend on whether this is added using an initial topography or through mesh deformation).

2. Could you modify this where the reference isostasy is calculated from a user-defined point? If lithostatic pressure boundaries are used, the way the isostatic rebound occurs should depend on where that reference point is (I think, or does it occur according to something else?), so the isostatic profile would need to also consider the same point.

3. I am not sure about the ordering of things when starting up ASPECT, but has the adiabatic pressure been calculated by the time required_initialize is run, and if so is it necessary to calculate the mass or could we only look at pressure differences at a reference compensation depth?

4. How does the isostatic profile work with increasing x-direction resolution? I wonder if it is necessary to find the topography at all surface points (though I could see this perhaps being a bit difficult with changing resolution with depth), or if a lower res topography can be established.

5. In relation to point 4, in current_surface.cc, we use the surface points to create a function to describe the surface topography using:

surface_function = std::make_unique<Functions::InterpolatedTensorProductGridData>(coordinates, topography_data_table);

I wonder if, to get around possible resolution issues, it might be easiest to find the topography points (or pressure profile compensation depth) and create a function like this so any possible surface point can be queried using something similar to:

double topography = surface_function->value(p);

[lhy11009]

@Djneu For your questions 1, 2, 3. I'll try to write more about how this works.

[lhy11009]

Here is my way of thinking about this once we have the initial inputs of composition and temperature like the subfigure on the left, basically continent and oceanic plates side to side with no topography in a box geometry.

Two pieces of information are extracted from this state, the integrated mass of the column below each surface point, and the total mass in the domain.

The idea is the total mass should be conserved in the isostatic adjustment, and the column-wise integrated mass should be equal below different surface points, as shown by the subfigure on the right. To achieve this, modify the topography by migrating some mass from some columns to others (blue arrow: from below the ocean to below the continent). In detail, this is achieved at the balanced depth.

Therefore, the initial topography at one point is computed by taking the deficit of the mass in the column relative to the expected average, and dividing it by the reference density at the balanced mass (h(x)).

This addresses parts of question 1,2,3, in terms of what information is needed. In addition, a reference point is not needed for the isostacy computation, as that is basically achieved through a lateral mass average of the whole domain.

[lhy11009]

Then, to your question 1. I also try to think about the differences between whether this is added using an initial topography or through mesh deformation. I don't really have a clear idea about how mesh deformation works. But I understand there will be a composition change due to adding initial topography.

Having a positive topography creates new space on top, and this won't be consistent, except if your composition layer extends to negative depth. On the other hand, having negative initial topography is effectively cropping out the compositional layer. I learned these from playing with the world builder files with ASPECT's initial topography. This is indeed different from deforming the mesh and moving things upward and downward while preserving its original layering.

It would be helpful if three of us could talk through possible alternatives to doing this in mesh deformation. My naive thought is that it might be feasible to deform the mesh once in the initiation stage.

The workaround considered but not implemented yet, is to have the initial composition module initiate a second time after the "required initiation" of the isostatic topography finished. This time, update the positions of compositions using the depth in the model

[lhy11009]

To your question 4.

The blue and orange are the same profiles as shown previously, while the green is a model run with lower resolution (refinement 4 vs 6, resolution is about 60 km vs 15 km), and its equilibrated topography at 2e5 yr.

The initial topography by the orange profile is independent of whether a coarser or finer mesh is picked by the model. The idea is that the inputs are not from the mesh points, but directly from modules they initiated from (e.g. initial temperature, composition, adiabatic conditions).

A little detail of the implementation is that there is no drawback in using very high res to initiate the topography by, for example, 100s meter resolution, because it doesn't loop the mesh points. This is only needs to be initiated once and is kept by the "topograph" vector in the plugin. Then, when values are queried at each query surface point later, the algorithm ensures O(1) cost.

And, maybe this is an odd point to make, but a model with low resolution would be problematic, as the equilibrated topography would deviate from higher-res models and could not be by captured by the initial isostatic model.

[lhy11009]

To your question 5

Yes, we might replace the current "topography" vector in the plugin with this data structure. This is also something we can discuss.

[lhy11009]

@Djneu and @anne-glerum , my location is in California. It would be great if we could 7 or 8 pm your time (9 or 10 am PST). Potentially tomorrow? We could talk through questions we had, and the planned progress on this PR

[Djneu]

Sure, how about 8 pm our time?

[lhy11009]

Yes, I'll send you a Zoom invitation soon.

…

On Wed, Jul 22, 2026 at 6:44 AM Derek Neuharth ***@***.***> wrote: *Djneu* left a comment (geodynamics/aspect#7035) <#7035 (comment)> Sure, how about 8 pm our time? — Reply to this email directly, view it on GitHub <#7035?email_source=notifications&email_token=AG2WMFOSE4IH3E3RWZACKOT5GDAMPA5CNFSNUABFM5UWIORPF5TWS5BNNB2WEL2JONZXKZKDN5WW2ZLOOQXTKMBUGY3DQNBWGE3KM4TFMFZW63VHNVSW45DJN5XKKZLWMVXHJLDGN5XXIZLSL5RWY2LDNM#issuecomment-5046684616>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AG2WMFKPIBU5VWRQKLKR46L5GDAMPAVCNFSNUABEKJSXA33TNF2G64TZHMYTSNJQGI2DKMR3JFZXG5LFHM2DSMZSGQZTONBYGGQXMAQ> . Triage notifications, keep track of coding agent tasks and review pull requests on the go with GitHub Mobile for iOS <https://github.com/notifications/mobile/ios/AG2WMFKRTJHBK53BPXDDWD35GDAMPA5CNFSNUABFM5UWIORPF5TWS5BNNB2WEL2JONZXKZKDN5WW2ZLOOQXTKMBUGY3DQNBWGE3KM4TFMFZW63VHNVSW45DJN5XKKZLWMVXHJKTGN5XXIZLSL5UW64Y> and Android <https://github.com/notifications/mobile/android/AG2WMFNLHTMR3TEYUCHSXOL5GDAMPA5CNFSNUABFM5UWIORPF5TWS5BNNB2WEL2JONZXKZKDN5WW2ZLOOQXTKMBUGY3DQNBWGE3KM4TFMFZW63VHNVSW45DJN5XKKZLWMVXHJLTGN5XXIZLSL5QW4ZDSN5UWI>. Download it today! You are receiving this because you were mentioned.Message ID: ***@***.***>

-- Haoyuan Li Postdoctoral Researcher Earth and Planetary Sciences, UC Davis

[lhy11009]

Haoyuan Li is inviting you to a scheduled Zoom meeting. Join Zoom Meeting https://ucdavis.zoom.us/j/96020571647 <https://www.google.com/url?q=https://ucdavis.zoom.us/j/96020571647&sa=D&source=calendar&ust=1785167845110402&usg=AOvVaw2pH-2hH0-uIRI81TFfkXSp> Meeting ID: 960 2057 1647

…

On Wed, Jul 22, 2026 at 8:53 AM Haoyuan Li ***@***.***> wrote: Yes, I'll send you a Zoom invitation soon. On Wed, Jul 22, 2026 at 6:44 AM Derek Neuharth ***@***.***> wrote: > *Djneu* left a comment (geodynamics/aspect#7035) > <#7035 (comment)> > > Sure, how about 8 pm our time? > > — > Reply to this email directly, view it on GitHub > <#7035?email_source=notifications&email_token=AG2WMFOSE4IH3E3RWZACKOT5GDAMPA5CNFSNUABFM5UWIORPF5TWS5BNNB2WEL2JONZXKZKDN5WW2ZLOOQXTKMBUGY3DQNBWGE3KM4TFMFZW63VHNVSW45DJN5XKKZLWMVXHJLDGN5XXIZLSL5RWY2LDNM#issuecomment-5046684616>, > or unsubscribe > <https://github.com/notifications/unsubscribe-auth/AG2WMFKPIBU5VWRQKLKR46L5GDAMPAVCNFSNUABEKJSXA33TNF2G64TZHMYTSNJQGI2DKMR3JFZXG5LFHM2DSMZSGQZTONBYGGQXMAQ> > . > Triage notifications, keep track of coding agent tasks and review pull > requests on the go with GitHub Mobile for iOS > <https://github.com/notifications/mobile/ios/AG2WMFKRTJHBK53BPXDDWD35GDAMPA5CNFSNUABFM5UWIORPF5TWS5BNNB2WEL2JONZXKZKDN5WW2ZLOOQXTKMBUGY3DQNBWGE3KM4TFMFZW63VHNVSW45DJN5XKKZLWMVXHJKTGN5XXIZLSL5UW64Y> > and Android > <https://github.com/notifications/mobile/android/AG2WMFNLHTMR3TEYUCHSXOL5GDAMPA5CNFSNUABFM5UWIORPF5TWS5BNNB2WEL2JONZXKZKDN5WW2ZLOOQXTKMBUGY3DQNBWGE3KM4TFMFZW63VHNVSW45DJN5XKKZLWMVXHJLTGN5XXIZLSL5QW4ZDSN5UWI>. > Download it today! > You are receiving this because you were mentioned.Message ID: > ***@***.***> > -- Haoyuan Li Postdoctoral Researcher Earth and Planetary Sciences, UC Davis

-- Haoyuan Li Postdoctoral Researcher Earth and Planetary Sciences, UC Davis

[lhy11009]

@Djneu Just, In case my previous message doesn't reach you.

https://ucdavis.zoom.us/j/96020571647

[lhy11009]

I have shifted things to mesh deformation, and implementation now works. This interface is indeed better, and there is no need to change core.cc anymore. The resulting initial topography (orange) also matches the evolution of an ASPECT model without an initial topography (blue).

[lhy11009]

Both this and the previous approaches have problems with changing the initial compositions. For example, here I show the particles with the test case with the updated approach. The initial composition of oceanic crust (figure on top) and continental crust (figure on the bottom) both follow the original 0 depth and don't align with the adjusted surface. Note that the gray dots are artifacts of Paraview.

[lhy11009]

This result is a little surprising. It seems to me that mesh deformation should assign constraints to the DOFs in the solving steps and therefore should deform the original layers and preserve their thickness.

@anne-glerum , @Djneu. Anyway, let me know what you think. And should we first stick to what we have here?

[anne-glerum]

Do you see the same effect when using fields instead of particles to track composition? I'll have a look at your PR tomorrow (today we have a half day off).

[lhy11009]

Do you see the same effect when using fields instead of particles to track composition? I'll have a look at your PR tomorrow (today we have a half day off).

Oh yes, the results are the same; both of these are not moving with the deformation mesh. I wonder if you checked something like this for the ASCII-file plugin? Or if someone has tested these with the initial constraints on mesh deformation in generl.

[Djneu]

Hi all, @lhy11009 @anne-glerum, I'm tagging @Minerallo and @MFraters here as well since they are going to look at mesh deformation/initial topography with worldbuilder, and I think they may run into a similar issue as the non-uplifting layers here (although perhaps since worldbuilder uplifts the composition it won't be an issue depending on ordering). Should we all maybe have a discussion on some point about this?

[lhy11009]

@anne-glerum @Djneu. This morning I had a discussion with Rene initially on another matter. Then I raised this question and he has a good point.

Currently, material by default doesn't move with mesh deformation. It's designed so as to couple the process in Fastscape. As with erosion or sedimentation, material is either added in or taken out, and the material in the mesh remains in its original places, rather than moving with the mesh. Moreover, this is currently done by changing the advection function with a compensation term that subtracts the mesh deformation velocity from the real velocity

In other words, the current implementation of mesh deformation might implicitly assume mass flow through the boundary rather than assume material movement with the mesh.

I'll do a test with a boundary composition option to check whether this is the case.

In the new test, I added an option for allow fixed composition on outflow boundary:

		subsection Boundary composition model
		  set Allow fixed composition on outflow boundaries = true
		  set Fixed composition boundary indicators = top
		  set List of model names = function
		  subsection Function
		    set Variable names = x,y
		    set Function constants = W=5500e3
		    set Function expression = if(x<W/2,1,0); if(x>=W/2,1,0)
		  end
               end

For both the field and the particle method, the composition field of the continent crust is fixed at the surface. However, the bottom limit of this composition remains unmoved.

(field method)

(particle method)

This test demonstrates that the issue of not preserving the composition at the top surface arises because material is not advected with mesh deformation.

[lhy11009]

Based on this understanding, I think we could proceed with this PR and handle the issue of composition elsewhere. This issue wouldn't change the calculated initial mesh deformation, which is the sole purpose here.

Ideally, we could make some progress with this PR to have it merged in the Hack. Then we could include a low-pass filter and the reference point in a following PR.

I could discuss with Rene again and see how we could solve the issue of composition elsewhere.

@Djneu @anne-glerum, Ideally, we could touch base on this. Would 8 pm your time still be a good option?

[Djneu]

@lhy11009, I agree with that plan and thanks again for implementing this! This implementation should work in most cases, and should be close enough to significantly help things even if the compositions are thickened instead of uplifted. For the reference point, did you test a model using a traction boundary? Perhaps it won't even be an issue, but I agree if it still is needed for those cases we can always add it later.

How much has this implementation helped with the initial timesteps velocity and runtime for models that need to rebound?

[lhy11009]

I have tested the model with initial topography myself, and that would help a lot with the initial rebound stage. velocity-wise, this appears to lower the velocity magnitude by an order.

I agree if the boundary condition changes, the situation would be different. In that case, this method would still help to bring things close to isostacy

In terms of mesh deformation not preserving the composition, we can talk with Rene to think of changing how advection is currently handled. Or add an option to just change it for the initial step.

For the current PR, there is still a hard-coded feature. For example, it now only works for a 2-D box model. Feel free to make some suggestions on what to extend before the PR is merged.

[lhy11009]

I do think I need to add a filter to this and an isostatic length scale. In this way, sharp topography would be filtered out. This seems to be needed to make it useful in application.

I have added the Gaussian filter dictated by the new variable of "Isostatic length scale". The natural meaning is that the low-wavelength topography is supported by the elastic strength of the lithosphere. While long-wavelength topography is expected to achieve isostacy (Airy?)

The green curve results from adding the filter and a 500-km isostatic length scale. It almost exactly matches model evolution from no initial topography at 2e5 yr (blue curve). In comparison, the case without this filter is shown with the orange curve.

This new case is also added as a test.

I have also added some comments to the variables and tests. This is now satisfactory for me. Would you like to give a review?

[lhy11009]

I did some study on mesh deformation implementation. I think that moving the material with the initial mesh deformation is not easy to achieve in the current structure.

@MFraters, I'd like to know whether this will be an issue with GWB's initial topography. When you have time, feell free to give us some feedback on we have tested here too.

I can see two route paths:

1. As in Arbitrary-Lagrangian-Eulerian (ALE) mesh and material velocity are separate matters. Following this thought, one has to update material locations as a separate step after initial mesh deformation. It also has to assume a Lagrangian framework at this step, or to say, make mesh and material velocity equal to each other.

2. Move the Initial mesh deformation to an earlier initiation workflow, or wrap the compute_initial_deformation_on_boundary function to earlier initiation workflow. Then have initial composition and initial temperature reinit after it.

It would be good to discuss these possibilities.

[Minerallo]

I did some study on mesh deformation implementation. I think that moving the material with the initial mesh deformation is not easy to achieve in the current structure.

@MFraters, I'd like to know whether this will be an issue with GWB's initial topography. When you have time, feell free to give us some feedback on we have tested here too.

I can see two route paths:

1. As in Arbitrary-Lagrangian-Eulerian (ALE) mesh and material velocity are separate matters. Following this thought, one has to update material locations as a separate step after initial mesh deformation. It also has to assume a Lagrangian framework at this step, or to say, make mesh and material velocity equal to each other.
2. Move the Initial mesh deformation to an earlier initiation workflow, or wrap the compute_initial_deformation_on_boundary function to earlier initiation workflow. Then have initial composition and initial temperature reinit after it.

It would be good to discuss these possibilities.

Hi @lhy11009, nice work! Having isostatic initial topography would be very useful. I have been calculating it externally for a while, but that is not an ideal solution.
However, I see that you encountered some issues when initializing the compositional and temperature fields and that you are therefore forced either to advect the material after applying the initial deformation or to prescribe some form of inflow condition at the bottom. For the same reasons, I had chosen to go for the second option from the beginning In #7155 GWB initial topography. Although I understand that the current trend is to handle more of this through mesh deformation, I think that approach becomes quite complicated in this case.
In #7155, I added the initial topography from GWB as part of the geometry model. This means that the deformed geometry is constructed earlier and can then be initialized directly with the initial compositional and temperature fields. Implementing this through the mesh-deformation would likely require considerably more extensive changes, but if you feel like you want to have it working from the mesh deformation like option 1 then fine for me, but it will be more complicate to get it right.

[Djneu]

I am not sure whether the same method would work here, as in world builder you already know the amount of uplift and have applied it to the composition and temperature as an input, where as here I assume people would write it assuming a flat surface not knowing how much it would uplift. Unless we first run this script to find the deformation and then update the position of all parameters with this topographic change as an input?

[lhy11009]

Hi Michael, I have looked through your PR, and I think it's a great structure to have in GWB utilities. Yes, there are slightly different ideas between the two approaches.

1. Using the initial mesh deformation, the initial temperature and composition do not know how the mesh is going to deform. The mesh deformation is initiated last, and the information only gets into the solver .
2. Using GWB, initial topography could be set consistently with initial temperature and composition. The initiation of initial topography also comes first in the simulator; therefore temperature and composition can be set accordingly.

And the approaches I suggest are indeed not easy to achieve. If it's up to me, advecting the ALE one step would be my choice, because:

1. that aligns with the concept of isostatic adjustment.
2. One takeaway from these tests is that the ASPECT topography is going to adjust that profile anyway in a 1e5 yr timescale. So it's reasonable to have a single step addressing that rather than having this initial stage if that's not the focus.

With GWB, I think the way one could pre-compute the isostatic topography, perhaps like you did, is:
a. Use the GWB in ASPECT to initiate a model without initial topography; this way the EOS are all consistent with P, T conditions to compute the density.
b. Write a post-processing step to extract the density structure, compute the isostatic topography, and put it back in GWB
I have something similar that functions like that.

There is another drawback of this format of isostacy that I just realized, which are the features that don't adjust to long-wavelength isostacy in nature, basically subduction zones and LLSVPs. For these, it's very likely that an initial isostatic topography might impact things, and perhaps in a bad way. I have noticed that with a near-isostatic topography given to the slab as initial topography, the trench develops a much deeper depth in later model runs than models initiated with no topographies.

For the reason mentioned above, I think after the PR here is merged, what we can do is:

I. In the GWB workflow, the post-processing for isostacy is not needed; one can at least take the topography from the isostacy calculation and put it into GWB.
II. Use this to keep sorting out a better way to do initial mesh deformation.

[lhy11009]

A second thought. It just came to me that advecting the material with the mesh won't work because the mesh doesn't move in a physical way.

This actually raises a question for me: would initial topography be a better interface than the initial mesh deformation? I think Michael's approach of using initial topography consistent with GWB could indeed be a better solution.

[lhy11009]

I have another idea, basically the following workflow, similar to what I have for the PR originally using the initial topography:

intial_topography->initialize()
initial_temperature->initialize()
initial_composition->initialize()
material_model->initialize()

if intial_topography->requires_EOS():
    initial_topography->initialize_with_EOS()
    initial_temperature->initialize()
    initial_composition->initialize()

Here, initial_topography->initialize_with_EOS() would take the density computed from the material model and update the isostatic topography, similar to what I currently have for this PR. Then the initial_temperature and initial_composition could initialize for the second time, similar to Michael's workflow with the GWB, in the sense that both are informed with the updated initial topography.

Okay, let's pause here a bit. Cause I really want this to be a ready solution. We'll need more inputs from everyone included in the discussion.

For this purpose, I would like to arrange for a meeting on this before we wrap up the Hack. @Djneu, can you help organize this for me in your evening time? Perhaps get those who are interested, and if possible, also get Rene.

[Djneu]

Hi @lhy11009, sorry for the late reply, I had left the hackathon on the 29th and wasn't available to set up a meeting! For whether this should be in mesh deformation or geometry, I am not sure how much a difference it makes and I think we should discuss it with @anne-glerum and @gassmoeller.

I agree that regardless of which is used, we would need to reinitialize temp/comp after we know the surface, and we would need to somehow feed in the new topography to each query point so the fields don't end up in the same position. I think this is the key thing we are currently missing that GWB does by default when it adds in this topography.

For doing this I think we would need to know the surface after deformation. There is a postprocessor (current_surface) that has this functionality for 2D boxes, and @Minerallo has an open pull request (#7150) for extending this functionality into the mesh deformation itself. Though perhaps neither of those are needed, since you already have the change in topographies maybe you could do like current_surface.cc, and create a public function that can be queried for any topography change given an X-Y point? Though this may then limit it to only working with initial isostasy.

If that works, then I think we would just need to, if mesh deformation/initial topography is active, have initial comp/temp call this function and add the topography to the depth component for every queried point when it is initialized the second time. But maybe it is best we still all set up a meeting to try to discuss the way to do this since it would require some changes to core.cc!

[lhy11009]

@Djneu @anne-glerum A slight change of plan and a fix to the previous issue after discussion with @tjhei and @naliboff this Monday:

1. A change of plan to switch back to using the Initial topography system rather than the mesh deformation system. The reasons are

• A. The initial topography is currently included in systems of initial temperature and initial composition systems. They can get information as "depth()" or "maximal_depth()" from the geometry model.
• B. The initial topography is also included in the mesh deformation model in some thoughtful way, which seems that I shouldn't think of introducing mesh deformation to other systems like the initial temperature and initial compositiion.

2. The rest of the solution lies in using depth instead of y in the initial composition function. However, this has to depend on another PR (cartesian_to_depth_coordinates utilities in the geometry model #7292) I just created.

The idea is that the isostatic topography will update the initial topography within the geometry model, and this information could be later used by the initial composition as long as it uses depth instead of y to assign compositions.

Now, in the test, the composition is migrated with the adjustment of the topography.

[lhy11009]

This is a small leftover. After the isostatic adjustment, the maximal depth changes from the value the adiabatic condition initiates and could lead to issues below region with positive topography and near the bottom boundary.

[lhy11009]

To start the review, I think we need to first agree upon the changes to the core.cc and the inital topography interface. This change aims to allow the initial topography to initialize again after other modules (e.g. initial temperature, composition, material model), so that we can use these modules to compute the isostatic topography.

[lhy11009]

@Minerallo I think this structure is also consistent with GWB's topography in that one could use both the GWB topography and the isostatic topography. Conceptually, it will still call the isostatic topography to adjust the input from GWB, with all compositional layers preserved.

[lhy11009]

@Minerallo, actually, I just tested this with GWB. Did you already modify something like:

      return world_builder->composition(Utilities::convert_point_to_array(position),
                                         -this->get_geometry_model().height_above_reference_surface(position),
                                         n_comp);

to

        return world_builder->composition(Utilities::convert_point_to_array(position),
                                          this->get_geometry_model().depth(position),
                                          n_comp);

?

I haven't rebased to the main branch, but after I modify these lines in my branch, the PR here becomes very consistent with GWB. So when we have GWB prescribe the structures and composition layers, we can still have the isostacy to adjust the topography and preserve it's structures.

[Minerallo]

@Minerallo, actually, I just tested this with GWB. Did you already modify something like:

      return world_builder->composition(Utilities::convert_point_to_array(position),
                                         -this->get_geometry_model().height_above_reference_surface(position),
                                         n_comp);

to

        return world_builder->composition(Utilities::convert_point_to_array(position),
                                          this->get_geometry_model().depth(position),
                                          n_comp);

?

I haven't rebased to the main branch, but after I modify these lines in my branch, the PR here becomes very consistent with GWB. So when we have GWB prescribe the structures and composition layers, we can still have the isostacy to adjust the topography and preserve it's structures.

Hi Haoyuan, sorry for the late reply I was in holidays in the Alpes for the last weeks ! :)
Thanks for testing this. You are right: height_above_reference_surface() still refers to the original undeformed surface, whereas depth() includes the initial topography, so that would help to conserve the right compositional layer thicknesses.
My PR #7150 is related, but it calculates depth below the evolving surface during mesh deformation. For the initialization here, geometry_model.depth(position) should be sufficient. cheers. MP