Patsa & Mandal present a paper to study return flow in accretionary wedges, a topic that has been widely studied in the recent decades and is still relevant today. This process has been modelled with theoretical solutions, numerical and analog models to explain field-based observations on pressure temperature conditions on metamorphic rocks along accretionary prisms. The authors provide a more generalized theoretical solution to include a non-parallel component on the slab to reproduce slab advance/rollback. In addition, they conducted analog models to enhance these results. Finally, an overview with natural observations is made.

The integration of analytical solutions with analog modelling represents a novel and valuable approach, with the potential to be further strengthened by incorporating existing numerical results. However, the current manuscript does not fully emphasize this novelty: the introductory section underplays the contribution, and the analytical and analog results are presented somewhat independently, without sufficient cross-comparison.

To improve the manuscript, I recommend (i) revising the introduction to more clearly articulate the novelty and significance of the combined approach, (ii) providing greater coverage and integration of the analog modelling, and (iii) addressing several technical issues within the analytical solution. If these major concerns are resolved, the manuscript would meet the standards for publication in Solid Earth and I look forward to see the revised version.

Major Comment 1

The use of a non-dimensional parameter to quantify the strength or weakness of return flow is both effective and intuitive. This approach allows the balance of incoming and outgoing material along the top boundary to be measured, corresponding to sediment influx and exhumed units in an accretionary wedge. In the simplest case—rigid walls with parallel subduction—the value must equal 1, as dictated by wedge geometry and mass conservation.

However, when a non-parallel component is introduced to the slab, the theoretical models also incorporate material flux along the slab boundary. This introduces bias in the calculated rates of burial and exhumation, since part of the return flow ratio (FR) is influenced by this artificial slab-sourced material. In the case of slab advance, FR appears anomalously high due to two factors: (1) the top-boundary influx of sediments, which the authors correctly identify, and (2) the additional, unaddressed influx of material along the slab boundary.

For the rollback scenario, the formulation permits material to exit through the slab boundary, which is physically unrealistic. With small rollback velocities, wedge geometry would still force material to return to the surface; however, in the present formulation, this is instead channelled out through the slab, leading to FR = 0. Furthermore, at high rollback velocities, the assumption of a downward-tapered wedge may break down, opening the system to the upper mantle. To remain consistent with the wedge geometry assumptions, I suggest restricting the analysis to small perturbations of the non-parallel component.

Finally, the comparison of slab advance, rollback, and normal subduction must be conducted under consistent assumptions regarding material influx. One possible way forward would be to fix the total incoming volume and instead vary return velocities, though I acknowledge that implementing this within the analytical framework may be non-trivial.

Major Comment 2:

I strongly recommend undertaking a systematic comparison between the analytical solution and the analog modelling, as this would greatly strengthen the validation of the theoretical framework. At present, the two sets of results are described independently, which makes the paper look unbalanced. Establishing a one-to-one correspondence between the analytical predictions and the analog experiments for identical setups would provide a more rigorous test of the model and highlight the novelty of the combined approach.

That said, some modifications and additional discussion will be required to enable such a comparison. As noted in lines 415–416, the analog experiments allow material to exit through the bottom boundary, a feature not incorporated into the analytical solution. This discrepancy must be explicitly acknowledged and its implications discussed, since it directly affects the comparability of the results. One option is to consider modified boundary conditions in the analytical framework, or alternatively to constrain the analog results so that they are evaluated under conditions more consistent with the theoretical assumptions.

In addition, I recommend expanding the analog modelling section. At present, the manuscript dedicates considerably more space to the analytical solution, leaving the analog results underdeveloped. A more balanced treatment would not only give greater weight to the experiments but also allow for meaningful side-by-side comparisons. Importantly, you could draw inspiration from the approach of Moulas et al. (2021), who validated their analytical solution against numerical models with a similar setup. Extending your study to include a three-way comparison of analytical, analog, and numerical would significantly increase the robustness and originality of the manuscript.

To facilitate these improvements, some reorganization of the manuscript structure is advised. For example, moving Section 4 earlier in the text, immediately after the presentation of the analytical results would allow for more direct comparisons between the different methods. This restructuring would make the narrative more cohesive and highlight the integrative character of the study, which is currently one of its main strengths but not fully emphasized.

Line to Line comments:

Line 20: Delete “the” or simply state “facilitate subduction.”

Line 21: References are missing for geophysical observations; see Abers (2005).

Lines 25–28: I recommend mentioning the P–T–t path for consistency with the rest of the sentence, and introducing the concept of recycling here.

Line 33: Retain only geochronological, since this sentence refers exclusively to exhumation rates and not geochemical constraints.

Line 39: Verify the reference “?, for review.”

Line 45: A reference is required. The corner flow model also accounts for both prograde and retrograde metamorphism; please mention this. Additionally, note that the model supports the possibility of sediments reaching mantle depths where partial melting may occur.

Lines 48–49: Add one or two sentences on the thermal regime of subduction zones, as this strongly influences eclogite formation.

Lines 56–57: Clarify the rationale of this sentence, or consider removing it.

Line 59: Quantify exhumation rates, providing values from numerical models and natural estimates. Restrict the discussion to the specific tectonic setting under study (wedge geometry)

Line 61: At present, the introduction does not clearly define the scientific gap. While the questions posed are valid, they appear abruptly. The gap would be clearer if you outlined: (i) the discrepancies between modelled and observed exhumation rates, (ii) the wide variability in return-flow models, and (iii) the influence of parameters such as channel width, rheology, and boundary conditions. Emphasize the mismatch between observations and models.

Lines 62–63: Clarify whether questions (1) and (2) are essentially identical.

Line 64: Revise to “theoretical and analog study.”

Line 65: Specify the depth range of both the theoretical and analog models. Although this is mentioned later, it should also appear here. State the main assumptions explicitly—for example, that the accretionary wedge is closed and material cannot enter the mantle—since this represents a special-case scenario.

Line 70: Integrate content from later lines: prior work has already tested different boundary conditions in numerical models (e.g., Gerya et al., 2002). The novelty here lies in the analytical treatment of non-parallel slab boundary conditions, which allows replication of slab rollback and advance. When combined with analog modelling, this provides a unique contribution.

Line 76: Clarify the phrase “oblique to the slab.” If it refers to the trench, rephrase as “a non-parallel component of slab velocity.”

Lines 86–88: This is the first mention of analog experiments, which are central to the manuscript. Introduce them earlier in the introduction and highlight the novelty of combining analytical solutions with analog modelling.

Line 96: Since some material may be dragged down, this is an important limitation—discuss explicitly.

Line 107: Even if non-linearity is not first-order, note that complexity may also arise from contrasting lithologies within the accretionary wedge.

Line 117: Clarify whether this component is oblique to the trench or simply non-parallel.

Lines 120–121: If the trench-oblique component is merely a reduction of subduction velocity, avoid presenting it as a trench-oblique term, since this implies a 3D model.

Lines 127–128: Replace “trench-perpendicular vertical plane” with “non-parallel component of slab velocity.”

Line 130: In Fig. 2 it seems only the oceanic plate’s fixed wall is constrained. Confirm whether the top boundary is also fixed.

Line 181: This derivation follows Moulas et al. (2021). Add a phrase such as “Following the approach of Moulas et al. (2021)...” You may shorten this section and direct readers to that reference until the non-parallel extension is introduced.

Line 214: Same as line 181.

Lines 229–230: Indicate that these models assume either a rigid overriding plate or very strong subduction channels, which generate extremely high overpressures.

Line 245: Clarify whether this is the discretization used to evaluate equations. If so, specify resolution and grid type.

Lines 261–263: Add references or case studies linking models to natural observations.

Line 267: In this simple case, the flow ratio (FR) should equal 1 due to mass conservation and wedge geometry. Clarify whether this depends on discretization.

Line 271: Is this because you are adding material through the slab or because of the "squeezing" of the wedge?

Line 286: Is this difference with the rigid case because of the viscosity ratio only? what if mu_r is even higher (i.e., 10^5 or 10^7), do you reach the rigid wall solution? Also, see later comments to define high/low FR.

Lines 290–291: Figures 4–5 show material still returning to the surface but further from the trench. Confirm whether FR is calculated only for the accretionary wedge (if so, specify at line 245).

Lines 292–294: Consider adding a figure similar to FR vs. obliquity (Fig. 7) to illustrate this result, and extend the same approach to other variables.

Line 323: Since deformation is not described, either remove the vorticity figures or move them to supplementary material with an explanation.

Line 344 and Fig. 7: FR = 0.5 is presented as a threshold between significant and negligible return flow. Explain how this value was determined, or move section 3.2 to the discussion. Kerswell et al. (2023) may provide guidance.

Line 351: Revise “wedge” to “downward-tapered wedge.”

Lines 356–357: Move this sentence to the discussion section.

Lines 392–395: Provide scaling for analog experiments. Do they correspond with analytical models or plate tectonic velocities? Indicate scaling parameters (e.g., Schellart & Strak, 2016) to demonstrate consistency with natural systems and analytical calculations.

Lines 411, 413: Replace “30%” with “0.3U” and “one-sixth” with “U/6” to align with line 415 (“0.6U”).

Line 413: Clarify whether oblique shortening with slab advance and oblique extension correspond to slab rollback. If so, use consistent terminology.

Line 419: Remove “grossly.” Add a one-to-one comparison with the analytical solution.

Lines 422–423: Gravity influences the analog model if the bottom boundary is open. Add a brief discussion of this effect.

Line 425: If slab rollback or oblique extension precludes return flow, clarify how comparisons were made, since the analog and analytical models differ in bottom boundary conditions.

Line 429: Change “crusts” to “crust.”

Lines 449–450: Clarify how comparisons between theoretical models and natural examples are made. Estimate taper angles and subduction dynamics (advance vs. rollback) for each exhumation case, then compare calculated FR values with reviews (e.g., Agard et al., 2009) or case studies (e.g., Franciscan Complex; Ring, 2008).

Line 453: Are there documented cases of absent accretionary wedges in modern rollback settings?

Line 454: Revise: note that HP units are present in Chile (Willner, 2005), with localized pressures of 2–2.5 GPa (González-Jiménez et al., 2017).

Lines 503–504: Mention the thermal regime, as it controls the brittle–ductile transition depth and influences viscosity.

Line 510: Replace “subduction” with “subduction zones.”

Line 519: Expand to “Multiple structural fabrics and fluid-assisted deformation (e.g., Muñoz-Montecinos & Behr, 2023).”

Lines 530–545: This section would benefit greatly from plotting the geological cases in Figure 8, enabling direct comparison between tectonic settings, model outputs, and natural data. Analog model results could also be added for completeness.

Line 590: Include rock strength in the brittle regime. In favorable conditions, rocks can sustain tens to hundreds of MPa before failure (Platt, 2019).

Line 625: Remove mu_r here, as it denotes viscosity ratio.

References in this review:

Abers (2005): https://doi.org/10.1016/j.pepi.2004.10.002 

Agard et al. (2009): https://doi.org/10.1016/j.earscirev.2008.11.002 

Gerya et al. (2002): https://doi.org/10.1029/2002TC001406 

González-Jiménez et al. (2017): https://doi.org/10.1127/ejm/2017/0029-2668 

Kerswell et al. (2023): https://doi.org/10.1029/2022GC010834 

Moulas et al. (2021): https://doi.org/10.1093/gji/ggab246 

Muñoz-Montecinos & Behr (2023): https://doi.org/10.1029/2023GL104244 

Platt (2019): https://doi.org/10.5194/se-10-357-2019 

Ring (2008): https://doi.org/10.1130/2008.2445 

Schellart & Strak (2016): https://doi.org/10.1016/j.jog.2016.03.009 

Willner (2005): https://doi.org/10.1093/petrology/egi035 

In the manuscript “On the criticality of return flows …”, the authors investigate scenarios of rock burial and exhumation during subduction and orogeny. They present analytical solutions for corner flow involving two adjacent deformable wedges: one representing the overriding plate and the other the accretionary wedge (or, in some cases, the subduction channel), where burial and exhumation may occur. The presented analytical model extends the model of Moulas et al. (GJI, 2021). While Moulas et al. considered only subduction velocities parallel to the wedge base, the present study also examines non-parallel velocities. The authors systematically analyze conditions under which significant return flow arises. They further compare the analytical velocity fields with velocity fields from analogue laboratory experiments for similar configurations and show first-order agreement between experimental and analytical velocity fields. Research on burial and exhumation in subduction zones remains highly relevant, as the controlling mechanisms are still debated. Comparing analytical predictions with analogue experiments is also very relevant. However, I have major concerns about the applicability of the presented analytical solution to natural subduction and orogenic burial-exhumation cycles.

Major comments

1) A key assumption of the analytical solution is that wedge geometry remains constant throughout the burial-exhumation cycle. Both the analytical model and analogue experiments yield essentially instantaneous velocity fields. However, the boundary conditions in the analytical model imply that wedge geometry must change over time, as observed in the experiments. If the velocity at the wedge base is not parallel to its base, the wedge will be squeezed or extended. Assuming average burial and exhumation rates of 10 mm/yr, a 50 km burial followed by 50 km exhumation would last ~10 Myr. In the scenarios presented (slab advance and rollback), non-parallel velocity components of just a few mm/yr would displace the wedge base by several tens of kilometers over this timescale. Such large geometric changes would strongly alter the internal corner flow field. Moreover, the issue of corner “stability” was examined by Moulas et al. (2021), but their findings are not fully considered here. For example, Fig. 4 includes results for viscosity ratios ≤1000, yet at such values the internal wedge boundary would deform significantly during the burial-exhumation cycle, effectively destroying the corner geometry. For these reasons, I am not convinced that the presented velocity fields for slab advance and rollback can reliably be integrated to predict a complete burial-exhumation path.

2) Analytical corner flow models can be useful for certain geodynamic scenarios. However, regarding (U)HP rock exhumation in subduction zones, most 2D thermo-mechanical simulations do not generate wedges with forced corner flow when the subduction zone develops self-consistently (i.e., without a pre-imposed weak zone or wedge). In such models, (U)HP exhumation is typically driven by buoyancy or plate divergence/extension, rather than by forced return flow. Thus, the presented corner flow model may be applicable to burial and exhumation at crustal depths, but for rocks buried deeper than ~35 km, exhumation is more likely controlled by buoyancy. Another limitation is the assumption of constant linear viscosity across a wedge spanning the entire crust or even deeper. In reality, significant temperature variations produce large variations in effective rock viscosity. Also, deeper ductile regions may localize strain into shear zones, allowing the subducting plate to slide beneath the wedge without initiating a distributed corner flow. Hence, the authors should clearer discuss the range of applicability of their model. In particular, I find the reference to both an accretionary wedge and a subduction channel in Fig. 1a problematic. An application to shallow accretionary wedges seems far more realistic than to subduction channels in sub-crustal depths.

Minor comments

The introduction is a bit confusing because the authors do not use a consistent terminology. They mention accretionary wedges in places but seem to refer to subduction channels. In Fig. 1a the corner flow seems to be representative for both accretionary wedges and subduction channels. Buoyancy may be ignored in a shallow accretionary wedge but is likely important for subduction channels at sub-crustal depths.

Line 79-90: In section 2.1. the authors state that they generalize the corner-flow model of Moulas et al. (2021), which is based on coupling two corner flow solutions to consider the deformation of the overriding plate. Also, dynamic pressure fields, velocity fields and a velocity profile showing burial and return velocities have all been shown and discussed in Moulas et al. (2021). Also, Fig. 2 is very similar to figures 2 and 3 in Moulas et al. (2021). Maybe the authors could clearer state in the Introduction that their model is a modification of the model of Moulas et al. (2021). In lines 79-90 the study of Moulas et al. (2021) is not mentioned.

Line 103: Buoyancy is excluded but may play the dominant role for rock exhumation from sub-crustal depths in many orogenies.

Line 236-237: I would argue that even more studies consider the positive buoyancy of subducted rocks as a more important mechanisms for the exhumation of (U)HP rocks that have been subducted to depth larger than the average crustal thickness.

Line 270: If I understood correctly, a value of F>1 implies that more material is exhumed by return flow than is subducted. How can this be applied to natural scenarios? It means that more material is exhumed than buried. How is this in agreement with mass balance and where is the additional material that is exhumed generated? F>1 also implies that the wedge disappears after some time. Likely I missed something, but the authors should better explain the meaning, implications and applications for models with F>1.

Line 274: The lower boundary of the wedge in Fig. 3 indicates the top of the subducting plate (see also Fig. 2a). When the velocity arrows are not parallel to this wedge boundary and point towards the subducting plate, then this implies that material from the wedge flows into the subducting plate, because the lower wedge boundary does not move with the velocity. This makes no sense to me. Keeping the wedge geometry constant but applying boundary velocities that would change the wedge geometry is not consistent.

Line 315-316: The mentioned interface migration has been studied by Moulas et al. (2021) with a 2D dimensionless “regime diagram“ distinguishing a “stable corner“ (with negligible interface migration) and an “unstable corner“ (with significant interface migration).

Section 4.1: Please provide a scaling analysis for the experiment so that it is clear how the laboratory experiments can be scaled to the natural situation.

Best regards,

Stefan Schmalholz