Lu et al. present a novel coupling framework within ALE-IB scheme, which couples the thermo-mechanical code Underworld 2 with the landscape evolution code Badlands. They suggest that this module could give an accurate material boundary and prevent numerical artifacts. The topic is timely and relevant to the geodynamic modeling community. I suggest publication after addressing the following comments.

Main concerns:

1) This new coupling module is based on the previous UWGeodynamics module, which is also coupled to the geodynamic code Underworld 2 with surface process codes such as Badlands. What are the detailed differences between your coupling module and the UWGeodynamics module? I speculate that you added the IB scheme, but the description in Section 2.2.2 is not clear about how it is implemented. Some details about this IB scheme need to be added. How exactly is the velocity evaluated on the S-surface and transferred to Badlands? What interpolation methods are used? The particle adjustment process based on surface changes (distinguishing sediment from air) requires more detail.

2) The benchmark results in Fig. 4 are not very clear to me (see details below).

3) The sediment distribution differences between schemes (Figure 10a,b) are subtle, yet the tectonic deformation is more different at the collision front (Fig. 10c,d). Why?

Minor:

38, typo: FasScape -> FastScape

Eqs 4-5, the C and phi are messed up. It’s better to define the starting cohesion (C_0) and friction angle (phi_0).

131, do you have a criterion for choosing the subtimesteps dt_sp?

138-139, Is the grid used in Badlands the triangular grid as shown in Figure 1? Although the ALE-IB scheme avoids interpolation from particles to grid nodes at the surface in Underworld 2, how is the processing carried out in the process of converting from surface rectangular grids in Underworld 2 to more refined rectangular grids and further triangular grids in Badlands? If interpolation is still used, what impact will it have on the result?

140-144, As mentioned above, a more detailed ALE-IB implementation scheme needs to be described.

167, the air density is set to be 0 kg/m3. Do you solve the energy equation, and does it have problems?

208, add units for Kd & Kf

Fig. 4, SP should add the abbr. for SP-surface process in the figure caption.

Fig. 4, why isn’t the deposition/sedimentation at the trough? Did the two models use the same random noise? Does the change of random noise have a significant impact? Besides, a figure of the differences in topography and erosion/deposition between these two schemes can better show the distinctions between the two schemes.

Fig. 5, drainage colorbar adds unit

Table 2, your rheological model has no pressure dependence, activation volume is not used, but you show it here. So, do you use the activation volume in your viscosity?

Fig. 12 is not used in the main text; you may want to refer to it at line 248 (Fig. 12e, not 11e?)

Table 2 defines the model names for two groups of experiments: Ex.1 vs. Ex. 2. It’s very confusing when referring to e.g., CM1 from Ex.1 or Ex.2. Please use a better approach for clarity.

“A novel ALE scheme with the internal boundary for coupling tectonic and surface processes in geodynamic models” by Lu et al., two-way couples surface processes to tectonic models using two different implementations, Eulerian and ALE-IB. They compare the two schemes using topography relaxation and continental collision models and find that ALE-IB is more accurate than the Eulerian method. I’ve often heard about the different implementation schemes, and I am very happy to see a comparison done! I think this comparison and, as mentioned by the authors, a future comparison to a pure ALE scheme are important contributions to the field. However, section 2.2 would benefit from a more detailed description of the two methods and how their implementation differs. Also, I think the manuscript would benefit from more quantitative discussion about the differences between the two, and what circumstances it may be necessary to consider one over the other.

Main comments:

The authors should expand section 2.2 to give a more in-depth description of the two implementations, as well as how they differ. In both descriptions, it sounds like the velocity is interpolated from the tectonic model to the surface processes model, and then badlands does the advection and SP routines. Is the difference only in how the interpolation is done or is there some other major difference between the two? I think updating figure 2 to show a flow chart of the Eulerian scheme next to the ALE-IB may help to understand the differences.

It would be helpful to give some discrete numbers on how much the two methods differ from each other, and in the topography relaxation model, from the analytical solution. For example, in figure 4ab when examining closely the ALE-IB does seem more accurate than the Eulerian scheme, but the Eulerian scheme already seems quite accurate. How much accuracy is gained? A similar quantitative comparison would strengthen the discussion of Figures 4e–f and Figure 10. The visual differences appear minor in some cases; quantifying them would clarify whether they are negligible or meaningful for certain applications.

I do not think additional models are entirely required for the paper; however, it would be very helpful to include one case with a high kf value using the Eulerian scheme. Sections 3.1.2 and 3.2.2 briefly discuss increasing kf, but only for the ALE-IB implementation. Does the advantage of ALE-IB increase as erosion and deposition intensify, or do both schemes perform similarly under these conditions? Addressing this would clarify when the ALE-IB approach is necessary.

Minor suggestions:

Line 141 – Should this be the that velocity is evaluated on the T-surface and then transferred to Badlands?

Line 210 – Could the authors put numbers on how close the Eulerian and ALE-IB agree with the analytical solution?

Line 212 and 244 - An arrow on figure 4 and 11 to indicate what the authors mean would be helpful. 

Does the accuracy between the two schemes change depending on the number of particles per cell or model resolution?

Is there a difference in model runtime when switching to ALE-IB, or is it getting more accurate results with no significant change in computation cost?

From these results it appears while having surface processes affects tectonics, changing between the two schemes does not have a large change. However, because of the change in deposition, especially in the continental collision model, the ALE-IB accuracy may be particularly important when comparing modelled basins with natural observations. In the authors opinion, when is it necessary to use ALE-IB over a Eulerian scheme?

Figure 2: Could the authors add a flowchart of the Eulerian scheme to highlight the difference?

Figure 4: the SP acronym should be defined within the figure. Sealevel should be indicated on a/b to make it clear. subplots e and f need more description, what does the black line and gray area explicitly show? I think this should be made clear for all the figures that plot this.

Figure 10: I notice in the viscosity plots that the ALE-IB has a smooth changes in slope at the surface (10c) compared to the more cell-defined one seen in 10d. Is this one of the changes from the new scheme?

Table 1: Are Ex. 1 and Ex. 2 supposed to refer to whether it is the time relaxation or continental collision models? I think the authors should describe this in the caption, and add something to the model name to distinguish between which set of models they belong to.