Reply on RC1

The manuscript aims to determine the orogeny-scale (northern Apennine) erosion pattern derived from multiple thermochronometers, and to reconstruct erosion rate variation with space and time. The authors process a large data set of already published apatite fission track and apatite (U/Th)/He data that are accompanied by new detrital AFT data from 7 catchments (modern river sand). Erosion rates have been calculated for each samples using AGE2EDOT code and by applying different values of geothermal gradients.

The paper is well written and in general clear to read. The new data are of high quality. The obtained erosion rates data set is particularly interesting and they worth alone to be published. The application of a kinematic model and an interesting discussion made this paper perfectly suitable to be published in Solid Earth, with a only minor corrections.
The reviewer provides an accurate summary of our study and we appreciate the positive feedback.
My main criticism is focused on the mechanism invoked to explain decreasing in erosion rate along the Ligurian side (the retrowedge of Apennine orogenic wedge). The change in trajectory in the retrowedge seems the first order raison to explain a decrease in erosion rate with time. I feel that the depth of this variation can have a strong impact in the change in erosion rate with time. This variation in trajectory should occur between AFT closure depth and the AHe closure depth. Therefore the closure depth for AFT and AHe systems should deeply control the erosion rate pattern in the retrowedge. In the text it not very clear how closure depth is calculated line 237). Moreover, I am wondering to see the impact of different closure depths in modeling results. This is an important point, and we agree that the methodology for calculating the closure depths should be included. We address this question in the "Kinematic Model" Section of the reviewer's comments.
Regional pattern of several data set (i.e. Ro, fg. 2, thermochronological ages in inset map of fig. 9) shows a clear variation along strike. In the manuscript this along strike variation is never discussed, although has been interpreted in literature as a first-order tectonic control on erosion and exhumation. I would like to know the raison and conditions to apply the same kinematic model of the entire of Apennine wedge.
We briefly explain the pattern of vitrinite reflectance across the orogen and have added additional text to explain the pattern along strike of the orogen.
"Ro values also decrease along strike of the orogen from NW to SE (Fig. 2), illustrating that maximum burial depths also decrease towards the SE. This pattern was in turn interpreted to reflect the shape of the Ligurian Unit as a wedge that thinned towards the east (Zattin et al., 2002), and thus resulted in shallower burial depths for the underlying Cenozoic Foredeep deposits. " We agree that the pattern of vitrinite and cooling ages reflects a first-order tectonic control related to rollback of the Adriatic slab and retreat of the hinge (Thomson et al., 2010) The timing and rate of rollback and hinge rate vary across the orogen, but the fundamental mechanism is the same. In our model, it is not possible to input a spatially variable slab rollback rate, although we do allow for a range of slab rollback rates that is consistent with estimates in our study area. Providing a single rate of rollback may be a model limitation, so we provide some additional text in the discussion to clarify these points.
"The acceleration of exhumation may be related to a change in the timing or rate of slab rollback, which has varied along strike and across the orogen (Faccenna et al., 2014;Rosenbaum and Piana Agostinetti, 2015) and is a first-order tectonic control on exhumation and erosion (Thomson et al., 2010). We allow for a range of rollback rates that are consistent with rates for the field area, although the kinematic model is not able to resolve variability in rollback rates in either space or time."

Line 100 to 102. Variation of Ro is clear to follow also a NW-SE gradient.
The reviewer brings up a good point. We do not discuss the along-strike variability in Ro from the Gottero to the Val D'Arno swaths, although it is clear that 1) Ro values decrease from NW to SE, but that there is also less variability in Ro values on the Adriatic side from NW to SE. We have added extra text to bring up these points.
"Ro values also decrease along strike of the orogen from NW to SE (Fig. 2), illustrating that maximum burial depths also decrease towards the SE. This pattern was in turn interpreted to reflect the shape of the Ligurian Unit as a wedge that thinned towards the east (Zattin et al., 2002), and thus resulted in shallower burial depths for the underlying Cenozoic Foredeep deposits." In the erosion rate result section, I found some difficulties to read the text following figures 7 and 8. Figure 8 is described before figure 7. To be fair, I do not understand the meaning of figure 7 and what information the authors want to explain. It could be useful to add the geographic orientation, i.e NW to SE or NE to SW The purpose of Figure 7 is to illustrate the along-strike differences in ages (related to the reviewer's comment above) and erosion rates for the Adriatic (Figure 7b,d) and Ligurian (Figure 7a,c) sides. This figure is the only example that illustrates the patterns of cooling ages and erosion rates from this orientation, whereas all other figures illustrate data along transects perpendicular to the strike of the orogen. However, we think that this figure is perhaps difficult to read because we have combined all thermochronometers in each panel. To make the figure easier to read, we have created two panels per row, one with the AFT data and one with the AHe data.
We are also more explicit on the difference between Figures 7 and 8 in the text when we introduce the erosion rate results: "Here, we present the erosion rate results for the Adriatic and Ligurian sides, given the two different methods used for constraining the final geothermal gradient, and by illustrating the data with two perspectives: (1) along a profile oriented parallel to the orogen strike (Fig. 7, location shown in Fig. 3) and (2) along swath profiles oriented perpendicular to the orogen strike (Fig. 8)."

Kinematic model.
In this section I suggest to add some lines to describe the code and the environment of modeling. Line 237: It is not clear how the closure depth are chosen.
To calculate the closure depths shown on line 237, we used closure temperatures from the literature: AHe = 70°C (Farley, 2000), AFT = 110°C , (Wagner and Van den Haute, 1992), and ZHe = 180°C (Farley, 2000). These temperatures were converted to a closure depth by assuming a geothermal gradient of 25°C/km. We think that this approach is likely too simplistic, given the constraints we have on final geothermal gradients derived from heat flow maps. Instead, we now use an average final geothermal gradient of 36.4 °C /km for all sample locations in our field area, calculated from the G F_heatflow estimates. Using a higher geothermal gradient will produce shallower isotherms. Given the temperatures listed above for each thermochronometer, this produces closure depths of 1.9 km (AHe), 3.0 km (AFT), and 4.9 km (ZHe). We added the following text to explain our calculations and procedure: "Closure depths were calculated using the closure temperature for each thermochronometer, divided by a spatially and temporally constant geothermal gradient. Closure temperatures are given as: AHe = 70°C (Farley, 2000), AFT = 110°C , (Wagner and Van den Haute, 1992), and ZHe = 180°C (Farley, 2000). Excluding the Alpi Apuane samples, we used the full set of unique sample locations in our field area (Tables 1-4) to calculate an average G f_heatflow = 36.4 °C/km and closure depths for the ZHe (4.9 km), AFT (3.0 km), and AHe (1.9 km) thermochronometers."

At line 372: the best -fit between what data? For large audience could be useful a short description how this model works. Moreover it is not very clear why the authors show this run.
We agree that the term "best-fit" could be clearer for the reader. We add text to both the Methods section and Results section to better describe the objectives of the model and how we found the "best-fit" model results.
Lines 238-242 "We used a range of kinematic and thermal parameters applicable to the Northern Apennines to characterize a kinematic model that aims to: (1) model the path of rock particles from accretion into the wedge to their erosion at the surface, (2) calculate uplift and horizontal rock velocities across the wedge, (3) predict reset cooling ages for AHe, AFT, and ZHe thermochronometers, and (4) calculate maximum burial depths across the model. Here, we describe the model geometry and the kinematic and thermal parameters used to constrain the model."

Lines 383-391
"To construct the best-fit model, our goal is to reproduce a realistic pattern of reset and non-reset thermochronometer ages across the orogen and uplift rates consistent with modern uplift estimates from geodetic releveling (D'Anastasio et al., 2006) for the prowedge (0.5 ̶ 1 km/My) and retrowedge . To this end, we adjust the slab rollback rate within the acceptable range for our field area (6 ̶ 11 km/Ma), and the AHe erosion rates within the range of values calculated from G F_heat flow (0.17 ̶ 1.9 km/My) ( We are not entirely clear about which run the reviewer is referring to. We include the full outputs for both the SCR and VER scenarios. We hope that the above explanation has clarified the reasons for which we included each scenario in our results.
Erosion rate pattern: 459, it could be interesting to specify what kind of tectonic control could be responsible for local high exhumation rate for the Apuane Alps, and to add a reference.
We have added the following text to specify in more detail the tectonic control on the Alpi Apuane, and add the reference from Molli et al. (2018): "An exception to this general pattern may be the Alpi Apuane massif, which represents a structural culmination exposing a deep section and where high exhumation rates from the latest Miocene to the Present likely reflect post-orogenic processes of crustal thinning (Fellin et al., 2007;Molli et al., 2018)." Thank you for pointing this out. This reference is an error and should in fact before Figure  10. We have fixed this mistake. In reference to the scale, we have kept the horizontal scale the same, but have enlarged the vertical scale only so that we can more clearly see the pattern of material motion in the wedge for the depths relevant to the lowtemperature thermochronometers (ZHe, AFT, AHe) that we included in the study .   Fig 9. To make the reading easier, it could be better to move the inset map within the panel 9b.
We agree that the position of the inset makes the figure more difficult to read, so we have moved the inset out of the figure panels and place it above as panel (a). The other panels have been relabeled accordingly and adjusted in the text.