An analytical solution for the exhumation of an orogenic wedge and a comparison with thermochronology data

An analytical solution for the exhumation of an orogenic wedge and a comparison with thermochronology data Elco Luijendijk1,2, Leo Benard3, Sarah Louis1,2, Christoph von Hagke2,4, and Jonas Kley1 1Department of structural geology and geodynamics, University of Göttingen, Goldschmidtstrasse 3, 37077, Göttingen, Germany 2Geological Institute, RWTH Aachen University, Wüllnerstr. 2, 52062, Aachen, Germany 3Mathematisches Institut, University of Göttingen, Bunsenstraße 3-5, 37073, Göttingen, Germany 4Department of Geography and Geology, University of Salzburg, Hellbrunnerstraße 34/III, 5020, Salzburg, Austria, Correspondence: Elco Luijendijk (elco.luijendijk@posteo.net)


Introduction
Thermochronology has increased our knowledge of the thermal and exhumation history of mountain belts (Reiners and Brandon, 2006). Numerical models that include faulting, erosion and heat flow enable using thermochronology data to quantify deformation (ter Voorde et al., 2004;Lock and Willett, 2008;Braun et al., 2012), and have for example been used to quantify 15 the geometry and slip rates of major faults (Herman et al., 2010;Robert et al., 2011). The combination of numerical models and deformation histories derived from balanced cross-sections has enabled more detailed reconstruction of deformation history (Erdös et al., 2014;McQuarrie and Ehlers, 2015). One drawback of using numerical models for the interpretation of thermochronometers in orogens is that they require significant resources in terms of computational cost, geological data and research time. Furthermore these approaches rely on the knowledge of the location of faults, which may not be complete and 20 depends on the quality of geological mapping and exposure. In addition, the models often assume that the fault blocks in between behave as rigid blocks, and may not include internal deformation of faults blocks by small-scale faults and folds inside these blocks.
Analytical solutions could offer an alternative as they require few computational resources and provide insight into the relation between input parameters and results without having to resort to computationally expensive sensitivity analyses. However, 25 to our knowledge no analytical solutions exist for the relation between deformation and exhumation of mountain belts. Interestingly, the problem of quantifying exhumation is somewhat analogous to the quantification of groundwater age, which deals with groundwater moving downward after infiltration instead of rock particles moving up to the surface. Analytical solutions for groundwater age have been published (Vogel, 1967), which suggests that the description of exhumation may be mathematically tractable. 30 Here we aim to develop a new closed-form analytical solution for the exhumation of an orogenic wedge that is in steadystate. The solution aims to include the main deformation mechanisms of orogenic wedges, internal deformation of the wedge, transport along a basal detachment, basal accretion and frontal accretion. In contrast to most numerical approaches we aim to represent a relatively simple system in which the orogenic wedge undergoes pervasive internal deformation, and we do not resolve individual faults. We evaluate how well this solution can represent thermal and deformation history by applying the 35 equation to a a well-studied cross-section in the Himalayas (Long et al., 2012;Coutand et al., 2014;McQuarrie and Ehlers, 2015), and by using the calculated exhumation histories to calculate thermochonometer ages and compare these to published data.

Conceptual model
Previous research has shown that at a first order the cross-sectional geometry of the brittle part of many orogens (i.e. until 40 the basal detachment becomes aseismic, typically at 300°C) can be approximated by a critical wedge (Buiter, 2012). After some amount of convergence the orogenic wedge reaches a steady-state that is formulated by the critical wedge theory (Davis et al., 1983;Chapple, 1978;Dahlen et al., 1984), whereby the topography of the wedge is a function of balance between the compressional force, the friction along a basal detachment and the gravitational force. Steady-state is maintained by the incorporation of new material by frontal (Platt, 1986), and basal accretion (Van Gool and Cawood, 1994), and the loss of 45 material by erosion at the land surface. In addition the rocks in the wedge itself deform as a result of compressional forces.
Here we derive a model that describes the motion of rock particles inside an orogenic wedge that is in steady-state, i.e., the geometry of the wedge is fixed and the flow of material into the wedge is balanced by erosion. The model covers three mechanisms for the deformation of an orogenic wedge: 1) Transport of the wedge along a basal detachment, 2) internal deformation of the wedge and 3) frontal and basal accretion. These three processes are shown schematically in Figure 1. The 50 geometry of the wedge and the parameters that are used in this manuscript are shown in Figure 2. The symbols in the figure are defined in Table 1.
Here we briefly summarize the way these processes are represented in the analytical solution that we derive below. The motion of material inside the wedge is described using the tip of the wedge as a fixed reference point.
1. Transport of the wedge along a basal detachment results in a uniform motion of the material in the wedge upward and 55 towards the tip of the wedge, with velocity vectors parallel to the wedge. Because the shape of the wedge is kept constant, transport along the detachment means that material is added on the right hand side of the wedge, and removed at the same rate at the top of the wedge. The material added at the right hand side and removed from the top of the wedge due to detachment transport are shown by blue and yellow areas respectively in Fig. 2a. 2. The internal deformation of the wedge is modelled as a uniform compression of the material inside the wedge. This is 60 the most simple way to represent the internal deformation of the wedge that in reality is accommodated by faulting and folding inside the wedge. The equations are based on the assumption that each rock particle in the wedge is compressed horizontally by exactly the same amount. This causes an expansion by the same amount for each particle. Note that we assume here that there is no volume loss. These assumptions lead to a total amount of expansion in the wedge that increases linearly from 0 at the tip of the wedge to a maximum at the right hand-side of the wedge. The linear correlation 65 between wedge thickness and exhumation is the result of more rocks being available to compress in a thicker wedge. The shape of the wedge remains constant, and therefore the deformation caused by compression is compensated by adding material on the right hand side of the wedge and by eroding an equal amount of material from the top of the wedge (Fig.   2b).
3. Frontal and basal accretion add material to the front and the base of the wedge, as shown in Fig 2c and Fig. 2d, tively. Given that we use the fixed geometry of the wedge as a reference frame, this causes a movement of rock particles inside the wedge to the right and upward. Note that by itself frontal accretion as described in our simple conceptual model will generate an unrealistic addition of material to the wedge above the land surface (Fig. 2c). Therefore frontal accretion should only be used in combination with rates of detachment transport, compression or basal accretion that are sufficient to balance out the material added with erosion at the land surface.

Development of analytical solution
We first develop solutions for the velocity of rock particles in the wedge (section 3.1), and subsequently use these equations to develop expressions for the trajectories of rock particles (section 3.2). We initially develop separate solutions for the transport of the wedge and basal and frontal accretions, and a solution for the velocity field due to compression of the wedge. The latter  Table 1.  and finally a full analytical solution for particle trajectories inside an orogenic wedge resulting from compression, transport along a detachment and accretion, which is given by equations 10 and 12. The full derivation of the solutions is provided in appendix A for the equations of particle velocity and appendix B for the equations for particle trajectories. The solutions are demonstrated for two example cases, which represent a wedge undergoing compression and transport along a basal detachment and a wedge undergoing compression and basal and frontal accretion. The parameters for these example cases are shown in 85   Table 2. The solution is also applied to a case study in the Himalayas, which is discussed in section 5.

Rock particle velocity
The combined horizontal and vertical particle velocity due to transport along a detachment and accretion is denoted as v xt and v yt , respectively: The horizontal and vertical components of the rock particle velocity due to compression are given by: and where and ζ are defined as: and v n denotes normalized compression velocity, which is defined as: The total velocity of rock particles can be found by combining the velocity due to transport (eqs. 1 and 2) and the velocity due to compression (eqs. 3 and 4). For the horizontal component of velocity this yields: And for the vertical component of velocity this yields: The combined velocity fields are shown in figure 3 and Fig. 4 for an example orogenic wedge. The parameters are shown in Table 2. The first case shows a wedge undergoing compression and transport along a detachment, which results in an overall rock particle motion towards the tip of the wedge and upward. Fig. 4 shows the velocity in a wedge undergoing compression and accretion, which results in an overall particle motion towards the interior of the wedge and upward. Figure 3. Particle velocity inside a wedge undergoing internal deformation and transport along a basal detachment, using parameter set 1 in Table 2. Panel a hows the velocity field resulting from compression of the wedge, panel b shows the velocity field resulting from transport of particles inside the wedge and panel c shows the velocity field due to both compression and transport.  Table 2.
Panel a hows the velocity field resulting from compression of the wedge, panel b shows the velocity field resulting from transport of particles inside the wedge and panel c shows the velocity field due to both compression and transport.

Particle trajectories 110
The equation for the horizontal position of particles over time follows from integration of the velocity equations as: Where x 0 is the horizontal position of the rock particle at t = 0.
The first solution for the vertical position is a simplified solution, where we take the depth below the surface for compression and for transport and add these two to get an expression for depth over time. This results in the following approximate solution 115 for the depth of rock particles (d) over time: Combining the expression for x (eq. 10) and the vertical velocity (eq. 9) and integration of the equation yields the full solution of the vertical particle position in the wedge: and 2) and the geometry parameters (eq. 5) and ζ (eq. 6), which are functions of α and β.

125
This equation can be adjusted to yield an expression for the depth of rock particles below the land surface (d), instead of the vertical position:

Validation of analytical solutions
The simplified analytical solution for rock particle trajectories (eqs. 10 and 11) and the full solution (eq. 10 and 12) were 130 verified using numerical particle tracking using the velocity fields described by the equations for particle velocity (eqs. 8 and 9). A number of particles were placed at the surface at t = 0. The velocity of these particles were calculated using these equations. In turn the velocities were used to calculate new particle positions at t − ∆t where ∆t is a time increment and the procedure was repeated. This procedure provides a history of depth versus time for each particle that can be used to verify the analytical solutions for particle depth over time, i.e. to explore whether the solution is consistent with the velocity fields that 135 were derived in section 3.1.
The comparison between the numerical model and the simplified and full analytical solutions shows that the full solution matches the numerical results perfectly for example cases dominated by transport along a basal detachment (Figs. 5) and basal and frontal accretion (Fig. 6). The maximum difference of the full solution and the numerical results is 0.12 %. The simplified solution deviates from the numerical and the full solution for both cases. The maximum difference with the full solution is 30 140 %. The error increases with depth ( Fig. 5d, 6d). This means that the simplified solution may still be useful to predict particle trajectories for the more shallow and low-temperature domains of the wedge. However, this is not the case for samples that are close to the tip of the wedge for models with significant basal and frontal accretion (Fig. 6d), where the deviation from the full solution is relatively large. 10 and 11), panel c shows the full solution (eq. 10 and 12), and panel d shows the difference between the simplified and full solution. The particle trajectories were calculated using parameter set 1 in Table 2. and 11), panel c shows the full solution (eq. 10 and 12), and panel d shows the difference between the simplified and full solution. The particle trajectories were calculated using parameter set 2 in Table 2.
5 Application to the exhumation of a mountain belt

Introduction
We used equations 10 and 12 to calculate particle trajectories for a case study based on a published cross-section in the Himalayas, the Kuru-Chu cross-section (Long et al., 2012;Coutand et al., 2014;McQuarrie and Ehlers, 2015). The crosssection is located in Bhutan in the central Himalayas and runs perpendicular to the strike of the Himalayas and the main faults.
The cross-section and the wedge geometry that was adopted based on this cross-section are shown in Fig. 7. The cross-section 150 was chosen because of the abundance of thermochronology data that covers both low-temperature (apatite fission track) and high temperature (muscovite Ar-Ar) thermochronology. In addition, the high convergence rate in the Himalayas makes it likely that exhumation is governed by the processes that are represented in the equations that we derived. Comparison of long-term (i.e. million year), intermediate (100-10 kyrs) and short term (yrs) erosion rate estimates indicates that the study area is in steady state (Ganti et al., 2016). This implies that other processes such as glacial erosion and isostatic compensation that are 155 not captured by the equations are likely to be relatively minor. Note that the goal here is not to provide a detailed geological case study, but to demonstrate the use of the equation to calculate deformation and exhumation rates.

Wedge geometry
The tip of the wedge was based on the location of the southernmost thrust fault in this part of the Himalayas, the main frontal thrust. The geometry of the wedge is based on published geological cross-sections (Long et al., 2012;Coutand et al., 2014;160 McQuarrie and Ehlers, 2015), with a surface slope (α) of 0.034 m m 1 , an average slope of the basal detachment (β) of -0.17 m m 1 and a total length of 200 km.
Particle trajectories were converted to thermal history using a uniform geothermal gradient of 15 C km −1 , which is based on a gradient of 10 to 20 C km −1 reported by Coutand et al. (2014). Surface temperature was calculated using a surface temperature at sea level of 24°C and a atmospheric lapse rate of 7°C km −1 (Coutand et al., 2014). Cooling ages were calculated using the thermal history and assumed resetting temperatures of 110°C, 180°C and 325°C for the apatite fission 170 track, for zircon (U-Th)/He and muscovite Ar-Ar thermochronometers, respectively (McQuarrie and Ehlers, 2018). Note that we aim for a simple model setup and do not include more complex calculation based either on a more formal definition of cooling ages (Dodson, 1973;Fox et al., 2014) or fission track annealing models or diffusion models of He or Ar, which we reserve for follow up work. Ages for samples that originated in the foreland of the wedge and that did not reach resetting temperatures in the wedge were calculated using an assumed exhumation rate of 1 × 10 −4 m a −1 . This was done to give these samples arbitrary high ages that approximate the typically higher ages found in foreland basins of orogenic wedges. Note that we did not calculate ages for samples that originated in the hinterland of the wedge, because it is difficult to make assumptions about their thermal history and initial age prior to entering the wedge.

Model calibration and parameter space exploration
The optimal values of the model parameters v c , v d , v xa , v ya were found using the downhill simplex algorithm (Nelder and 180 Mead, 1965), as implemented in the Python module scipy (Virtanen et al., 2020), to minimize the mean absolute error (MAE) of the calculated thermochronometer ages.
In addition to model calibration we also ran a set of models with different parameter combinations to explore the sensitivity of

Results
The the best-fit parameter values show a relatively good fit of the calculated thermochronometer ages to the measured ages in the cross-section (Fig. 8a), with a coefficient of determination of 0.75. The calibrated parameters values indicate that transport along the basal detachment and internal deformation are equally important mechanisms driving exhumation of the wedge, with a calibrated value of the transport (v d ) and compression velocity (v c ) of 5.0 × 10 −3 m a −1 and 4.8 × 10 −3 m a −1 , respectively (Table 3). Frontal and basal accretion play an insignificant role in this cross-section, with a best-fit value of 1.7 × 10 −6 m a −1 for the horizontal component (v xa ) of the accretion vector, and a value of 0.0 for the vertical component (v ya ).
The exploration of parameter space show a relatively small area with a good fit (Fig. 9). A reduction in detachment transport 195 velocity can to some extent be taken up by an increase in compression velocity and vice versa. However, both processes results in fundamentally different exhumation and cooling age curves, with a straight line for models dominated by transport over a detachment only (Fig. 10a, b), and a curved line for models that include compression (Fig. 10c, d). This means that the degree of curvature of cooling ages along a cross-section and the rate of decrease in ages towards the interior of the wedge can be used to quantify the rate of internal deformation of orogenic wedges.

200
Note that Coutand et al. (2014) report that the geothermal gradient varies systematically in their model experiments from 10 C km −1 at the tip of the wedge to 20 C km −1 in the hinterland. This is caused by high advection rates in the hinterland and the cooling effects of the subducting footwall plate in the foreland. Our results do not include a spatially varying geothermal gradient, and could therefore have underestimated geothermal gradients and overestimated thermochronometer ages for the samples furthest away from the wedge. However, the magnitude of this effect is difficult to constrain given the lack of 205 subsurface temperature or heat flow data in this part of the Himalayas. (2015) shows a better fit of the analytical solution to the data (Fig. 11). Note that some of the misfit of the McQuarrie and Ehlers (2015) data may be due to errors introduced by digitizing the model results from their publication. Nonetheless the comparison shows that the 210 new analytical solution fits as well or better than a more sophisticated numerical model that included the activity of individual thrust faults based on balanced cross-sections. This suggests that, for this cross-section, deformation of the wedge can be better represented by uniform deformation.

Comparison of our results with published thermochronometer model ages by McQuarrie and Ehlers
The sum of the best-fit values of internal deformation and transport along the detachment results in a total shortening rate of 10 mm a −1 . This is comparable to previous findings. Inverse numerical models by Coutand et al. (2014) found the best fit 215 an overall convergence rate of 22 mm a −1 . Roughly half the convergence rate is accommodated by underthrusting, and the convergence taken up by the wedge itself was found to be 10 to 11 mm a −1 , which is equal to our estimate. McQuarrie and Ehlers (2015) found slightly lower total convergence rate of 18 mm a −1 . However, the rate of over and underthrusting were not reported separately. Coutand et al. (2014) and McQuarrie and Ehlers (2015) also included models with variable convergence rates over time. However, in Coutand et al. (2014) the reported misfit for models with variable convergence rates was similar to 220 models with steady convergence rates. McQuarrie and Ehlers (2015) reported a highly variable convergence rate through time.
However, the extent to which these inferred changes in convergence rate improved the model fit statistics was not reported.
In contrast, the good fit of the thermochronometer ages predicted by the steady-state deformation and exhumation equation (Fig. 8) suggests that instead this part of the Himalayas may be in steady-state. This is corroborated by the finding that erosion rates determined on shorter averaging timescales reproduce thermochronological data within error, as expected for steady state 225 landscapes (Ganti et al., 2016). The distributed nature of internal shortening of the wedge is supported by the deep exhumation of Himalayas southern slope as demonstrated by erosional remnants of the once continuous MCT thrust sheet overlying rocks of the Lesser Himalaya (Yin, 2006). This supports the idea that, over long timescales, internal shortening is distributed in nature and does not only occur by large offsets on a few major faults. Panel a shows calculated and published apatite fission track (AFT), zircon (U-Th)/He (ZHe) and muscovite Ar-Ar ages (MAr), and the exhumation rate at the surface. Panel b shows the particle trajectories that were used to calculate the cooling ages, along with the isochrons.
The calibrated parameter values used for this model run are shown in Table 3. The parameters for these model runs can be found in Table 3.

230
We have developed an analytical solution for rock particles trajectories in an orogenic wedge that is in steady-state and that undergoes uniform internal deformation, transport along a basal detachment, frontal accretion and basal accretion. The solution was validated using a numerical solution based on calculated velocity fields and numerical particle tracking. We also derived a simplified solution that provides a good approximation for the upper half of the orogenic wedge, and may therefore be useful to quantify exhumation in the shallow and low-temperature part of orogenic wedges. The solution predicts that exhumation rates 235 follow a curved line along the wedge, with an increase in exhumation rate and a decrease in thermochronometer ages towards the interior of the wedge. The degree of curvature depends on the rate of internal deformation of the wedge. Application of the equation to a published cross-section in the Himalayas show that the equation provides a very good fit to thermochronology data, with an coefficient of determination (R 2 ) of 0.75. This indicates that the Himalayas may be in steady-state and that, at a large scale, the exhumation of mature mountain belts may be approximated by a relatively simple model of uniform and 240 steady-state deformation, accretion and transport. The equations enable a direct estimate of deformation and shortening rates from thermochronology data in orogens.
Code availability. Jupyter notebooks that reproduce the workflow and the figures presented in this study have been published at Zenodo (Luijendijk, 2021) and are available at GitHub (https://github.com/ElcoLuijendijk/wedgex). The code is distributed under the GNU General Public License, version 3. The notebook collection include a notebook that reproduces the figures that compare the equation to numerical 245 results (Fig. 5, 6) and a notebook that can be used to calibrate parameters and compare the results to thermochronometers (Fig. 8, 9). Wedgex depends on the Python modules Numpy (Harris et al., 2020), Scipy (Virtanen et al., 2020), Matplotlib (Hunter, 2007), Pandas (McKinney, 2010;Reback et al., 2021), scikit-learn (Pedregosa et al., 2011) and Scientific colour maps (Crameri, 2021).
Appendix A: Derivation of equations for rock particle velocity A1 Velocity due to wedge transport and accretion 250 For transport of the wedge along the basal detachment the additional horizontal velocity is constant everywhere and defined as v d . The vertical transport velocity is the product of the slope of the basal detachment and the horizontal transport velocity: Exhumation due to basal and lateral accretion is described by a constant term for the horizontal and vertical velocity of rock particles inside the wedge that represents the lateral and upwards movement of rock particles caused by accretion. These are 255 denoted by v xa and v ya , respectively. The combined horizontal and vertical particle velocity due to transport along a detachment and accretion is denoted as v xt and v yt , respectively: A2 Velocity due to internal deformation 260

A2.1 Horizontal velocity
The horizontal velocity due to internal deformation is represented by uniform compression that is assumed to increase linearly from 0 at the left-hand side to v x = v c at the right hand side. This yields the following expression for velocity: A2.2 Vertical velocity at the top of the wedge 265 For compression without any volume loss the volume balance of the wedge as shown in Figure 2 can be written as: The velocity at the top of the wedge is assumed to increase linearly from v ys = 0 at the left hand side to v ys = v ym at the right hand side (x = L): The area that the wedge extends at the top should be equal to the area that is compressed at the right hand side: which can be combined with the previous eq. to form: which can be reorganized as to write yield an expression of the particle velocity at the top of the wedge as a function of x, 275 the height of the wedge (H), the length of the wedge (L) and the compression velocity (v c ):

A2.3 Vertical velocity inside the wedge
We assume that the vertical velocity due to wedge compression (v yc ) at any point (x,y) increases linearly from v b at the bottom of the wedge (y=b) to v ys at the top of the wedge: which can be rearranged as: where f is the relative vertical position of the rock particle in the wedge and is defined as The equation for v ys that was derived previously is: combining these two equations yields: which after reorganizing becomes: reinserting the original expression for f = (y − b)/h results in: and inserting b = βx and h = (α − β)x = γx and v yb = βv xb = β x L v c yields: which after rearranging yields an expression of the vertical particle velocity due to compression: We simplify this equation by introducing two new terms for the geometry of the wedge: and introducing the normalized compression velocity v n : Inserting these three new terms into the equation yields: A3 Velocity due to compression, transport and accretion The total velocity of rock particles can be found by combining the velocity due to compression (eqs. A4 and A22) and the velocity due to transport (eqs. A2 and A3). For the horizontal velocity this yields: And for the vertical component of velocity this yields: Appendix B: Particle trajectories B1 Trajectories due to compression B1.1 Horizontal particle position We integrate the equation for horizontal velocity due to compression of the wedge (eq. A4) to yield an expression for the 315 horizontal particle position over time. The equation for horizontal velocity can be rewritten as: Integration of this equation with x = x 0 at t = 0 yields:

320
Taking the equation for vertical velocity due to compression (eq. A22) and replacing v yc = ∂y ∂t results in: This can be simplified as where the new constants C 1 and C 2 are defined as: and Integration of this equation with boundary condition y = y 0 at t = 0, yields an expression for y: y = ζv n v n − v n x 0 e vnt + y 0 − ζv n v n − v n x 0 e vnt (B7) 330 which simplifies to: The term ζ 1 − equals β, which yields: and taking y 0 = αx 0 results in: 335 y = βx 0 e vnt + (α − β) x 0 e vnt (B10) B2 Trajectories due to transport, compression and accretion

Horizontal position
Using the equation for the horizontal velocity of rock particles (eq. A23) and recasting the velocity term as the derivative of position over time results in: Integration of this equation yields: Where K is a constant. Taking x = x 0 at t = 0 gives a solution for the horizontal particle position over time:

B2.1 Simplified solution for vertical position
The first solution for the vertical position is a simplified solution, where we take the depth below the surface for compression and for transport and add these two to get an expression for depth over time. For deformation due to compression the vertical position of rock particles over time is given by equation B10. If we define depth due to compression as we can combine this with eq. B10 to yield: which can be rearranged as: The equation for depth over time due to the transport of rock particles in the wedge follows: Integration of this equation yields: rewriting the equation in terms of depth due to transport (d t ) gives: Combining the depth due to transport and due to compression results in an approximate solution for the depth of rock particles (d) over time: x 0 e vnt − e vnt (B20)

B2.2 Full solution for vertical position
Inserting the expression for x derived earlier (eq. B13) into the solution for the vertical velocity (eq. A24) yields: which can be shortened to ∂y ∂t = C 1 y + C 3 e vnt + C 4 (B22) where which can be recast as: The solution of this integral yields: where K is a constant and λ(t) is: For the boundary condition y = y 0 at t = 0 this yields: Reinserting the expressions for the constants C 1 , C 2 , C 3 and C 4 yields the following equation: well. EL and LB derived the equations. EL led the model experiments and the writing. CvH and SL drafted figure 1 and 2. All authors contributed to the discussion of the results and writing the manuscript Competing interests. The authors declare that no competing interest are present.