Insight into collision zone dynamics from topography: numerical modelling results and observations

Abstract. Dynamic models of subduction and continental collision are used to predict dynamic topography changes on the overriding plate. The modelling results show a distinct evolution of topography on the overriding plate, during subduction, continental collision and slab break-off. A prominent topographic feature is a temporary (few Myrs) basin on the overriding plate after initial collision. This "collisional mantle dynamic basin" (CMDB) is caused by slab steepening drawing, material away from the base of the overriding plate. Also, during this initial collision phase, surface uplift is predicted on the overriding plate between the suture zone and the CMDB, due to the subduction of buoyant continental material and its isostatic compensation. After slab detachment, redistribution of stresses and underplating of the overriding plate cause the uplift to spread further into the overriding plate. This topographic evolution fits the stratigraphy found on the overriding plate of the Arabia-Eurasia collision zone in Iran and south east Turkey. The sedimentary record from the overriding plate contains Upper Oligocene-Lower Miocene marine carbonates deposited between terrestrial clastic sedimentary rocks, in units such as the Qom Formation and its lateral equivalents. This stratigraphy shows that during the Late Oligocene–Early Miocene the surface of the overriding plate sank below sea level before rising back above sea level, without major compressional deformation recorded in the same area. Our modelled topography changes fit well with this observed uplift and subsidence.


Introduction
In this study we aim to look at the evolution through time of topography on the overriding plate at a collision zone. 2-D numerical models of lithosphere-mantle interactions at subduction and continental collision zones are used. These generic modelling results are compared to specific observations from the Arabia-Eurasia collision zone. The study aims to illustrate that topography can be used as an indicator of the dynamics associated with continental collision and slab break-off.
To a first order, the topography of the Earth's surface can be explained by isostasy causing regions of variable crustal thickness to equilibrate at different heights. These variation in crustal thickness are caused by a variety of local processes including tectonic thrusting and rifting. However crustal thickness is just one contributor to the Earth's topography and to fully account for topography we need to consider all forces that act on the Earth's surface. The forces that influence topography can be broadly categorised into isostatic forces, from thermal and crustal buoyancy, flexure forces associated with the lithosphere, crustal strength as well as stresses imposed at the base of lithosphere due to mantle dynamics. Topography driven by flow in the mantle is referred to as dynamic topography (Lithgow-Bertelloni and Silver, 1998). Many modelling studies have found a change in the dynamics of mantle flow during the transition between subduction, collision, and oceanic slab break-off (Gerya et al., 2004;Andrews and Billen, 2009;Duretz et al., 2011). These changes in flow in the mantle would be expected to affect the topography generated due to changes in the stresses at the base of the lithosphere (Faccenna and Becker, 2010). Subduction zone topography is characterised by a long linear oceanic trench flanked by an outer rise on the subducting plate side and raised topography on the overriding plate (Melosh and Raefsky, 1980;Hager, 1984;Gephart, 1994). The origin of this topography is the result of the sum of dynamic forces and isostatic forces (Forte et al., 2010;Husson et al., 2011). Linking of subduction dynamics to topography (Husson, 2006) has shown how features such as the back arc basin depth can be correlated to subduction velocity. Numerical modelling (Husson, 2006) and analogue models (Husson et al., 2011) further show predictions of topography due to mantle dynamics.

Published by Copernicus
Crustal shortening and thickening at collision zones produces topographic signals, via isostasy, in addition to those caused by mantle processes. The transformation of a subduction zone to a fully formed collision zone is thought to be categorised by two important events: initial collision of the continental material and later, break-off of the subducting oceanic slab. Slab break-off has been proposed to occur when the arrival of continental material at the subduction zone slows and eventually stops subduction (Davies and Blanckenburg, 1995). This resistance to continued subduction by the buoyant continental crust could then cause a build up of stress where the oceanic material joins the continent. The slowdown of subduction is also thought to allow the subducting slab to heat up, further weakening it. A combination of the build up in stress and thermal weakening is thought to cause the slab to neck and detach. Recent numerical modelling work (Gerya et al., 2004;Andrews and Billen, 2009;van Hunen and Allen, 2011) has shown that both these mechanisms play important roles in the ability of slabs to detach. Estimates for the lag time of break-off after initial collision vary between 5-40 Myrs (Andrews and Billen, 2009;Gerya et al., 2004;van Hunen and Allen, 2011) and have been found to be dependent on the age and associated strength of the subducted lithosphere.
Many studies have produced numerical models of the dynamics of slab break-off (Buiter et al., 2002;De Franco et al., 2008;Faccenda et al., 2009;Andrews and Billen, 2009;van Hunen and Allen, 2011). Various estimates have been made of the change in topography due to slab break-off, ranging from 1-10 km (Gerya et al., 2004;Andrews and Billen, 2009;Duretz et al., 2011) of uplift due to the loss of the slab. This large variation in estimates clearly shows the need for future work.
One of the challenges with topography is to separate the contributions to topography from different sources. Collision zones produce elevated topography due to the shortening and thickening of the crust as well as uplift associated with slab break-off, plate flexure and dynamic effects from the mantle. Recent work (Duretz et al., 2011) has also highlighted that different break-off dynamics would produce different topography changes. Numerical modelling results have also been related to observation for South America (Shephard et al., 2010) showing that the modelled topography changes can be linked to observations. The study by Shephard et al. (2010) focuses on the proposed drainage reversal in South America during the opening of the south Atlantic and proposes that the subsidence and uplift patterns required can only be explained by changes in subduction dynamics influencing the topography.
The Arabia-Eurasia collision zone ( Fig. 1) offers an area for the study of topography changes associated with collision and slab break-off (Agard et al., 2011). Collision occurred here after the closure of the Neo-Tethys ocean basin. Collision is still active in the region with GPS measurements putting the current north-south convergence rate at ∼2.5 cm yr −1 (Sella et al., 2002). The time of initial collision is still debated, but estimates vary between >65-5 Ma (Berberian and King, 1981;Philip et al., 1989;McQuarrie, 2003;Ghasemi and Talbot, 2006). A common estimate is late Eocene ∼35 Ma (Vincent et al., 2007;Allen and Armstrong, 2008;Ballato et al., 2011;Mouthereau et al., 2012;Agard et al., 2005) based on deformation on both sides of the suture at that time, and a shutdown of arc magmatism. In this study we will use the estimate of 35 Ma for the initial collision in this region.
Local slab break-off has been estimated to have occurred at 10 Ma (Ghasemi and Talbot, 2006;Omrani et al., 2008) for the Arabia-Eurasia collision giving a delay time of ∼25 Myrs since initial collision. This estimate for the timing of slab break-off comes from observations of collisional magmatism (Omrani et al., 2008) which is used as an indicator of slab break-off. In that study, collisional magmatism is proposed to be produced by slab break-off due to the descent of the slab into the mantle draws hot material into the mantle wedge region. There is also evidence from mantle tomography that slab break-off has occurred at the Arabia-Eurasia collision zone (Lei and Zhao, 2007), which shows low velocity regions where the slab would be expected. These low velocity regions can be interpreted as areas where the slab is no longer present and has been replaced by hot mantle material. Neo-Tethys opened in the Permian (Ş engör et al., 1988) making the subducted oceanic plate entering the subduction zone just before collision approximately 200 Myrs old. This old oceanic crust will have a thick lithosphere and high strength, so slab break-off would be expected to take longer than for young weak slabs. These properties support the case for using the upper estimate for the delay time between collision and break-off, of around 25 Myrs, as shown in numerical modelling studies (Andrews and Billen, 2009;van Hunen and Allen, 2011). The current and past topography of the Arabia-Eurasia collision zone shows some interesting features that may allow further understanding of the dynamics of the transition from a subduction zone to a collision zone. Much of central Iran, north of the Zagros suture and within the present Turkish-Iranian plateau ( Fig. 1) also has the advantage of being relatively unaffected by compressional deformation for approximately 10-15 million years after initial collision (Ballato et al., 2011;Mouthereau et al., 2012). Instead, the region underwent marine carbonate-dominated sedimentation, without major fault control (Morley et al., 2009). This history means that possible dynamic mantle effects on topography are expressed without an overprint of crustal shortening and thickening. This in turn will allow direct comparison of modelled dynamic topography with topography inferred from the sedimentary record in the region.

Methodology
The numerical modelling is done with a finite element geodynamical code, Citcom (Moresi and Gurnis, 1996;Zhong et al., 2000;van Hunen and Allen, 2011). Citcom uses a cartesian grid, assumes incompressible flow and makes the Boussinesq approximation. Non-dimensional governing equations are as follows (van Hunen and Allen, 2011): where symbols are defined in Table 1 and mantle temperature T m = 1350°C (mantle temperature). The code used solves for conservation of mass, momentum, energy and composition (van Hunen and Allen, 2011), using quadrilateral finite elements with bi-linear velocity and constant pressure.

Model setup
The models in this study simulate the closure of a small oceanic basin leading to continental collision and subsequent slab detachment. The initial model setup is shown in Fig. 2. The modelling domain is 660 km by 2640 km giving a 1:4 aspect ratio. The grid is refined in the top 200 km and between 1700 km and 2200 km in the horizontal. This gives a grid resolution, over the collision zone, of 4 km by 5 km, in the x and z directions, respectively. Boundary conditions for the model are free slip on the top and sides with no slip on the base. The no-slip boundary condition at the base of the model is designed to simulate the interaction of a slab with a higher viscosity layer such as that proposed for the phase change at 660 km. The thermal boundary conditions are 0°C at the surface and mantle temperature (1350°C) at the bottom and left edge. The right edge has a zero heat flux boundary condition. Subduction is initiated by a hanging slab and facilitated by a zone of weak material between the subducting plate and the overriding plate ( Fig. 2; Appendix 1). This numerical 13.2 km wide weak zone is maintained in the same shape throughout the model run and kept at a constant mantle viscosity (1 × 10 20 Pa·s) to allow permanent decoupling of the two plates. The coupling of the two plates in a real subduction zone is likely to change over time. As the coupling evolution during continental collision is not well known, we have chosen to keep the weak zone constant over time. The model also contains a 100 km deep, 200 km wide mantle wedge that is also held constantly at mantle viscosity (1 × 10 20 Pa·s), thereby mimicking a weak, hydrated area above the subducting slab. These two features allow decoupling of the subducting and overriding plates without the complexities of slab dehydration and wedge hydration.  serves to decouple the plate from the left boundary to allow the subducting plate to move freely. The model setup initially has a 60 Myrs old oceanic lithosphere partially subducted under a continental overriding plate with a zero age mid-ocean ridge at the left edge. The initial thermal structure of the oceanic lithosphere is calculated using the half-space cooling model. This plate age is representing generally "old" oceanic lithosphere, since beyond 70 Myrs there is little change in plate thickness and "apparent thermal age" (Ritzwoller et al., 2004;van Hunen et al., 2005). The subducting plate has a 700 km long 40 km thick continental crustal block embedded in it. The overrid-ing plate also has a 40 km thick continental crust and is fixed to the right edge of the model. The thermal structure of the continental regions is set as a linear geotherm from 0°C at the surface to mantle temperature at 150 km depth. The age distribution of the oceanic lithosphere surrounding the continental block is linearly increasing with distance from the MOR, and does not reflect any rifting events that might have isolated the continental block. This simplification does not significantly affect the presented modelling results.

Solid Earth, 3, 387-399, 2012
www.solid-earth.net/3/387/2012/ Compositional variation (crust versus mantle material) is advected using particle tracers (Di Giuseppe et al., 2008). The continental crust in the model resist subduction due to its compositional buoyancy. Oceanic crustal buoyancy is ignored in the models, as the transformation of basalt to eclogite occurs at 30-40 km, which would remove most compositional buoyancy (Cloos, 1993), making the buoyancy of the subducted oceanic lithosphere an almost purely thermal effect.

Rheology
The strength of the material in our model is governed by temperature and stress dependent rheology. Four different deformation mechanisms are used (diffusion creep, dislocation creep, a stress limiting rheology, and a model maximum viscosity). For both diffusion creep and dislocation creep the viscosity is described by Rheological pre-exponents of 10 20 and 2.8 × 10 13 and rheological power law exponents of 1 and 3.5 are used for diffusion creep and dislocation creep, respectively (Karato and Wu, 1993). Adiabatic compression is ignored in the Boussinesq approximation which allows simplification of the model by using zero activation volume. Rocks have limited strength and yield under high-enough stresses, both close to the surface, as described by Byerlee's rule (Byerlee, 1978), and deeper using Peierl's creep (Guyot and Dorn, 1967;Kameyama et al., 1999). We impose the following stress limiting rheology to simulate both these mechanisms The effective viscosity for each element is simply defined as the minimum of the individual viscosity components. The same rheology is assumed for both crustal and mantle material. For a further description of rheology see van Hunen and Allen (2011).

Model Topography calculation
Model topography is calculated with and without an "elastic filter" to illustrate the effect of an elastic strength. The unfiltered topography is purely the subaerial isostatic response from the normal stress (q 0 ) at the top free-slip surface. For this case (without any elastic strength) the topography is calculated by assuming that the surface is in direct isostatic equilibrium and the normal stresses generated by the model are balanced by a column of crustal material. The height needed for this column at each nodal point in the model gives a first order estimate of the topography.
Although this does give an estimate of the topography due to both dynamic forces and the isostatic buoyancy of the material, it does not account for any elastic properties of the overriding plate. To account for the lateral elastic strength of the lithosphere, the normal stress is filtered using the flexure equation and symbols are defined in Table 1 for an elastic material (Turcotte and Schubert, 2002). Equation (7) is solved using a finite difference technique. The boundary conditions used in this are (1) the left edge of the model at -2.7 km to simulate the depth of the mid ocean ridge; (2) the right edge is at its expected isostatic height having corrected its height by the same amount required for boundary condition 1; (3) there is no change in topography gradient at either boundary, e.g. the second differential of topography with respect to distance is zero. An effective elastic thickness (Te) of 30 km was chosen as representative for the region based on elastic thickness estimates for the whole of Africa and the Middle East (Pérez-Gussinyé et al., 2009). We have also produced model results with 20 km, 40 km and 50 km elastic thicknesses (Appendix 2). Results for the pre-collisional back-arc depth basin depth are comparable to those produced by He (2012) for a visco-elastic numerical model of the mantle wedge and overriding plate. Instead of assuming two decoupled plates, we have chosen to use one elastic plate for the whole model: the topography of interest postdates initial collision, such that the two collided plates are assumed to behave mechanically as one. Although this is not likely to offer the best solution during on-going subduction, the coupling between the two plates will be complex and vary over time. Therefore, we choose to not add extra complexity to the model by trying to estimate the change in plate coupling over time.
This secondary calculation of elastic deformation offers an alternative to the use of a model with fully free surface and viscoelastic rheology. Models with a full free surface are computationally expensive (Schmeling et al., 2008) though do simulate true elastic deformation. Although the elasticity in our models has to be calculated as a secondary process we would argue that we still resolve the important topographic features.

Results
This section describes how model dynamics are reflected in the surface topography and how these topographic features evolve over time.
As subduction slows down after the onset of collision (Molnar and Stock, 2009), the subducted slab steepens before finally slab break-off occurs. Figure 3 shows four stages of the collision process, the resulting viscosity profile and the topographic expression. The two continental regions in the model are associated with relatively high elevation (Fig. 3a), created by the buoyancy of the continental blocks. The The features that are of particular of interest here are those that are on the overriding plate. This is because this area, in a collision setting, is most likely to preserve evidence of topography changes. The other reason to focus on the overriding plate is that convection in the mantle wedge is proposed to be responsible for dynamic topography. The results for this region of the model show a depression on the over-riding plate at 300 km from the trench both before and after initial continental collision, but not after break-off. Before collision this feature is dynamically produced by flow in the mantle wedge (Husson, 2006). This feature is purely dynamic since the model contains no mechanism for slab rollback and associated spreading and thinning of the overriding crust. At 7 Myrs (Fig. 3b), this collisional mantle dynamic basin (CMDB) has deepened during the start of collision. At 10 Myrs (Fig. 3c) it has become shallower again, until it has almost disappeared by the time slab break-off occurs at 17 Myrs (Fig. 3d). The viscosity plots also show evidence of the basin in the weakening of the overriding plate at the position of this dynamic basin. This weakening is due to the stress imposed on the plate, caused by mantle wedge material being drawn away from the base of the plate. This stress causes weakening due to the stress and temperature dependent rheology.
Post-collisional uplift also occurs on the overriding plate between the trench and CMDB. At 7 Myrs and 10 Myrs, this uplift is restricted to the region close to the trench where subduction of buoyant continental material is taking place. The uplift predicted in this area should be considered a maximum, since we do not model any imbrication mechanisms. This means that initially large quantities of buoyant continental material are subducted, whereas we might expect some of this material to be removed during the subduction process. This uplift is also coincides with the expected position of the arc (from the subduction angle) at 100 km behind the trench and 200 km in front of the CMDB.
After slab break-off at 17 Myrs the uplift moves further into the overriding plate. This pattern of topography change is further displayed in Fig. 4 by a topography-time map. This topography presented is elastically filtered with an effective elastic thickness (T e ) of 30 km. This filtered topography is the same as that shown in green in Fig. 3. The elastic filter has the effect of removing a lot of the short wavelength variation in these topographic features, and reduces the overall amplitude of the signal. Topography time plots are useful as they show which topographic features are persistent for a significant time. They also allow easy linking of the dynamics of the model to topographic changes. Figure 4a contains three highlighted areas of particular interest and also shows the subduction velocity over time. The subduction velocity during this initial onset of subduction varies over a large range. This is caused by the initial large slab pull and lack of interaction with the 660 km discontinuity (Quinquis et al, 2010). This initial high velocity should not affect the results, as subduction has slowed to a more geologically realistic rate by the time collision occurs. Region 1 shows a basin on the overriding plate, approximately 300 km back from the trench, that is present between initial collision and break-off. We will refer to this post collisional basin as a collisional mantle dynamic basin (CMDB). This region is due to slab steepening as it does not correlate with subduction velocity, unlike the back arc basin feature in the same position during on-going subduction ( Fig. 4a and d).
The effect of slab steepening can be more clearly observed by comparison between Fig. 3 panels A (5.5 Myrs), B (7 Myrs) and C (10 Myrs). These show how the slab, defined by the region of high viscosity, steepens during collision. The topography profiles for these time slices show an increase in depth of the CMDB. The steepening slab causes a deepening of the basin because it draws material away from the overriding plate causing the surface to sink. This hypothesis is further supported by the velocity fields for the three time slices in Fig. 3. At 7 Myrs the velocity field under the CMDB region has steepened to almost vertical position compared to the velocity field for on-going subduction at 5.5 Myrs. At 10 Myrs the velocity field is also vertical but reduced in magnitude. This reduction in magnitude of the velocity field between 7 Myrs and 10 Myrs explains the reduction in size of the CMDB. The combined effect of slab steepening over time is shown in Fig. 4a area 1. Region 2 in Fig. 4a marks the initial uplift after collision. This is due to the partial subduction of the continental block in the subducting plate. This creates a region with double continental buoyancy and hence of large uplift. Figure 3 shows the continental material (white contour line) in the subducting plate thrust under the continental material on the overriding plate. This partial subduction of continental material creates areas close to the trench on the overriding plate with crustal thickness of up to 80 km. This large crustal thickness creates uplift of 3 km after elastic filtering which is potentially unrealistically high. This should be treated as a maximum uplift value due to the lack of an imbrication mechanism in the model.
Region 3 in Fig. 4a shows migration of the uplift further into the overriding plate after slab break-off. This migration of uplift is also clear from comparison of Fig. 3 panel C (10 Myrs) and D (17 Myrs). This change in uplift pattern in the short-term is caused by the redistribution of stress after slab detachment allowing the whole region to uplift. The removal of the slab also allows some continental material to be educted back up the subduction channel as well as being underplated on to the overriding plate. The effect of this on the overriding plate is a slightly increased crustal thickness in the region due to under plating which accounts for the long-term uplift seen.
To test the sensitivity of model results to various subduction model parameters, a number of model calculations with different initial setups were performed. For comparison of these models, the maximum depth of the post collision basin was recorded. Appendix 1 elaborates how different input parameters affected the depth of the collisional mantle dynamic basin. It was found that neither the wedge width nor viscosity affected the ultimate depth of the post collision basin by a statistically significant amount. The depth of the mantle wedge and the weak zone viscosity were found to have an effect on the basin depth. This is because these start-up parameters affect the coupling between the two plates, which in turn affects the model dynamics. The model was also found to be sensitive to the elastic thickness (T e ) (explored further in Appendix 2) and surface yield strength. Although several model parameters potentially could affect the topography, the startup conditions chosen in our preferred model gives geologically reasonable subduction velocities once the slab interacts with the 660 km discontinuity.
We have also tested the model at double the resolution presented here. The result from these models qualitatively show very little difference to the standard model. Quantitatively, the maximum depth of the CMDB is positioned 16 km closer to the trench in the high resolution model compared to the standard resolution model. The CMDB maximum depth, in the high resolution model, was also found to reduce by 16 % compared to the standard resolution model. This shows that the features identified in this study are numerically robust.

Discussion
The modelled topography evolution in this study shows a clear pattern of topography change through the processes of subduction, collision and slab break-off. Identification of these topography changes in collision zones on Earth provides insight into continental collision dynamics. A schematic overview diagram of the main features of the topography change is shown in Fig. 5.
To illustrate the applicability of the presented model results, we compare them to geological observations from the Arabia-Eurasia collision zone. The Turkish-Iranian plateau is typically 1.5-2 km above sea level at present. To achieve such high topography, isostatic theory suggests the need for a large crustal thickness (Turcotte and Schubert, 2002). Part of the crustal thickening is caused by thrusting within the crust, but active thrusting seems limited to areas with elevations below 1250 m (Nissen et al., 2011). Wholesale underthrusting/subduction of the northern side of the Arabian plate beneath Eurasia has been imaged on deep seismic lines (Paul et al., 2006), which is a plausible mechanism for generating overall crustal thickening, and hence isostatic uplift, in regions with elevation above the limit of seismogenic thrusting. Our modelling highlights the potential contribution to crustal thickening, and hence topography of such underthrusting (Fig. 4a). The spatial extent of the plateau in the model output is more localised to the overriding plate compared to the much broader uplift observed in Turkey and Iran (see region 3, Fig. 4a). This is because of the complete Fig. 6. The extent of Upper Oligocene -Lower Miocene carbonate strata present on the overriding Eurasian plate in the Arabia-Eurasia collision zone. These carbonate deposits were laid down after initial collision. Compiled from (National Iranian Oil Company, 1977aNational Iranian Oil Company, 1977bNational Iranian Oil Company, 1978;Ş enel, 2002;Reuter et al., 2007;Morley et al., 2009). decoupling between the plates in our models after collision, whereas the coupling in a real subduction zone is likely to change throughout the subduction process. The uplift predicted by the model is larger than the observed uplift in the region because of complete subduction of buoyant crustal material and the lack of imbrication to remove some of this material.
A second topographic feature of the Arabia-Eurasia collision is related to a sequence of Upper Oligocene -Lower Miocene carbonate sedimentary rocks found on the overriding Eurasian plate, in modern day Iran and southern Turkey ( Fig. 6; Reuter et al., 2007). The unit is known in central Iran as the Qom Formation. The Qom Formation and its equivalents are typically 500-1000 m thick in central Iran (Gansser, 1955;Morley et al., 2009), and lie sandwiched between terrestrial clastic strata of the Lower Red Formation and the Upper Red Formation. This relationship implies that the area on the overriding plate was above sea level during the Early Oligocene, then subsided below sea level during the Late Oligocene -Early Miocene, before returning to above sea level in the late Early Miocene (Burdigalian stage). These carbonate deposits also extend laterally along most of the collision zone, suggesting that they are intimately associated with the collision process. Importantly, the carbonates indicate there was little or no compressional deformation within or adjacent to a large portion in the overriding plate in the collision zone, for up to 15 Myr after initial collision, assuming that initial collision was at 35 Ma. The carbonate deposits occupy a region 200-300 km from the suture zone and have a trench-perpendicular width of approximately 200-500 km (Fig. 6) (Morley et al., 2009). The geometry and timing of the carbonate basin in Iran and Turkey (Fig. 6) fits well with the identified collisional mantle dynamic basin from the modelling results (Fig. 4a, region 1). The CMDB width in our model is around 300 km (Fig. 4), which fits well with the area over which carbonate sediments are found across Iran (Fig. 6).

Solid
The CMDB is present in our modelled results for a period of 10 Myrs between collision and slab break-off. Such a pattern of topographic change fits well with the change in elevation needed to deposit the carbonate sediments that are currently found across Iran and southern Turkey. We therefore propose that slab steepening after collision (Fig. 3c) would have caused the subsidence necessary to deposit these carbonates in a CMBD.
This pattern of topography change calculated from our models is similar to that produced in (Duretz et al., 2011), but the absolute topography range, we produced using the elastic filter, is closer to actual examples of subsidence or elevation observed in collision zones. The models also reproduce well the scale of the dynamic topography modelled and observed at subduction zones: Husson (2006) predicts a dynamic back arc basin of 1.5-2 km during on-going subduction similar to the models in this study before collision when subduction is at a normal rate of around 10 cm yr −1 .
Another geographic area that offers possibility for observation of topography change produced during continental collision is the Italian Apennines. Here, slab detachment is proposed to have started around 30 Ma and a tear propagated north to south along Italy (Wortel and Spakman, 2000). van der Meulen et al. (1999) observed a period of basin formation followed by uplift using the sedimentary record. Migrating depocentres were interpreted as evidence of a slab tear propagating north to south. Using modern estimates of the position of the suture in Italy, these depocentres seem to be mainly located on the overriding plate with the maximum observed depression around 100 km from the suture (Ascione et al., 2012). Tentatively, these observed depocentres could have been created by the same slab steepening mechanism discussed here, and be analogous to our modelled CMDB. This comparison is complicated by large roll back rates (Carminati et al., 1998;Jolivet et al., 2009) for subduction adjacent to Italy. This roll back would put the overriding plate into extension thinning it and forming basins pre collision. However, dating of sediments from the basins (van der Meulen et al., 1999) puts the time of deposition after initial collision, around the time of slab break-off, when the overriding plate would be expected to be in compression. This therefore does still fit with our model but could also be explained by delamination proposed for this region (Channell and Mareschai, 1989).
In our models there are a number of assumptions that will affect the modelled topography. One of the main assumptions is that the overriding plate does not shorten or thicken during subduction or collision. Substantial shortening and thickening obviously occur at collision zones, which must influence the overall topography. Our model does not incorporate these features, so our resultant topography does not include these effects. Therefore, although our study might not fully model the expected topography in some regions, it does allow us to understand the contribution of dynamic topography through time. This lack of crustal shortening and consequent crustal thickening is particularly relevant to the subduction of stretched and thinned continental material during the early stages of collision (Ballato et al., 2011).

Conclusions
Our modelling work emphasises that changes in surface elevation form a useful tool to study the process of continental collision. Dynamic topography is expected in the form of a back-arc basin during ongoing subduction, whose depth is correlated to subduction speed. As collision starts, this basin deepens due to the steepening of the slab. Surface uplift is expected between this basin and the trench, caused by subduction of continental material. After slab detachment, the uplift migrates into the overriding plate, where the basin had previously been.
These modelling results fit well with the sedimentation record and topography on the overriding plate for the Arabia-Eurasia collision zone. Upper Oligocene -Lower Miocene carbonates deposited in between terrestrial clastics show that a basin around 300 km-wide contained a shallow sea during the Upper Oligocene -Early Miocene for a period of around 8 Myr (Reuter et al., 2007). Present-day high elevation of the region also fits with the expected evolution of uplift after slab break-off on the overriding plate, although there is also a contribution from internal shortening and thickening.

Sensitivity testing
This Appendix discusses the sensitivity of modelling results to the choice of subduction model parameters. Figure A1 illustrates a suite of model runs with different input parameters. For each model the maximum depth of the collisional mantle dynamic basin was picked in time and space. Here we can see how each input parameter affects the depth of this depression. The value for our preferred model presented in this study is highlighted in red. All models share the initial set-up as described in the method section of this study, apart from the setup condition under investigation. In Fig. A1a, we show how a reduction in surface yield strength reduces the depth of the basin. This is due to a smaller surface yield strength allowing material at the top boundary of the model to more easily deform. This has the effect of reducing the magnitude of all topography by lowering the normal stress placed on the boundary. Figure A1b shows how an increase in the weak zone viscosity minorly increases the depth of the basin. In our model the weak zone viscosity defines the coupling between the subducting and overriding plate. As the weak zone increases in viscosity this decreases subduction speeds but also holds the slab to the overriding plate with a greater force. This means that when steepening occurs it has more effect in models with a stronger weak zone. It should also be noted that for the higher weak zone viscosities the basin was formed much later in time and that subduction speeds were much slower.
In Fig. A1c, we show how the mantle wedge viscosity has little effect on the basin depth. This is due to the basin be-ing due to slab steepening which creates a flow over a much larger area than just the mantle wedge. This means that even with a stronger or weaker mantle wedges the stress is still transmitted to the overriding plate, thereby creating the topography. In Fig. A1d we show how the width of the mantle wedge has little effect on the basin depth. This is for the same reason as the mantle wedge viscosity in that the flow responsible for the basin is much larger that the mantle wedge. The lack of sensitivity to the mantle wedge width also shows that the topography and in particular the CMDB is not due to the placement of the mantle wedge. Which implies that our model setup is not controlling our topography signal. Figure A1e shows how the depth of the weak mantle wedge affects the basin depth. Here, we can see that, as the wedge is made shallower from our preferred model, the basin becomes shallower. If made deeper than our preferred model, little effect is observed. This is due to the mantle wedge being placed below the lithosphere in our preferred model. If the wedge is placed shallower than this, it allows the slab to partially steepen before collision, so a smaller basin depth is created after collision. Placing the wedge deeper has little effect, since the material in this area is already weak due to being at mantle temperature. Figure A1f shows the effect of the chosen effective elastic thickness on the collisional mantle dynamic basin depth. This shows that with larger effective elastic thickness, the depth of the basin is reduced. The choice of effective elastic thickness and its effects are further discussed in Appendix 2.
Overall, this Appendix demonstrates that the topography features are robust and not a product of the initial conditions. Although there are a number of initial setups that do affect the topography, the initial setup chosen in our preferred model offers a geologically reasonable subduction and collision setup. Our preferred model setup is not specifically tuned to the Arabia-Eurasia collision though does offer a good model for comparison to this collision zone.

Elastic Thickness
The effective elastic thickness T e of the lithosphere around the Arabia-Eurasia continental collision zone is relatively poorly defined given the uncertainty in plate thickness (Mckenzie and Fairhead, 1997;Maggi et al., 2000;Watts and Burov 2003) and the uncertainty in the strength of the newly formed suture. The effect of increasing the effective elastic thickness can be seen in Fig. A2 were we show the topography time map from our preferred model with 20, 30, 40 and 50 km effective elastic thicknesses. This figure shows that greater elastic thickness gives topographic features a greater lateral extent as well as reducing their amplitude. Despite quantitative differences changing the elastic thickness produces, the discussed topographic features of this work are present over the full range of T e values explored here. The exact T e value representative of the collisional tectonics settings discussed here is not the focus of this study, though our results show the need to include elastic thickness to produce realistic topography estimate.