General comments

The article “The effect of sediments on the dynamics and accretionary style of subduction margins” presents 2-D compositional models of free subduction with a layering of the subducting plate including a sediment cover (of variable thickness and viscosity) over a stronger mantle layer. The authors calculated model diagnostics during the simulation to categorize the subduction simulations either as accretionary margins or as erosive (non-accumulative) margins.

I find the model results and comparison to natural subduction zones relevant to discuss how the interface properties feedback with convergence rate through plate viscous coupling, despite a concern regarding the simulations featuring very-low viscosity sediments (see point 5 below).

I suggest below major corrections regarding

(i) the “angle of attack” of the paper (points 1 and 2 below)

(ii) choice and representation of the diagnostic parameters (point 3)

(iii) additional reference to previous studies and interpretation of some results (points 4-5)

(iv) a new repartition of figures between main article and supplementary (point 6)

These corrections do not require any new simulation to be run.

Major comments

1- The title and abstract do not at present capture what in my opinion are the main results of the models. Indeed, only late in discussion (line 440) do the authors acknowledge that “it is unclear whether sediment influxes affect convergence rate, or convergence rate affect sediment accumulation at trench, or both”.

I suggest that this should be very clear from right from the introduction that the sediment layer properties will be both an input in the models (viscosity, thickness) and also a diagnostic output (width, angle, volume) because feedbacks are expected between interface properties, viscous coupling, convergence rate and sediment accumulation.

Hence the paper could not only bring insights on “how sediment fluxes influence subduction dynamics and plate coupling” but also on “how the margins form in response to subduction dynamics” and how these feedbacks lead to unstable, time-dependent subduction dynamics and sediment fate (see point 2).

The abstract should also contain more explanations on the processes, since it is at the moment too descriptive, listing the features of the different margins.

2-Beyond the end-member classification between accretionary vs. erosive margins, I think that the model results could be used to highlight the temporal evolution of sediment fate in subduction zones. The gigantic error bars in Figure 4 (is a “mean” value relevant for e.g. radius of curvature or wedge width?) suggest the need for a more detailed post-processing of some simulations, to go beyond the mild sentence “In any given system, both processes may be occurring simultaneously, either in time and space or at the same time in different parts of the subduction zone” (lines 447-448).

Therefore I strongly suggest that the “unstable” mode is given more attention. This mode, briefly described in lines 303-307, suggest that “the accretionary wedge quickly reaches a maximum size and critical angle and instead of moving laterally the wedge, material is being expelled down the subduction channel.” The time-dependency of the sediment wedge geometry is in my opinion key to analyze for two reasons:

- it could explain why the models predict erosive margins at low convergence rates (Figure 4) whereas natural statistics show the opposite (Figure S2).

- the variation of convergence velocities through time because of “avalanche” effect of sediment wedge (or because of temporal variation in the sediment supply rate due to other factors) offers a new perspective on the spatio-temporal variations of plates velocities.

3- The 2-D volume of the accretionary wedge (triangle ABC in Figure 1c) should be an additional diagnostic parameter since it allows a clearer distinction between “tectonic erosion” vs. “accretionary margin” that the current criteria described in lines 287-290.

For example, I do not see on Figure 3 that only “low-angle accretionary margins have increasing radii of curvature” (high-mangle margins do too) or that low-angle margins have “fairly constant” whereas high-angles ones have “irregular” convergence rates (this is too qualitative).

Furthermore, tracking the evolution of this 2-D volume through time should highlight the temporal evolution of the sediment fate, for example by plotting:

- Figure 6 for unstable simulations

- co-evolution of 2-D volume (x-axis) and convergence rate (y-axis) as a function of time (marker color in the x-y diagram).

4- Please consider the following paper for the introduction dans discussion (first 3 are mandatory given the article topic):

- recent paper by Menant et al., Nature Comm, 2020 on the spatio-temporal evolution of forearc topography from transient stripping events of slab surface (basal accretion or erosion)

https://doi.org/10.1038/s41467-020-15580-7

For example, their discussion on topography evolution (see their Figure 3) is quite complete with a higher resolution than the present study.

- Cizkova & Bina, EPSL, 2019 - how plate interface viscosity plays a major role in controlling trench retreat https://doi.org/10.1016/j.epsl.2018.12.027

- Arcay, Solid Earth 2007 – influence of rheological properties of the weak layer on viscous coupling https://doi.org/10.5194/se-3-467-2012

- models of mass fluxes (and temporal evolution) of sediment fluxes in a subduction system – Beaumont et al., 1999 https://doi.org/10.1029/1999JB900136

- review on suduction erosion – Straub et al. 2020

https://www.nature.com/articles/s43017-020-0095-1

- analogue sandbox experiments on largins, and importance of convergence rate vs. sediment rate - Lallemand et al., 1994 https://doi.org/10.1029/94JB00124 or Gutscher et al., 1998 https://doi.org/10.1016/S0191-8141(97)00096-5

- strong interface coupling for a thin sediment layer and implications for seismogenesis:

Bangs et al., JGR, 2020 https://doi.org/10.1029/2020JB019861

- accretion at collision margins – Selzer et al., 2008 https://doi.org/10.1029/2007TC002169

- papers on how a weak crust may detach from the slab (in these papers in the transition zone, but perhaps similar to the “expulsion” of weak sediments down the subduction channel in unstable mode SubdSed04_100 on Figure S7), e.g. van Keken et al. GRL 1996 (https://doi.org/10.1029/96GL01594) and Yoshida et al. JGR 2012 (https://doi.org/10.1029/2011JB008989)

- papers discussing the order of magnitude of sediment viscosity since it appears critical an input in the models. For example, discussing the effect of sediment hydration/fluids on viscosity, the expected strength of detritic (continental) sediments vs. oceanic sediments.

The method section needs to be clearer on the choice of viscosity parameterization for the weak crust. At present, the sentence “The crustal thickness represents a parameterization of the strength weakening of the lithosphere with depth due to hydration and weak sediment cover” (line 155) is confusing and needs to be clarified: why is the depth-dependent parameter related to the sediment cover? Does hydration increased with depth?

The comment “high viscosity sediments are representative of a stronger crust and mantle component at the slab surface” on line 392 is too light.

5- I am worried that the sediment body labelled “high-angle accretionary wedge” (in Figure 2) is not a good analogue for the natural object. Indeed the sediments seem to accumulate only downward: does the sticky-air layer( viscosity > 10^18 Pa.s?) prevents the upward growth when the sediment viscosity is too low (10^18 Pa.s)?

The authors should comment in section 3 upon the thickness of the accretionary wedge, reaching 100 km in Figure 2 !!! The wedges simulated in Menant et al. 2020 (https://doi.org/10.1038/s41467-020-15580-7) are much thinner and appear closer to the observed geometries on Earth.

Thus, the “flat slab subduction” of the high-angle accretionary wedge result from the irrealistic accumulation of low-viscosity sediments downward, with the wedge forcing the slab to flatten at very shallow depths.

I therefore strongly suggest that the authors put the emphasis on “small-volume wedge” (analogue to erosive margins) vs. “large-volume wedge” (accretionary) rather than over-analyzing the differences between high- and low-angle accretionary margins.

6- Along with the change of focus of the article main results (see points 1-2-3 with suggestion of new figures), I suggest that:

- Figure S2 (and maybe S1) is moved to the main article to be directly compared with Figure 4. Beware to use the same scale and units for model (Figure4) and nature (Figure S2) graphics.

- addition of a figure (and movie?) showing an unstable simulation through time ( SubSed04_100?)

- Figure S7 is moved to the main article

- Figure 5 (topography) is moved to the supplementary

Intermediate comments

- the authors need to better describe the concept of tectonic erosion at non-accretionary margins, to discuss both the variation of volume of the arc crust and the removal of upper-plate material by basal erosion. The current figures do not show any erosion: on Figure 2, the velocity field is unreadable because too small. I suggest that the authors provide in the supplementary material a zoom of the interface/upper plate base to more clearly see whether ‘the sediment layer is an integrated part of the subducting slab and is eroding the upper plate.’ (line 275).

- Same thing about the viscous coupling that need to be better defined and better tracked in models outputs. Ex on lines 297-300: the localization of the decoupling motion between lower and upper plates is not visible from the too small velocity field vectors in Figure 2. The authors could either plot zoomed-in view of the velocity field at the interface, of plot a velocity profile along a transect through the interface to visualize the velocity transition.

- the authors need to clarify when they talk about correlations or when they refer to causal links (a process has been identified).

For example on lines 35-37: the very assertive sentence “On a regional scale, sediments influence patterns of deformation by controlling the morphologies of subduction interfaces, accretionary prisms and forearc basins” refers to three observational papiers and only one modelling paper (Simpson, 2010).

Example on lines 5-6 and 440: it is not clear whether the sentence refers to a descriptive correlation (“Accretionary margins are dominated by accretion of thick piles of sediments”) or to a causality (“tectonic erosion is favored in regions where the sedimentary cover is <1 km”).

- section 2: clarify whether “weak crust” (throughout the text) refers to both green layers on Fig.1 (“magmatic” crust and sediments) or only to the dark green layer (“magmatic” crust).

The sentence “15 km combined weak crust and sediments” is not clear: does the thickness of the magmatic crust vary opposite to the sediment thickness variation? (the sentence on lines 158-160 should appear earlier)

Did you run model with the same sediment thickness but different magmatic crust thickness to check that this did not have also an influence on the model dynamics? Indeed, when the thickness of the sediment layer is varied is the different simulations, the thickness of the magmatic weak crust does not remain constant.

- Figure 3: I see the point of the authors’ concern about variable timestep magnitude (section 3.1) however I am much concerned about plotting diagnostics parameters’ evolution as a function of (time-dependent) timestep and not of time or dimensionless time. It is very hard to understand subduction dynamics if the time-scale is dilated (in the current version, the spacing between different timestep may not be equal to the same absolute duration since we do not know the dispersion of timesteps around the mean value). Hence, I demand that Figure 3 plots show the x-axis with dimensionless time (e.g. time divided by the mean timestep shown in Fig. S3a).

Minor comments

-line 1: I would not refer to sediments as “continental material” which refers first to silicic magmatic rocks (e.g. granite). I suggest you keep the term “sediments”, or refer directly to the chemical elements.

-lines 2-3: if the “recycling of volatiles” links well with the subduction of sediments (full of water and carbon), the “petrogenesis of continental crust” and the recycling of “continental material” are not obvious to relate to the paper topic during the first read of the abstract > I suggest you rephrase into “… is important to the understanding of both recycling of water and carbon back into the mantle, and the petrogenesis of continental crust through water-induced partial melting.”

- line 4: perhaps rephrase the unclear sentence “When sediments are considered, convergent margins appear to fall into one of two classes”

-lines 6-7: what is the “geometry of the global subduction system”? It is not clear whether you refer so the accretionary wedge or to the slab.

- line 9: “how sediment fluxes influence subduction dynamics and plate coupling” could be more accurately rephrased as “how sediment fluxes influence subduction dynamics through plate coupling”

-line 11: you refer to “modes” of the subduction interface, however you only categorize them depending on the thickness of the sediment cover (erosion vs. accretion) and on the angle (high vs. low). Hence I suggest you refer to them as “geometries of the subduction interface”, that you can later relate to the system dynamics (how did these different geometries form over time?)

-line 13: I would be more precise regarding “the extent of viscous coupling” with “the lateral extent and the depth of viscous coupling”.

-line 14-15: you should clarify how the properties of the sediment layer modulate viscous coupling, ie. “When the viscous coupling is increased” or “when the viscous coupling is reduced” = how is the thickness/viscosity of the incoming sediment layer different? You need to make the link with the previous sentences categorizing the styles of interfaces (high viscosity, low viscosity).

- lines 16-17: the sentence “Diagnostic parameters are extracted automatically from numerical simulations to analyze the dynamics and differentiate between these modes of subduction margin.” should go higher up in the asbtract, right after the description of the models (before “Our results...”).

- line 18: the authors should specify at which depth is the “radii of curvature” considered, and they should expand on which “observations of present day SZs” are the model predictions confronted.

- lines 25-26: the point “ii) whether large volumes of existing continental crust are ever recycled back into the mantle over long periods of geologic time” is less obvious than the others two (“petrogenesis of continental crust” and “cycling of volatiles”). I suggest that the authors eithr drop it or clarify it (I personnaly think of continental subduction – generating HP eclogitic igneous rocks - when I read this sentence, not of sediment subduction).

- line 34: “The lubricating effect of sediments […] is critical for the mechanism of plate tectonics” is a bit an overstatement > maybe mention only the importance a decoupling interface.

- line 36: what is the “morphology” of a subduction interface = dip, thickness, lateral extent, maximal extent? Same for the accretionary wedges and forearc basins, please specify which properties you are referring to.

- line 43-44: it is a bit unfair to state that “Most of the effort in subduction dynamics studies focused on quantifying dissipation due to slab bending” given the diversity of focus of subduction models, including the ones co-authored by the authors themselves. You should replace by “a lot of effort”...

- line 47-48: the papers cited do not refer to “analytical models” but to “laboratory and numerical models”.

- the introduction is quite lengthy and somewhat repetitive and could be more concise, especially for lines 63-97.

- line 63: it is an overstatement to write “in reality, subduction interfaces are […]” since there is a poor direct access to these deep features. The paper by Agard et al. refers to either indirect geophysical monitoring or to fossil SZ exhumed rocks, which are proxies/models but not “real” samples.

- l. 65-68: this statement is oversimplified and the assertion should be less bimodal (it is not one type of margin or the other). The authors cite for example earlier the paper of Simpson (2010) who states that “rather than viewing compressive plate margins as accretionary versus erosive, the dominant mode may repeatedly switch back and forth through time”.

- the definition of accretion and tectonic erosion should occur earlier in the introduction (presently on line 69 and 80), i.e. before discussing the controlling factor

- when mentioning the influence of convergence rate, it would be relevant to mention that trench retreat (a diagnostic of subduction models), part of convergence rate, is likely to play a role in the margin accretion or erosion.

- line 103: the question “Why some margins accrete sediments while others do not?” is also related to the issue of sediment flux from continents to the margin -depending on topography, climate, river system, oceanic currents – hence the question could be rephrased as “for the same sediment flux at the margin, why do some margin accrete sediments and other not?”

- line 104: the question “How do sediment fluxes influence subduction dynamics and back?”, especially the “and back?” should be better phrased.

- line 104-105: the last question “How should the subduction interface be treated in numerical models, while relaxing the assumption of an interface with constant thickness?” is not really discussed later on in the paper. Please remove.

- line 114: for clarity, the authors should mention some of the “dependent and independent variables” compared between present-day SZs and model outputs.

- line 120: it is confusing that the authors specify that thermal diffusion is neglected when solving the equations… whereas there is no equation of energy conservation (the models are isothermal).

- section 2:specify that the subduction models are mechanical (or compositional) with not temperature variations.

- section 2.1: what are the properties of the upper plate? The upper plate models an oceanic plate, so why is its layering and rheological structure so different from the one of the subducting plate?

Fig.1: it makes no sense to refer to the upper plate mantle material as “slab”, please use a different wording.

- for clarity and conciseness, the introduction of section 2.2 should be directly placed into the subsections “input parameters” and “diagnostic parameters” except lines 180-183.

- line 152: please rephrase “mimicking mid-ocean ridge margins at the tails of the slab” since formation of an oceanic plate take plate at the ridge (whereas “tail” suggest the end), since “slab” only refers to the sunken plate portion and since “margin” often refers to an ocean-continent transition rather than to an oceanic ridge. Suggestion: “mimicking mid-ocean ridges at the trailing edge of the plate away from the trench”

- line 189: specify that the variable thickness mimicks variable ages/thermal structure of the upper plate.

- line 189-191: the sentence “Since the sediment construction is a parameterization of the strength of the oceanic lithosphere, a no-sediment case is still considered when the sediment viscosity is high (higher proportion of crust at the interface).” is incomprehensible, please reformulate.

- line 191-193: the two sentence describing the effect of sediments on mechanical coupling do not belong to the method section but rather to the results or discussion section.

- Line 200-203: the sentence “These variables are compatible with parameters derived from statistical analyses of present-day subduction zones (i.e., Clift and Vannucchi (2004); Lallemand et al. (2005);Wu et al. (2008); De Franco et al. (2008); Heuret et al. (2012)).” is not relevant in the methods section but belongs to the discussion.

-line 207: you mean “Marker 2” instead of “Marker 3”

- line 216: I do not understand the assertion “However, considering that slabs bend elastically, radius of curvature is more appropriate to describe slab deformation with depth, while slab dip represents only the tangent to curvature close to the surface.”, please explain in more details the link with elastic deformation. The real issue seems to me the depth interval up to which the radius of curvature is calculated or over which a mean dip is calculated.

- lines 232-235: it is puzzling that the authors chose to characterize similar bent geometries by two different geometrical diagnostics: a radius of curvature in the case of the slab and a dip in the case if the wedge. I am not sure I understand since slab dip seems a diagnostic also available from SZ data.

- line 482 belong to the methods section 2.1 to justify the range of sediment thickness tested.

- Figure 3e: please specify whether a positive trench rate means advance or retreat.

This manuscript describes a series of 2D geodynamic models that focus on the role of sediment thickness and viscosity on subduction dynamics and accretionary/erosive style of the margin. The study is well-organized and a systematic modeling approach is adopted. The manuscript is well-written and outcomes are relevant to better understand how interface properties affect subduction dynamics. However, I think a more detailed explanation on the physical processes controlling the models behavior is missing. 

Major comments: 

1. Abstract should be more detailed on the results of the study and include an explanation on the processes
2. Since viscosity of sediments is critical, the method section should include more details on the choice of the orders of magnitude of sediment viscosity. Where does the choice of using 15 km-thick weak crust come from? More importantly, since the thickness of the weak crust remains the same in all the simulations and only sediment thickness is varied, have you checked that the thickness of the magmatic weak crust (dark green layer) has no influence on the results? 
3. Lines 279-280: I am not sure that this is a good way to classify a margin as erosive. A probably better way would be to measure the area of triangle ABC and check if this increase or decrease through time. This could be also useful to understand the behavior of unstable accretionary wedge (see following point)
4. The behavior of the model with highly unstable accretionary wedge (SubSed04_100) is really interesting and probably deserves more attention. For example, figure S7 shows the temporal evolution of the convergence rate. It would be interesting to understand if the decrease in convergence rate corresponds to an increase of the sediment subducted and/or variations of subduction geometry (cf. Brizzi et al. 2021)
5. I think the manuscript should include a figure more informative on the role of sediment thickness and viscosity (and upper plate thickness) on convergence velocity, radius of curvature, etc. At present, figure 4 shows correlations among diagnostic parameters  colorcoded by margin style (erosive or accretionary), but I am not fully convinced of the classification criteria used to distinguish them (see point 3).  As the authors acknowledge, the thickness of the sediments is less relevant compared to their viscosity for the margin style (table 2 shows that an erosive margin can form also when the sediment thickness is 10 km, given that viscosity is high). Similarly, the upper plate thickness, which controls the downdip extent of the viscous coupling, seems to be important for e.g. plate velocity (see also Brizzi et al 2021). For example, the high-angle accretionary wedge case makes me think to the influence of the plate interface length, which influences the integrated shear stress. When the wedge is growing, the plate interface is smaller, resistance to subduction is lower and convergence rate is higher. Perhaps plotting the diagnostic parameters as a function of sediment and upper plate thickness, and sediment viscosity could help defining a hierarchy.

Minor comments: 

• Lines 6-7: sentence is unclear. What is the geometry of the global subduction system?
• Lines 14-15:  This could be more specific. How does the viscous coupling relates to the thickness and viscosity of sediments?
• Line 16: this sentence should be probably moved up, before writing about the results. Also, I would suggest to specify what these diagnostic parameters are.
• Line 36: what is the morphology of the subduction interface? Please clarify 
• Line 47-48: Funciello et al 2008 is based on analog experiments. 
• Line 105: This doesn’t seem to be addressed in the paper. I would remove. 
• Line 151: Why did you choose to simulate ocean-ocean subduction instead of ocean-continent? What are the properties of the upper plate and why are they different from the subducting plate?
• Line 155: This sentence is confusing. Please clarify 
• Line 163-164: I would add the tested sediment thicknesses and viscosities here. Also, what is the density of sediments and crust?
• Line 169: What are the subducting plate parameters? Please specify
• Line 188: In my opinion, a sediment thickness of 5 km as the lowest end member is a bit too high. Natural subduction zones can have a (trench) sediment thickness < 1km. Is this choice related to the model resolution? 
• Line 189-191. The sentence is unclear. Please clarify
• Line 207: Marker 2 instead of marker 3
• line 249: are you sure that slab buoyancy is constant? Brizzi et al. 2021 shows that the amount of sediments subducted influences slab buoyancy. 
• Line 291: lubricates and promotes
• Line 297-300: it is really difficult to see what the authors state about interface dynamics from the velocity field in figure 2. The colorbar for velocity is also quite impossible to read. I suggest to use a different colorbar and include a zoomed version of the interface in figure 2
• Line 321: Parentheses for the reference are missing 
• Lines 394-395: the finding that a bigger accretionary wedge promotes high radius of curvature is consistent with results from Brizzi et al 2021