Exhumation and erosion of the Northern Apennines, Italy: new insights from low-temperature thermochronometers

Analysis of new detrital apatite fission-track (AFT) ages from modern river sands, published bedrock and detrital AFT ages, and bedrock apatite (U-Th)/He (AHe) ages from the Northern Apennines provide new insights into the spatial and temporal pattern of erosion rates through time across the orogen. The pattern of time-averaged erosion rates derived from AHe ages from the Ligurian side of the orogen illustrates slower erosion rates relative to AFT rates from the Ligurian side and 10 relative to AHe rates from the Adriatic side. These results are corroborated by an analysis of paired AFT and AHe thermochronometer samples, which illustrate that erosion rates have generally increased through time on the Adriatic side, but have decreased through time on the Ligurian side. Using an updated kinematic model of an asymmetric orogenic wedge, with imposed erosion rates on the Ligurian side that are a factor of two slower relative to the Adriatic side, we demonstrate that cooling ages and maximum burial depths are able to replicate the pattern of measured cooling ages across the orogen and 15 estimates of burial depth from vitrinite reflectance data. These results suggest that horizontal motion is an important component of the overall rock motion in the wedge, and that the asymmetry of the orogen has existed for at least several million years.


Introduction
The Apennine mountains of Italy are an active orogen characterized by contemporaneous extensional and compressional tectonics. In the Northern Apennines, these features are linked to rollback of the Adriatic slab beneath Eurasia, suggested to 20 be active since the Oligocene (Malinverno and Ryan, 1986). The interplay between extension and compression has affected the overall tectonic evolution of the Northern Apennines and, in particular, its exhumational and topographic evolution. Lowtemperature bedrock and detrital thermochronology studies have constrained the timing and rates of exhumation at the orogenscale (e.g. Thomson et al., 2010;Malusà and Balestrieri, 2012), and at the regional scale along the extensional retrowedge (Ligurian side) of the orogen (e.g. Fellin et al., 2007) and in the frontal fold-and-thrust belt (Adriatic side) (Balestrieri et al., 25 1996;Carlini et al., 2013;Zattin et al., 2002). Age-elevation profiles and multiple thermochronometers have revealed spatially variable exhumation across and along strike of the orogen, and temporal variability in exhumation rates. While these findings have improved our understanding of the evolution of the Northern Apennines, a comparison between basement and detrital https://doi.org/10.5194/se-2021-96 Preprint. Discussion started: 10 August 2021 c Author(s) 2021. CC BY 4.0 License. data and between exhumation patterns across the primary drainage divide within the frame of an orogenic-scale analysis is lacking. 30 In this paper, we use an updated analysis method that derives long-term erosion rates from cooling ages. We derive timeaveraged erosion rates for individual samples, using two different methods for constraining the initial and final geothermal gradients. We additionally calculate erosion rates through time for existing paired AFT and AHe samples, to compare with results from age-elevation transects (Thomson et al., 2010) that illustrate a change in erosion rates at 4 Ma and were interpreted 35 to reflect an orogen-wide increase in erosion at this time. Our results suggest that the increase in exhumation is restricted to the Adriatic side of the orogen, and may have occurred later (~1-3 Ma), whereas exhumation rates decreased on the Ligurian side at ~1-5 Ma. To understand how this pattern of regional erosion rates relates to orogen-scale kinematics of the Northern Apennines, we propose an updated kinematic model that allows for crustal accretion from both frontal accretion and underplating, and variable temperature at the base of the crust. We find that the pattern of AFT, AHe, and ZHe cooling ages, 40 and the pattern of vitrinite reflectance across the orogen are broadly consistent with the wedge kinematics for an asymmetric orogen that is dominated by frontal accretion and has slower erosion rates on the Ligurian side by a factor of two.

Structural evolution
Development of the Apenninic wedge began at ~30 Ma, due to convergence and southwest-directed subduction of the Adriatic microplate beneath Eurasia. From the late Oligocene, sediments supplied largely by the Alps were deposited as turbidite 45 sequences into a series of northward-migrating foredeep basins (Macigno, Cervarola, and Marnoso-Arenacea Basins) (Fig. 1), which were eventually deformed and thrust during the Neogene (Ricci Lucchi, 1986). Until the Pliocene, these Tertiary foredeep basins were overridden by the Ligurian Unit (Fig. 1), a non-metamorphosed, allochthonous accretionary complex that was thrust upon the Tertiary foredeep deposits as a surficial nappe (Merla, 1952;Pini, 1999). Eocene-to-Pliocene basins formed on top of the Ligurian Unit (Epi-Ligurian Unit) (Ori and Friend, 1984;Cibin et al., 2001), which record discontinuous 50 deposition of shallow-marine and continental sediments (Ricci Lucchi, 1986), and presently exist as denudational remnants above the Ligurian Unit. Today, the Ligurian and Epi-Ligurian Units are the highest structural units exposed in the Northern Apennines.

60
The onset of near-surface exhumation is constrained by the present extent and depositional ages of the Epi-Ligurian Units on the Adriatic side ( Fig. 1), which are commonly not younger than the Tortonian in the NW of the study area; however, to the ESE, near Bologna, they can be as young as the Pliocene (Cibin et al., 2001). The onset of near-surface exhumation in the Northern Apennines is suggested to have begun earlier than 14 Ma, during the Tortonian (Ventura et al., 2001), although the timing of the onset is debated. However, it is clear that rapid exhumation began at 8-9 Ma on the Ligurian side (Balestrieri et 65 al., 1996), and at 4-7 Ma near the divide between the Ligurian and Adriatic sides, based on the ages and younging trend in AFT and AHe thermochronometers <10 Ma ( Fig. 2) (Thomson et al., 2010). https://doi.org/10.5194/se-2021-96 Preprint. Discussion started: 10 August 2021 c Author(s) 2021. CC BY 4.0 License.

Figure 3 (a) Location map for detrital samples.
Detrital AFT samples from this study are illustrated as red diamonds, and published detrital AFT samples are shown as yellow diamonds. Numbers in the gray squares correspond to the age population plots shown in (b). (b) Peak distribution curves (black curves), total PDFs (gray curves), and peak ages for all sampled Ligurian catchments.

Thermochronology data compilation
Cooling ages in the Northern Apennines are primarily limited to apatite fission track (AFT) and apatite (U-Th)/He (AHe) 90 methods (Balestrieri et al., 1996;Ventura et al., 2001;Zattin et al., 2002;Balestrieri et al., 2003;Fellin et al., 2007;Thomson et al., 2010;Malusà and Balestrieri, 2012;Carlini et al., 2013;Balestrieri et al., 2018), as the region is dominated by sedimentary rocks (Fig. 1) that have experienced relatively low burial temperatures of less than 200-250C (Reutter et al., 1983). Non-reset clastic rocks in the Northern Apennines are commonly exposed close to the frontal thrust zone, near the mountain front, and at high elevations (e.g. Zattin et al., 2002;Thomson et al., 2010;Carlini et al., 2013). 95 Maximum burial depths of rock across the Northern Apennines are constrained most commonly from vitrinite reflectance data (Fig. 2;Reutter et al., 1983;Ventura et al., 2001;Botti et al., 2004;Carlini et al., 2013), which is a proxy for burial temperatures recorded by organic particles, and is expressed as Ro (%). Higher Ro values generally reflect higher burial temperatures; Ro increases steadily from NE to SW in the Northern Apennines, with maximum Ro values near the Ligurian coastline, as shown 100 along swath profiles in Fig. 2. This pattern of Ro values was interpreted to reflect NE-directed Miocene thrusting of the Ligurian Unit, which buried the underlying Tertiary foredeep deposits (Reutter et al., 1983). https://doi.org/10.5194/se-2021-96 Preprint. Discussion started: 10 August 2021 c Author(s) 2021. CC BY 4.0 License.

Detrital AFT thermochronology
Bulk samples of modern sand were collected from six rivers on the Ligurian side of the Northern Apennines (Fig. 3a) and are 105 representative of the Macigno, Alpi Apuane, and Ligurian Units (Fig. 1). As some Ligurian catchments (Magra and Serchio Rivers) contain basins with Pliocene sediments, additional samples were collected in tributaries above these basins to avoid sampling the younger, post-orogenic sediments.
Samples were processed according to the external detector method for AFT dating, using standard methods. Bulk samples 110 were sieved and heavy minerals were separated using standard techniques, involving the use of the Wilfley table, heavy liquids, and the Frantz magnetic separator. Apatites were mounted in epoxy and were subsequently polished to expose the internal surfaces of the apatite grains. Multiple mounts per sample were produced to maximize the number of datable grains, and we aimed to date at least 100 grains per sample. However, only 37, 87, and 77 apatites were countable in samples Lima1 (6), Bisenzio (7), and Pescia (8), respectively, whereas the high number of countable apatites in samples Vara (1) and Magra1 (3)  115 allowed us to date 150 grains in each sample.
Apatites were etched in 5.5 N HNO3 for 20 seconds at 21 °C. AFT ages were measured and calculated using the externaldetector and the zeta-calibration methods (Hurford and Green, 1983) with IUGS age standards (Durango and Fish Canyon apatites) (Hurford, 1990). The analyses were subjected to the  2 test (Galbraith, 1981) to assess whether the sample age 120 distributions were over-dispersed; a probability of less than 5% denotes mixed distributions.
We determined age populations for detrital samples based on dominant age peaks identified with the Binomfit program (Brandon, 2002), which is well suited for AFT data with low spontaneous track density. In order to estimate the degree of resetting of the detrital age populations relative to the Apenninic orogenic event, we compared the detrital cooling ages with 125 minimum depositional ages of the Tertiary foredeep units exposed in the drainage areas ( Fig. 1).

Erosion rate analysis
We compiled ages from new and existing detrital AFT samples (23), bedrock AFT samples (139), AHe samples (135), and ZHe samples (26) (Tables 1-4) (Abbate et al., 1994;Balestrieri et al., 1996;Ventura et al., 2001;Zattin et al., 2002;Balestrieri et al., 2003;Fellin et al., 2007;Thomson et al., 2010;Malusà and Balestrieri, 2012;Carlini et al., 2013). As the Alpi Apuane 130 have an erosional history different from the rest of the Northern Apennines (Balestrieri et al., 2003;Fellin et al., 2007), we removed these samples in our compilation of thermochronometric ages.   We converted ages to erosion rates using a half-space cooling model and a closure temperature concept (Willett and Brandon, 2013). This model has the advantage of including an accurate representation of the transience associated with whole lithosphere 150 geotherms. Reset ages were converted to erosion rates using the closure temperature concept (Dodson, 1979), with closure temperatures specific to each thermochronometer, although this is a simplification of diffusional daughter product loss that neglects effects associated with complex cooling histories. For monotonic cooling histories, the measured age of the sample is represented by the time needed for a rock to move from the closure depth to the surface (e.g. Reiners and Brandon, 2006).

155
The conversion to erosion rates was performed using the AGE2EDOT program, (Willett and Brandon, 2013), which estimates an erosion rate through solution of the 1-D thermal advection, diffusion problem for a lithospheric column subjected to a constant rate of erosion. Thermochronometric data required for the calculations include the measured ages and kinetic parameters from which a closure temperature is calculated. In addition, the thermal initial and boundary conditions, as well as thermal parameters, must be specified for each sample site. 160 For the kinetic parameters for AHe, we assumed grain sizes of 45 µm, given that sizes of dated grains are not reported by previous studies, and that a grain size of 60 µm is larger than the mean size of detrital apatites that are typically dated in the Northern Apennines. Thermal parameters (see definitions in Table 5) include an estimate for the age of onset of erosion (t1); https://doi.org/10.5194/se-2021-96 Preprint. Discussion started: 10 August 2021 c Author(s) 2021. CC BY 4.0 License. the sample elevation, given as an elevation above a regional mean (h); surface temperature (Ts); and either an initial geothermal 165 gradient (G0) or the final geothermal gradient (Gf). Only one estimate of the geothermal gradient is needed, but we took several approaches, as will be discussed below. To calculate Ts, we adjusted a base temperature value for the elevation of each sample, given a lapse rate of 5C /km. For the base temperature, we used a modern surface temperature of 13.8C, which represents the calculated yearly average for an elevation of 53 m at Bologna from 1813-2004 (NOAA Global Temperature Summary of the Year dataset). 170 Table 5 Definitions of thermal parameters used in the erosion rate analysis.
The geothermal gradient is the most important parameter incorporated into the erosion rate analysis and is also the largest 175 source of uncertainty. It can be specified either as a final geothermal gradient (Gf), which is the present geothermal gradient at the surface, or as an initial geothermal gradient (G0) that is assumed to be constant with depth at the onset of exhumation (Willett and Brandon, 2013). We calculated and compared erosion rates derived using two approaches. In the first method, we imposed a spatially constant G0 of 25C/km (G0_25) (Balestrieri et al., 2003;Ventura et al., 2001;Zattin et al., 2002). In the second method, we assumed that the present-day geothermal gradient (Gf) matches the geothermal gradient calculated from 180 geothermal heat flow measurements. We converted the heat flow measurements to a Gf (Gf_heatflow) using a spatially constant thermal conductivity value for sandstone (2 W/mK). Heat flow values were extracted from contour maps that interpolate geothermal well data (Pauselli et al., 2019;della Vedova et al., 2001) The modelling procedure described above was applied to all ages, assuming that erosion initiated over the entire region at 10 190 Ma. The resulting erosion rate applies from the onset of exhumation at 10 Ma to the present and reflects the time-averaged https://doi.org/10.5194/se-2021-96 Preprint. Discussion started: 10 August 2021 c Author(s) 2021. CC BY 4.0 License. erosion rate constrained to pass through the closure temperature at the age and with a cooling rate commensurate with the average erosion rate. Thus, this method is limited to a single, average erosion rate. However, changes in exhumation rates through time in the Northern Apennines are supported by several lines of evidence, particularly by age-elevation transects (AETs). In fact, AETs from the existing literature illustrate differences along the age-elevation slope for a single 195 thermochronometer (as in Balestrieri et al., 1999) or among age-elevation slopes for multiple thermochronometers (as in Thomson et al., 2010).
It is possible to use AGE2EDOT in an incremental manner, allowing us to use paired thermochronometers analyzed from a single sample. In this case, the temporal range of exhumation is bracketed by the AFT and AHe ages, with independent erosion 200 rates determined from each age, thus resolving two time intervals (Willett et al., 2020). In principle, this violates the assumption of a constant rate of cooling implicit to use of the closure age concept, but provided that the transition between erosion rate intervals is not close to either age, the error will be small. We analyzed 30 available paired ages to detect temporal changes in erosion rate. For the paired ages analysis, the exhumation path is divided into two segments: the first segment extends from the onset of exhumation to a specified transition time (ttrans) after cooling through the AFT system, and the second segment 205 extends from this transition time to the present, thus passing through the AHe age in this second interval (Fig. 4). We derive an erosion rate for each of these time-segments by analyzing each segment with AGE2EDOT, linking the two solutions at ttrans.
The solutions are matched by noting the depth and temperature of the sample at ttrans, based on the erosion rate in the second interval, and using this and the geothermal gradient at ttrans as the boundary conditions for calculations of the first interval (Fig.   4). 210 The difference in age between some of our paired ages is less than 1 Ma, but larger than 0.5 Ma, so we set the transition at 0.5 Ma before the AHe closure for all samples (i.e. the AHe cooling age plus 0.5 Ma), in order to allow the onset of advection to precede the AHe closure. Calculation of the erosion rate over the second interval requires the modern surface temperature (TS); the sample elevation above the regional mean (h); the AHe cooling age (tau); the final geothermal gradient as derived from 215 heat-flow measurements (Gf_heatflow); and the length of the time interval (ttrans), calculated as the AHe age plus 0.5 Ma. The erosion rate is then solved from these data and the kinetic parameters.
To calculate the erosion rate for the first interval, we require in addition, the temperature at ttrans (T0) (Fig. 4); the sample elevation above the regional mean (h); the sample AFT age (tau), and the time of onset of erosion, taken as 10 Ma. To match 220 solutions at ttrans, we simply reduce the age by ttrans, and reduce the elevation by the amount of exhumation that occurred during the second interval. We take the initial geothermal gradient obtained from the model for the second interval as the final condition for the first interval.

Figure 4
Schematic of paired ages erosion rates analysis, illustrating the theoretical temperature path through time for a sample with paired AFT and AHe cooling ages, given the thermal parameters described in the text. ė represents the erosion rate for each interval.

Kinematic model
The kinematic model presented here approximates the Northern Apennines as a doubly tapering, asymmetric wedge, given the 230 geometric parameters illustrated in Fig. 5. The Adriatic and Ligurian sides of the orogen are defined as the accreting prowedge and non-accreting retrowedge of the orogen, respectively (e.g. Willett et al., 2001). The geometry of the wedge is defined by surface and basal angles for the prowedge (P and P) and retrowedge (R and R). The lengths of the prowedge (LP) and retrowedge (LR) are 60 km and 40 km, respectively, based on average widths measured from an SRTM 90 Digital Elevation Model (DEM). The maximum crustal thickness is 56 km (Spada et al., 2013), the maximum elevation is 2 km, and the thickness 235 of the accreted crust is 20 km, partitioned between frontal accretion (h0 = 10 km) and prowedge basal accretion (h1 = 10 km).
Material is accreted to the wedge through thrusts slices in the upper plate (frontal accretion) or is offscraped from the subducting plate at depth (underplating). Material motion is constrained by balancing frontal and rearward fluxes, underplating, 240 and erosion. We prescribe a compressional prowedge and an extensional retrowedge, where horizontal velocities decrease along the prowedge and increase along the retrowedge as a function of distance. The vertical rock velocity is also variable with depth, and is defined as the sum of the erosion rate and a component of crustal thickening driven by accretion.
The velocities in the model are defined as follows: plate subduction velocity (VP), prowedge underplating velocity (UP), 245 prowedge erosional velocity (eP), and retrowedge erosional velocity (eR). The plate subduction velocity, or convergence rate, for the Northern Apennines is suggested to be driven entirely by slab rollback, so we used estimates of slab rollback to parameterize the convergence rate. Slab rollback rates are on the order of 6-10 km/My in this region of the Apennines https://doi.org/10.5194/se-2021-96 Preprint. Discussion started: 10 August 2021 c Author(s) 2021. CC BY 4.0 License. (Faccenna et al., 2014;Rosenbaum and Piana Agostinetti, 2015), so we run the model using these minimum and maximum values as end-member scenarios. We also vary the spatial pattern of erosion rates in the model using two model assumptions: 250 1) a constant erosion rate across the orogenic wedge (SCR), and 2) a higher erosion rate on the prowedge relative to the retrowedge (VER).

Figure 5 Kinematic model of the Northern Apennines as an orogenic wedge with internal deformation driven by frontal and basal accretion and surface erosion. Mass is balanced to maintain a steady size, and internal deformation is calculated to be consistent
255 with boundary conditions.

Detrital AFT cooling ages
New detrital AFT (8) sample ages are given in Fig. 3b and Tables 6-7. Central ages vary from 5.4  0.6 Ma to 10.5  0.7 Ma, 260 and single grain ages show a wide range of values from 5.1 to 145.3 Ma. All samples except Lima1 show at least two distinct age populations, with minimum age peaks between 5.1 and 8 Ma. All minimum age peaks are younger than the stratigraphic ages of units within the catchment (Cita Sironi et al., 2006;Delfrati et al., 2002;Pialli et al., 2000), with the exception of Magra1. This site contains Plio-Pleistocene deposits exposed in its catchment that are younger than the minimum age peak, but these are locally derived sediments, and the sedimentary bedrock has stratigraphic ages older than the young peak. 265 In the five southern samples (Serchio, Lima1, Lima2, Pescia and Bisenzio), the youngest peaks represent the largest age populations and are similar to the sample central ages. The three northern samples (Vara, Magra1, Magra2) show two common age populations at 5-6 Ma and at 12-13 Ma and have central ages older than the minimum peak ages, due to a large proportion of older grains. 270

Geothermal gradients and erosion rates
We report initial geothermal gradients (G0) and final geothermal gradients (Gf) using the two approaches described in the methods. Given a G0 = 25 C/km common to all samples, Gf_25 ranges from 27.4 to 55.2 C/km for AHe samples (Table 8) and ranges from 31.2 to 49.7 C/km for AFT samples (Table 9). Using the second method based on modern heat flow 280 measurements, G0_heatflow for AHe samples ranges from 9.3 to 42.2 C/km (Table 8)  AFT samples (Table 9). Relative to the Gf_25, Gf_heatflow derived from Pauselli et al. (2019) are consistently lower, where all samples lie left of the 1:1 trendline for AHe samples (Fig. 6e) and all but one lie left of the 1:1 trendline for AFT samples (Fig.   6b). In contrast, Gf_heatflow derived from della Vedova et al. (2001) are highly variable, although the majority lie to the right of the 1:1 line for both AHe (Fig. 6e) and AFT (Fig. 6b) samples, indicating that these values are higher relative to Gf_25. 285 Erosion rates calculated using the two methods for estimating geothermal gradients also illustrate different trends for the della rates derived from G0_25, but always by a factor of less than two (Fig. 6c,f).

310
Here, we present the erosion rate results for the Adriatic and Ligurian sides of the orogen. Calculated with G0_25, erosion rates inverted from AFT bedrock ages (Table 9; top row, Fig. 8) vary between 0.41 km/My and 1.19 km/My on the Adriatic side ( Fig. 7d) and between 0.36 km/My and 0.84 km/My on the Ligurian side (Fig. 7c). The highest erosion rates on the Ligurian side are located in the Macigno Unit near the Alpi Apuane, whereas the highest rates on the Adriatic side are located near the drainage divide (Cervarola Unit; top row, Fig. 8). AFT bedrock erosion rates across the divide are similar or slightly higher on 315 the Adriatic side ( Fig. 7c-d, top row, Fig. 8). AFT detrital erosion rates are similar to bedrock AFT rates, which also exhibit erosion rates that are higher on the Adriatic side (top row, Fig. 8). Erosion rates derived from AHe ages (Table 8) range between 0.14 and 1.39 km/My on the Adriatic side and between 0.14 and 0.74 km/My on the Ligurian side. AHe erosion rates on the Adriatic side are more variable relative to the Ligurian side, particularly in the southeast region of the field area ( Fig. 7c-d).
https://doi.org/10.5194/se-2021-96 Preprint. Discussion started: 10 August 2021 c Author(s) 2021. CC BY 4.0 License. Similar to the bedrock AFT erosion rates, the highest AHe erosion rates are found on the Adriatic side near the drainage divide 320 and are lowest near the Ligurian coastline (top row, Fig. 8).
Calculated with Gf_heatflow, the pattern of erosion across the drainage divide is similar to the pattern for erosion rates calculated with G0_25 (middle row, Fig. 8). On the Adriatic side, erosion rates inverted from AFT bedrock and detrital ages vary between 0.34 and 1.63 km/My and between 0.88 and 1.44 km/My, respectively. On the Ligurian side, AFT bedrock erosion rates vary 325 between 0.26 and 1.28 km/My, and detrital AFT erosion rates vary between 0.34 and 0.78 km/My. Erosion rates derived from AHe ages range from 0.17 to 1.92 km/My on the Adriatic side and from 0.10 to 1.02 km/My on the Ligurian side. Detrital AFT erosion rates on the Adriatic side are higher relative to the Ligurian side, regardless of the method used for constraining the geothermal gradient. However, calculated from Gf_heatflow , detrital AFT erosion rates on the Adriatic side are up to a factor of two higher than erosion rates calculated with G0_25 (Fig. 8). 330 https://doi.org/10.5194/se-2021-96 Preprint. Discussion started: 10 August 2021 c Author(s) 2021. CC BY 4.0 License. Fig. 3a. In (b), individual AHe symbol dated to <1 Ma is illustrated with a double gray slash.

Figure 7 Cooling ages for all thermochronometers plotted along the Adriatic side (a) and the Ligurian side (b). Erosion rates for the (c) Adriatic side and (d) Ligurian side. Profile location for all plots is illustrated in
https://doi.org/10.5194/se-2021-96 Preprint. Discussion started: 10 August 2021 c Author(s) 2021. CC BY 4.0 License. Fig. 3a. In the Bologna swath profile (right column), one AHe sample could not be resolved for an erosion rate with Gf_heatflow.

Paired ages 340
Of the 30 paired samples analyzed here, erosion rates for two samples could not be resolved (Table 12), due to the similarity in ages between the AFT and AHe thermochronometers (sample C16) or due to an AHe age that is older than the AFT age (sample C22). Six paired samples are located on the Ligurian side of the orogen, and the remaining 22 samples are located on the Adriatic side of the orogen (Fig. 9).

345
Erosion rates from samples on the Adriatic side vary from ~0.3 to 5.2 km/My (Table 12). Twelve samples illustrate an increase in erosion through time (Fig. 9a), with the shift generally occurring between ~1 and 2 Ma. Ten samples from the Adriatic side illustrate a decrease in erosion rate through time (Fig. 9b), with the shift generally occurring between ~2 and 4 Ma. With the exception of three samples (AP53, AP57, and C29), decelerating sites are located in the headwaters of the Reno River (inset map, Fig. 9c), which drains the Cervarola and Macigno Units east of Mt. Cimone (Fig. 1). For the six paired samples from 350 the Ligurian side, the range of erosion rates is similar and varies from ~0.3 to 5.1 km/My (Table 12). However, only one sample on the Ligurian side illustrates an increase in erosion rate, occurring at ~3 Ma (Fig. 9c). All other samples illustrate a decrease in erosion through time, with the shift occurring between ~6 and 2 Ma (Fig. 9d).

Figure 9 Erosion rates through time from paired AFT and AHe age samples on the Adriatic side for (a) increasing erosion rates and (b) decreasing erosion rates, and on the Ligurian side for (c) increasing erosion rates and for (d) decreasing erosion rates. Inset map in (c) shows the locations of paired thermochronometer age samples (circles). Red circles indicate locations where the erosion rate through time is increasing, and blue circles indicates locations where the erosion rate through time is decreasing.
360

Kinematic model
The orogenic wedge model shows how the spatial pattern of exhumation rates relates to the polarity of accretion and the pattern of horizontal and vertical motion. Figure 10 illustrates the predicted horizontal velocities, uplift rates, and material paths 365 through the wedge for the spatially constant erosion rate setup (SCR) and the spatially variable erosion rate setup (VER).
Horizontal velocities at the toes of the wedge are equal to the rate of slab rollback and decrease to a minimum at the drainage divide between the prowedge and retrowedge (Fig. 10a). Extension in the retrowedge results in higher horizontal erosion rates towards the Ligurian coast.

370
For the SCR setup, we ran the model with a single, orogen-wide erosion rate varying between 0.4 and 1.0 km/My, which are values consistent with average AHe-derived erosion rates for the Adriatic side of the orogen. Here, we illustrate the best-fit model results for an orogen-wide erosion rate of 0.8 km/My and a slab rollback rate (VP) = 10 km/My. Horizontal velocities decrease to a minimum of 2.0 km/My at the drainage divide (Fig. 10a). Uplift rates on the prowedge vary from 0.89 to 1.1 km/My, and from 0.3 to 0.66 km/My on the retrowedge (Fig. 10b). Predicted reset cooling ages across the orogen for AHe 375 samples (triangles) vary from 2.8 to 4.4 Ma, from 4.2 to 7.6 Ma for AFT (diamonds), and from 6.5 to 11.5 Ma for ZHe samples https://doi.org/10.5194/se-2021-96 Preprint. Discussion started: 10 August 2021 c Author(s) 2021. CC BY 4.0 License.
(circles) (Fig. 11a, c). Maximum burial depth increases almost linearly across the orogenic wedge, ranging from 2.1 to 6.1 km on the prowedge and from 6.9 to 13.1 km on the retrowedge (Fig. 11b).
For the VER setup, we illustrate the best-fit model results for a prowedge erosion rate of 0.8 km/My and a retrowedge erosion 380 rate of 0.4 km/My. Horizontal velocities decrease to a minimum of 2.0 km/My at the drainage divide (Fig. 10D). Uplift rates on the prowedge are the same as in the SCR (0.87-1.1 km/My), but uplift rates on the retrowedge are lower (-0.1-0.26 km/My) (Fig. 10e). Predicted reset cooling ages range from 3.8 to 5.8 Ma for AHe samples, from 6.2 to 8.3 Ma for AFT samples, and from 12.1 to 12.8 Ma for ZHe samples. For samples that reach the surface on the prowedge side of the range, predicted reset cooling ages decrease towards the primary divide. On the retrowedge, ages initially increase away from 385 the divide, but young slightly at the model boundary (Fig. 11d). Maximum burial depth increases almost linearly along the prowedge (2.9-6.1 km) up to a kink at the drainage divide, which reflects the shift towards a shallower slope and smaller range of maximum burial depths along the retrowedge (6.9-9.4 km) (Fig. 11e).

Detrital versus bedrock ages
Previous studies place the onset of exhumation between 8 and 14 Ma (Balestrieri et al., 1996;Ventura et al., 2001), so it is not 400 clear whether the 12-13 Ma old population in the Vara and Magra samples represent partially or completely reset cooling ages.
However, these ages are consistent with high-elevation samples west of the Vara catchment that record slow cooling prior to 8 Ma. The 8.2 Ma age peak is present in the Magra1 sample, which drains the extensional intermontane basin within the catchment, but is absent in the Magra2 sample, which drains only small tributaries upstream of the basin. Thus, the peak at 8.2 Ma in Magra1 likely reflects exhumation ages of the nearby Macigno Unit, which would have been eroded and redeposited 405 https://doi.org/10.5194/se-2021-96 Preprint. Discussion started: 10 August 2021 c Author(s) 2021. CC BY 4.0 License. into the Pliocene basins (Balestrieri, 2000;Fellin et al., 2007). The pulse of exhumation in the Northern Apennines between 6 Ma and 4 Ma, when most of the Northern Apennines became sub-aerially exposed (Zattin et al., 2002;Balestrieri et al., 2003;Fellin et al., 2007), is consistent with the youngest peaks shown in the Vara, Magra, Lima, and Bisenzio Rivers.
The pattern of detrital AFT ages on the Ligurian and Adriatic sides (Malusà and Balestrieri, 2012) is consistent with the pattern 410 of bedrock AFT ages, which illustrate younging exhumation ages towards the northeast, regardless of whether cooling ages are corrected or uncorrected for topography (Fig. 2). Overall, we find consistent results between the detrital AFT and bedrock AFT ages, reinforcing that the reset detrital ages illustrate a true exhumation signal, rather than an artefact of the technique.
Fertility analysis of sediment from sampled Adriatic catchments also confirm that the detrital samples are representative of the eroded bedrock (Malusà et al., 2016), in the absence of hydraulic sorting effects. Finally, we note that using an initial 415 exhumation age of 14 Ma would not change the pattern of exhumation, but would proportionally decrease all erosion rates in the Northern Apennines if the geothermal gradient increases in response to a longer period of erosion. As we incorporate only minimum reset ages from each sample, inclusion of the 13 Ma-year age population would have no effect on the erosion rate calculations.

Inferred erosion rate and relationship to the geothermal gradient 420
The geothermal gradient is the most important external parameter for converting cooling ages into erosion rates (Willett and Brandon, 2013). The two major studies of regional heat flow (della Vedova et al., 2001 andPauselli et al., 2019) show large variations in heat flow across the region, but are also internally inconsistent by up to 50 mW/m 2 in the regions where they overlap. It is thus unclear how much of the spatial variability is real and how much is due to local effects or local errors in heat flow measurements, a large source of uncertainty for the geothermal gradients and ultimately, the erosion rates. To address 425 this uncertainty, we took two approaches. First, we assumed that the initial geothermal gradient in the region was uniform and all variations in the modern geothermal gradient are due to advection in response to erosion. Second, we constrained the thermal model to be consistent with the modern heat flow measurements and inferred an initial geothermal gradient that was spatially variable.

430
The uncertainties in the modern heat flow measurements are evident in the erosion rate analysis, particularly when comparing the range of G0_heatflow inferred from Gf_heatflow inputs that are calculated with heat flow measurements (Fig. 6a, d). However, it is unclear whether the large range of G0_heatflow values represents how the geothermal gradient may have varied in either space or time at the onset of erosion. As the Northern Apennines evolved, sediments were accreted to the accretionary wedge shortly after being deposited in a subsiding foreland basin, whose modern equivalents are the Po Plain and the Adriatic Sea. There, 435 modern heat flow values are generally low (≤ 50 mW m -2 ) (della Vedova et al. 2001), although with significant spatial variations (Pauselli et al., 2019), indicating that the present geothermal gradients in the foreland should be not higher than https://doi.org/10.5194/se-2021-96 Preprint. Discussion started: 10 August 2021 c Author(s) 2021. CC BY 4.0 License. about 30 ºC/km. Given the uncertainties on the modern heat flow measurements, the erosion analysis approach based on a common, G0_25 can be considered a viable alternative to the approach based on an input Gf_heatflow.

440
The erosion rates resulting from these two analyses differ significantly: erosion rates from one analysis are a factor of two different from the alternate analysis (Fig. 6). However, the two sets of results projected along swath profiles from SW to NE show little difference in their spatial patterns across the main Apenninic divide (Fig. 8). The main differences are that the erosion rates derived from Gf_heatflow vary over a larger and higher range than those derived from G0_25, and the maximum rates are higher from Gf_heatflow. In particular, the youngest detrital age populations give much larger rates on the Adriatic side than 445 on the Ligurian side with the analysis based on Gf_heatflow. These observations suggest that the erosion rates derived from G0_25 may be more conservative estimates overall. However, the most important observation for the scope of this contribution remains that the large-scale spatial pattern of erosion rates along the swath profile does not change with the employed analysis method.

Erosion rate patterns 450
Bedrock cooling ages on the Ligurian side of the Northern Apennines generally vary between 4 and 10 Ma, with only a few ages younger than 4 Ma. On the Adriatic side, bedrock cooling ages younger than 4 Ma are a large component of the age distributions, especially among the AHe ages (Fig. 7b). Similarly, the youngest populations for detrital AFT samples on the Ligurian side are nearly all older than the youngest detrital AFT populations on the Adriatic side ( Fig. 7a-b).

455
Erosion rates derived from both bedrock and detrital thermochronometric ages suggest a difference between the Ligurian and Adriatic sides that is valid at the regional scale, regardless of the method used for constraining the geothermal gradients. An exception to this general pattern may be the Alpi Apuane massif, which represents a structural culmination exposing a deep section. It is likely that this exhumation reflects structural control that is unique to the Northern Apennines. On the Ligurian side, erosion rates derived from bedrock AFT ages tend to be higher than erosion rates obtained from bedrock AHe ages (Fig. 460 7C) reflecting a regional decrease in erosion rate. This is particularly evident in the region east of the Alpi Apuane, at the main drainage divide north of Florence and in the Val d'Arno ( Fig. 1 and 2). In contrast, on the Adriatic side, erosion rates derived from AHe ages tend to be higher than erosion rates obtained from AFT ages (Fig. 7D) and suggest an increase in erosion rate over the last 5 Ma.

465
Paired thermochronometers on the same sample (as for instance, AFT and AHe) or age-elevation transects (AETs) also indicate changes in erosion rate. The majority of the paired age samples (12) from the Adriatic side illustrate an increase in erosion rate through time (Fig. 9a), although 10 of these samples illustrate a decrease in erosion through time (Fig. 9b). Of these 10 samples, seven are from one region of the upper Reno River Valley (inset map, Fig. 9c). The headwaters of the Reno River Valley extend farther south than the adjacent basins of the Serchio River and Bisenzio River, which flow to the Ligurian Sea (Fig. 470 https://doi.org/10.5194/se-2021-96 Preprint. Discussion started: 10 August 2021 c Author(s) 2021. CC BY 4.0 License. 3a). Interestingly, the exhumation rates from the upper Reno River are similar to rates from the Serchio River, suggesting that the upper Reno River presents an erosion rate signal akin to Ligurian Rivers, rather than to Adriatic Rivers, and are thus resolving a consistent pattern of erosion rate in space, but not restricted to catchment boundaries. We also note that modern erosion rates from cosmogenic nuclide concentrations in the upper Reno tributaries are at least a factor of three lower than rates for the entire basin (Cyr et al., 2010), suggesting that this trend of lower erosion rates in this area has continued to the 475 present.
Paired ages from the Ligurian side show the opposite trend. With the exception of one sample (C34) (Fig. 9c) located within the Magra2 catchment area (Fig. 3a), all other samples consistently illustrate a decrease in erosion rate through time (Fig. 9d).
Thus, the results from the paired thermochronometer ages on the Adriatic and Ligurian sides of the orogen confirm 480 the regional trends observed from the simple erosion rate analysis method.
The results from our paired ages analysis can also be discussed in the context of the AETs from Mt. Falterona, Mt. Cimone, and Val d'Arno (see Fig. 1 for locations). Between 4 and 2 Ma, these AETs have previously been interpreted to reflect an orogen-wide increase in exhumation and erosion rates, although there are notable differences between the results from the 485 profiles on the Adriatic side (Mt Falterona and Mt. Cimone) and on the Ligurian side (Val d'Arno) (Thomson et al., 2010).
The Mt. Falterona and Cimone AETs illustrate a two-fold increase in erosion rates between 4 and 5 Ma, from 0.29 ± 0.1 km/My to 0.58 ± 0.23 km/My and from 0.22 ± 0.09 km/My to 0.58 ± 0.16 km/My, respectively (Thomson et al., 2010). Excluding the samples from the upper Reno River Valley that illustrate a decrease in erosion rate through time, our paired ages analysis reflects erosion rates that increase from 0.68 ± 0.42 km/My to 1.31 ± 0.70 km/My, given 1s uncertainties. Our average erosion 490 rates are higher than the average erosion rates calculated for the Mt. Falterona and Cimone AETs, although, given the high uncertainties in our values, they are within the range of the AETs.
The results from the Val d'Arno AET are less straightforward, due to the fact that some samples illustrate a decrease in erosion rate through time, while others illustrate an increase in erosion rate through time. When corrected for topographic and advection 495 effects, this AET shows a negative slope that was previously interpreted to reflect post-cooling tilting of the footwall block of an extensional fault (Thomson et al., 2010). On the Ligurian side, cooling ages and erosion rates vary locally as a function of elevation and of fault activity, and extensional faults can control differences in the exhumation pattern. However, in light of our results from the simple analysis of erosion rates and the paired ages, this negative slope could also be interpreted as a decrease in erosion rates, and thus would reflect a regional signal, rather than local tectonics. We infer that such regional scale-500 differences must be controlled by first-order features of the Northern Apennines. In order to address the question of what could control such differences, we compare two different kinematic models for the Northern Apennines orogenic wedge.

Kinematic model
The orogenic wedge kinematic models illustrate differences in cooling ages, maximum burial temperatures, and material paths across the Northern Apennines, assuming simple continuum accretion and mass balance (Fig. 11). Using a spatially 505 constant erosion rate across the orogenic wedge (SCR) predicts that reset ages decrease from northeast to southwest and are youngest on the retrowedge model boundary. In contrast, the variable erosion rate setup (VER) predicts minimum reset ages near the drainage divide and maximum reset ages on the retrowedge, close to the model boundary. The VER model is consistent with observed cooling ages, which are youngest near the drainage divide in the core of the Northern Apennines ( Fig. 2). 510 Vitrinite reflectance (VR) data provide an additional estimate for maximum paleotemperature and burial depth and thus, an additional calibration of the kinematic model. In the Northern Apennines, Ro values reach 5.1% at the Ligurian coastline along the Mt. Gottero swath profile (Fig. 2a). With the exception of this profile, maximum VR values are generally within the range of 1.5-2.5 % for the Mt. Cimone, Bologna, and Val d'Arno profiles ( Fig. 2b-d). Maximum paleotemperatures from 515 VR are estimated at 200-250C in the core of the range and along the Ligurian coastline in the northwest (Fellin et al., 2007), whereas paleotemperatures are 150-190C in the Cervarola Unit (VR = 1.0-1.7%), and are 100-110C in the Ligurian Unit (VR = 0.5-0.6%) (Ventura et al., 2001;Botti et al., 2004). Maximum paleotemperatures should correspond to maximum burial depths; thus, we expect to find the maximum burial depths along the Ligurian coastline and near the drainage divide in the Cervarola Unit. Both the SCR and VER models predict maximum burial depths near the Ligurian coastline, consistent 520 with the trends observed in the Mt. Gottero and Val d'Arno profile. Near the drainage divide on the prowedge, predicted maximum burial depths are 6.1 km for both models.
We can estimate maximum burial depths at the drainage divide, given the generalized relationship between Ro and burial depth (Suggate, 1998), VR values, and a modern geothermal gradient. To estimate the modern geothermal gradient at the 525 drainage divide, we use Gf_heatflow rates from the simple erosion analysis for AHe samples in the Cervarola Unit (Tables 1 and   7). Given these parameters. we estimate maximum burial depths in the range of 4.0-5.5 km near the drainage divide; the upper bound of this range are similar to the predicted maximum burial depths (6.1 km) for both the SCR and VER models.
Collectively, our erosion rate analysis and kinematic model illustrate that the east to west particle trajectories, combined with lower erosion rates on the retrowedge, are consistent with the spatial pattern of cooling ages and maximum 530 paleotemperatures estimated from both vitrinite reflectance and thermochronometric data.
The particle paths in the VER kinematic model, combined with the lower erosion rates in the retrowedge, suggest an explanation for the apparent decrease in erosion rates with time on the Ligurian side of the Apennines. As rocks are advected from prowedge to retrowedge, the vertical component of their motion decreases (Fig. 10f). Particle paths where the AFT age 535 https://doi.org/10.5194/se-2021-96 Preprint. Discussion started: 10 August 2021 c Author(s) 2021. CC BY 4.0 License. is set in the prowedge, but where the AHe is set in the retrowedge, will record this change as a temporal deceleration of cooling rate. However, rather than representing a change in surface erosion rate, the change in cooling rate reflects the motion of the rock from the fast erosion rate prowedge into the low erosion rate retrowedge. The fact that the decelerating sites are all found to the southwest, regardless of drainage basin (inset, Fig. 8c) supports the idea that this is a tectonically controlled spatial pattern. 540 The acceleration of exhumation observed in the prowedge is not explained by the kinematic model. Although the apparent increase in exhumation rate might be explained by spatial changes in tectonic uplift and an associated increase in erosion rate, there is no strong evidence for this in the spatial pattern of ages or in the geomorphology, which should show higher uplift rates in the range interior. In contrast, the highest uplift rates are more often observed at the mountain front (Picotti and 545 Pazzaglia, 2008). It is more likely that this is a true temporal increase in erosion rate, which could be associated with an increase in accretionary flux as the mountain front advanced into the Alpine sediments of the Adriatic foreland. The foreland basin fill thickened as Miocene alpine sediments filled the foredeep and again in the Quaternary as glacial sediments filled the Po plain and parts of the Adriatic Sea. The increase in accretionary flux would lead to an increase in wedge size and in erosion rate, processes that our kinematic model does not include. The increase in erosion rate could also be associated with 550 a direct, externally-driven increase in surface erosion rate associated with Quaternary climate change. Although the Apennines were not significantly affected by alpine glaciation, the cooling and strong cyclicity of the Quaternary climate may have led to an increase in erosion rate through the efficiency of periglacial processes and hillslope processes such as landsliding (Amorosi et al., 1996;Borgatti and Soldati, 2010;Simoni et al., 2013;Wegmann and Pazzaglia, 2009).

Conclusion 555
We present evidence from multiple thermochronometers that the spatial and temporal pattern of erosion rates in the Northern Apennines orogen differs at the regional scale. Time-averaged erosion rates from individual thermochronometers predict faster erosion rates derived from AHe ages on the Adriatic side relative to the Ligurian side. These results are consistent with erosion rates derived from paired AFT-AHe thermochronometer samples, which illustrate an increase in erosion rates through time on the Adriatic side, but a decrease in erosion rates through time on the Ligurian side. The pattern of erosion rates across the 560 orogen is also consistent with a kinematic model for an asymmetric orogen that includes both frontal accretion and underplating modes of crustal accretion, a slab rollback rate of 10 km/My, and prowedge erosion rates that are a factor of two higher than retrowedge erosion rates. This model suggests that that observed decelerations on the retrowedge are the result of the spatial advection of rock to the SW, although the observed acceleration of erosion rates on the prowedge requires external forcing, either through an increase in accretionary flux or through more erosive conditions linked to climate change.
All data used in this study are included in the text. Related codes for the erosion rate analysis and kinematic model are available upon request to the corresponding author.

570
Author Contribution EDE, MGF and SDW collected the detrital AFT samples in Northern Apennine catchments. EDE and MGF processed and analyzed the detrital AFT samples. SDW wrote the kinematic model. EDE, MGF, and SDW all contributed to the interpretation of the data. EDE wrote the manuscript and created the figures, with input from MGF and SDW.