Interpreting the interseismic deformation of the Altotiberina Fault ( central Italy ) through 2 D modelling

The Altotiberina low-angle normal fault in central Italy has been a focus of many recent studies. Although the existence of this fault has long been known, its seismicity and relationship to other faults are still debated. We present a 2D elastoplastic finite-element model that reproduces the interseismic deformation of the Altotiberina Fault. The model predictions are compared to observed geodetic velocities, stress orientations and geological data. The influence of the Altotiberina Fault on interseismic evolution is tested by building several models with different boundary conditions. The best model is 180 km long, 40 km deep and contains two layers with different rheological parameters, two ramps, two faults and four freely slipping segments. The main factors contributing to the largescale interseismic deformation include basal traction, rheology and the Altotiberina Fault itself, whereas the local, small-scale variations are due to two secondary high-angle faults.


Introduction
The interseismic phase is the period in the seismic cycle following post-seismic relaxation and preceding the next large earthquake.During the interseismic phase, it is useful to identify locked faults because they have the potential to be seismogenic [e.g., Chamot-Rooke and Le Pichona 1999, Simoes et al. 2004].The locked portion of a fault can accumulate stress, which is later released when the fault slips during an earthquake [e.g., Subarya et al. 2006, Konca et al. 2008, Moreno et al. 2010].
Several studies have described the interseismic phase, mainly by using geodetic or geophysical data [e.g., Mazzotti et al. 2000, Simoes et al. 2004, Brooks et al. 2011], which have different time scales.Geodetic data record deformation at time scales that are much shorter than the intervals between large earthquakes [e.g., Wei et al. 2010].The stress field derived using borehole data and the strain field derived using earthquakes describe present-day tectonic activity [e.g., Zoback 1992].Geological observables such as vertical velocities from uplift and subsidence rates [e.g., Cook andRoyden 2008, Huang et al. 2010] record several seismic cycles and help constrain tectonics over longer periods of time.
Interseismic deformation is often studied through analytical and numerical methods.Faults are embedded into an elastic medium while allowing the deeper portion of the fault to slip aseismically [Savage 1983, Ruegg et al. 2009].Alternatively, the fault can be embedded into a crustal structure that is modelled as elastic upper crust overlying viscous lower crust or elastic crust overlying the upper mantle [e.g., Vergne et al. 2001, Correa-Mora et al. 2008, Huang et al. 2010, Doglioni et al. 2011].Typically, these models are used for convergent or transform plate boundaries [e.g., Marshall et al. 2009, Vigny et al. 2009], where the alongstrike length and the depth extent of the fault are both very large.These models are able to reproduce the interseismic deformation for large thrust and strike-slip faults, but they are not able to reproduce deformation induced by low-angle normal faults (LANFs).
Indeed, LANFs activity is not predicted by traditional Anderson-Byerlee fault mechanics.Many authors have explained the existence of the LANFs by using mechanical theories with deformation localised at the interface of elastic and viscous layers [e.g., Melosh 1990] or by using an additional horizontal force [Yin 1989, Westaway 1999].
Among LANFs, the Altotiberina Fault (ATF) is an example of a basal detachment in the extensional area of central Italy as shown in the CROP03 deep-crust reflection profile [Barchi et al. 1998].The existence of the ATF is well documented in geological [e.g., Brozzetti 1995, Brozzetti et al. 2000] and geophysical data [e.g., Mirabella et al. 2011, and references therein].The ATF borders the west side of the Tiber Basin with a NNW-SSE strike for approximately 70 km (Figure 1) and dips 10°-30° eastward to a depth of approximately 12 km [Barchi et al. 1998], which suggests a downdip extent of 35 km or more.The displacement along the ATF ranges from 5 to 8 km [Mirabella et al. 2011] and has developed since 3 Ma (Late Pliocene) [Brozzetti 1995].The Tiber Basin is a 6-km-wide tectonic depression filled by Plio-Pleistocene deposits and Holocene lacustrine and fluvial deposits [Boncio et al. 2000].At the surface, the ATF produces a number of low-angle (30°) synthetic splays [Boncio andLavecchia 2000a, Brozzetti et al. 2009].A series of antithetic high-angle (30°-60°) normal faults (HANFs), such as the Colfiorito Fault (CFF), the Gubbio Fault (GF), and the Umbria Valley Faults (VUFs), are located in the hanging-wall block of the ATF (Figure 2).These faults are related to shear stresses at the base of the brittle upper crust beneath the Apennines, which is caused by eastward mantle flow [Doglioni et al. 1998].As the ATF belongs to the large Etrurian Fault System, which extends approximately 350 km, LANFs may be widespread in central Italy [Boncio et al. 2000].
Earthquakes of moderate magnitude, such as the Norcia 1979 (M=5.9),Gubbio 1984 (M=5.2), and Colfiorito 1997 (M=5.9)earthquakes, have occurred within the hanging-wall of the ATF [Barba and Basili 2000, Boncio and Lavecchia 2000b, Boncio et al. 2000].The present-day activity of the ATF is characterised by small-magnitude earthquakes [Boncio et al. 1998, Carminati et al. 2001, Chiaraluce et al. 2007], mostly in its southern section [Piccinini et al. 2003].Based on the lack of microseismicity in the 20-km-long northern section of the ATF, Piccinini et al. [2003] suggest the existence of a locked portion along the ATF, which may be considered a seismic gap and thus may be hazardous.Other authors have assumed an Andersonian state of stress and model the ATF as creeping [e.g., Collettini and Barchi 2002] and nearly aseismic.The choice between these different interpretations is debated and is critical as it affects the hazard estimates of the area.
The ATF has been modelled for various purposes.Using a two-dimensional viscoelastic model of Apennine subduction along the CROP03 seismic profile, Carminati et al. [2001] modelled the ATF as a structural and rheological discontinuity embedded within viscoelastic and elastoplastic stratifications.They computed the stress directions and intensities surrounding the ATF and the dynamic topography along the modelled section.By comparing the model predictions with the dip of the T axes of earthquakes and the earthquake distribution at depth, Carminati et al. [2001] proposed that the ATF re-orients the v3 stresses at depth, increasing their dip and becoming nearly par- allel to the ATF.These results provide a locking mechanism for the ATF and suggest that the local stress field is non-Andersonian.Pauselli and Federico [2003] developed a two-dimensional elastic stratified finite-element model in which the faults are modelled as areas of decreased Young's modulus.This model allows the stress directions to be computed along the CROP03 seismic profile for different stratifications and assuming different hypotheses about the Gubbio Fault.Pauselli and Federico [2003] found that the inclusion of the Gubbio Fault along with the ATF in the model is essential to reproducing the spatial pattern of seismicity at depth and the moderate-magnitude earthquakes.Hreinsdòttir and Bennett [2009] modelled the ATF as a variable-dip fault in an elastic half space.By using a constrained random search algorithm, they determined the locking depth of the ATF and the slip rate that best reproduce the vertical and fault-perpendicular components of the horizontal velocities of continuous GPS sites in the region.
In this paper, we assumed that shear stresses at the base of the plastic lower crust are the cause of fault motion, following Westaway [1999].The aim is to demonstrate how the shear basal tractions contribute to large-scale interseismic deformation.We applied this model to the ATF and associated HANFs, and we implemented the main characteristics highlighted in previous studies.We developed an elastoplastic finiteelement model where the ATF is embedded as a discontinuity as proposed by Carminati et al. [2001], and introduced antithetic normal faults, as implied by Pauselli and Federico [2003].We compared the model predictions with the interseismic fault-perpendicular component of the horizontal velocities of continuous GPS sites in the region, similarly to Hreinsdòttir and Bennett [2009], and with stress directions and geological data.We adopted two layers with elastoplastic rheology and a realistic subsurface geometry.Given the time scale of interest (a few years), we neglected the poroelastic effects and focused on crustal heterogeneities.This allowed us to determine the locked (or unlocked) state of the faults.
In this paper, we present different models, and we discuss the geological implications of the best fit.

Data and model
We aimed to reproduce the interseismic data through finite-element numerical modelling.Starting from an initial model, we varied the rheological parameters, fault parameters and boundary conditions, and we chose the combination that best reproduced the geodetic velocities, stress orientations and surface geology.In all models, an unlocked fault is represented by a freely slipping segment.

Data
A dataset of 182 continuous GPS velocities was analysed by Caporali et al. [2011].Forty-seven velocity data points are within the study area and are expressed with respect to the Eurasia-fixed reference system defined by Devoti et al. [2008].The horizontal GPS velocity data were projected along model section A-A', shown in Figure 1, and the section-parallel component of motion (hereinafter referred to as V h ) was considered.By selecting subsets of GPS data located at different distances from the section (20, 40, 50, 70, and 150 km), we verified that the GPS velocities are largely homogeneous.Thus, we selected all data within 150 km of the section (Figure 1).
The orientations of the maximum horizontal stress (S Hmax ), available from the World Stress Map 2008 (WSM) [Heidbach et al. 2010], were projected along model section A-A', and data were plotted at the same ranges used for the GPS velocities (20, 40, 50, 70, and 150 km).We used data within 70 km of the section, which includes 150 S Hmax values (Figure 1).
We filtered both the stress azimuth and the sectionparallel components of the horizontal GPS velocities using spatial Gaussian smoothing, i.e., by convolution with a Gaussian function of a certain width.We used a 50-km width for both the GPS and WSM data (Figure 3).The V h values increase to the NE, and the S Hmax data indicate extension corresponding to the Tiber Basin.
The interseismic model should predict the location of the subsiding area of the Tiber Basin (Figure 1).Thus, we used the data of Mirabella et al. [2011] to constrain the direction of vertical displacement in the model.

Geological section and model geometry
From the geological sections shown in Figure 2, we built the model based on a simplified configuration of the principal faults and stratigraphy.We divided the model into a lower layer (LL), which includes most of the acoustic basement, and an upper layer (UL) representing the Triassic-Miocene stratigraphic sequence of the Umbria-Marche.The interface between these two layers represents the elastoplastic transition (EPT).The EPT is shallower than the brittle-ductile transition below the ATF and corresponds approximately to the brittleductile transition elsewhere.This geometry allows the ATF to be within the brittle (elastic-frictional) layer.Its depth increases to 15 km to the east, in agreement with Brozzetti et al. [2009] and Carafa and Barba [2011].
The ATF is modelled by a low-angle (20° east-dipping) fault along the EPT that is either locked or unlocked.The locked case corresponds to an interface with elastic-frictional properties (normal friction; e.g., µ = 0.6-0.8).The unlocked case corresponds to an interface with free-slip properties to represent low-friction aseismic slip.Following Doglioni et al. [1998], we included the regional monocline (RM) between the Apennines and the Adriatic Sea (Figure 2c).
We modelled two of the principal HANFs (represented as faults A and B in Figure 4) and a series of tentative faults that vary in dip and are considered freely slipping segments.The HANFs cut through the UL and do not offset the LL.

Numerical model
We developed a 2D finite-element model using MSC.Marc 2010 [MSC.Software Corporation 2010].The displacement was computed assuming elastoplasticity and plane strain.The model was run for the 11 years, corresponding to GPS time span.Our model is 180 km long and 40 km deep.The two layers have elastoplastic rheology, and the yield stress in the UL is greater than in the LL (Table 1).This choice lets the UL reach the plastic limit at localised points along the slipping segments and lets the LL be close to the plastic limit to exhibit a small amount of diffuse plastic strain.
The grid includes four-node quadrangular elements.The number of elements ranges from 7000 to 8000 (the area of each element is ~1 km 2 ), and there are 7500 to 8500 nodes, depending on the faults considered.The faults are modelled through duplicate nodes on both sides of the interface [Melosh and Raefsky 1981].In the unlocked case, the fault nodes allow relative movement with zero friction between the footwall and hanging-wall.Reaction forces prevent the interpenetration of the footwall and hanging-wall (see e.g., Megna et al. [2005]) except for a negligible amount that is required for numerical reasons.
The boundary conditions for the right and left sides are containment edges, i.e., the material cannot leave the model, but it is free to move inside it.The nodes along the bottom of the model are locked vertically and freely slip horizontally.A shear basal traction is applied to the bottom of the model with values from 1 to 10 MPa.These tractions generate flow in the plastic lower continental crust, which is transmitted to the upper crust.The model is subject to gravity (9.8 m/s 2 ).The Poisson's ratio (h) is 0.3, and the density (t) is 2300 kg/m 3 for both layers.The yield stress (Y) ranges from 10 to 200 MPa, and the Young's modulus (E) ranges from 1 to 50 GPa.We developed a series of models by varying the rheological parameters, geometries, and faults to test the effects of each parameter on the surface velocities and vertical displacements.We varied Y from 30 to 200 MPa in the UL and 10 to 40 MPa in the LL.The values of E ranged from 1 to 50 GPa for the UL and from 1 to 10 GPa for the LL.We then varied the dips of the HANFs from 30° to 60° (Figure 4), imposing different widths of free slip along each fault.We also tried different values of the model length (100 to 200 km), the model thickness (20 to 40 km) and the distance of the ATF ramp from the left edge of the model (10 to 60 km, see Table 1).For simplicity, we show geometries (Figure 5) and results (Figure 6) of four end-member models.In Model A, we tested a simple crustal structure incorporating two plane-parallel layers (the UL and LL), the containment edge boundary condition, and no faults.In Model B, we introduced a low-angle ramp along the EPT at the left edge of the model.We built Models C and D to verify the effect of the opposing fault dip (NE vs. SW).These models include a free-slip portion along the EPT and RM.We calculated the root mean square (RMS) for each model.The RMS is an indicator of the difference between the filtered GPS data and the model predictions.The lowest RMS indicates the best-fitting model.
We obtained the horizontal velocities at the surface from the horizontal displacements versus the time factor at the top edge of the models and compared these curves with filtered V h and S Hmax data.The time factor corresponds to the duration of the GPS data (11 years).The vertical displacement was calculated at the surface of the model and was compared with the position of the area of subsidence.Thus, the best-fitting model was chosen, and the model curve that is closest to the data represents the best model.To better evaluate the significance of the model results, we computed the vertical and horizontal displacement fields in the model section due only to the boundary conditions, in the pre-stressed model, and for the total stress and plastic strain, which includes the gravity load.

Results
By varying the model parameters (Table 1), we ran approximately 200 models and chose the one that best fit the data.We focused on the locking state of the faults in the best model, together with the rheology in the interseismic period.Of all the tests, we show only the most relevant outcomes and describe the contributions to the surface velocities, the vertical displacements of the EPT, the Altotiberina Fault (ATF), detachments, unlocked faults, and the rheological parameters.
Model A shows that the horizontal velocities increase in the SW half of the section (Figure 6); this is due to the presence of the EPT (the interface of the two layers), together with the shear basal traction (4.3 MPa).We observe a wide area of slight downward displacement at the left edge of the model (from 0 to 60 km; Figure 6b).In Model B, the presence of the ATF controls the area of subsidence and the horizontal velocities (Figure 6a,b).We observe a step in the horizontal velocity and subsidence, in the left portion of the section.In Model C, the dip of the faults was reversed (from SW to NE), causing a reversal in the patterns of local subsidence and uplift (Figure 6b) that does not produce observable variations in the horizontal velocities (Figure 6a).The lowest RMS (0.1 mm/yr) was obtained for Model D, which is deemed the best model.
Model D shows that the free slip on the two highangle faults (faults A and B) produces steps in the horizontal velocities near each fault (at 20 km and 60 km from the origin) (Figure 6a).In the other tests, the rheological parameters control the surface velocities and vertical displacements.By changing the yield stresses for the upper and lower layers, we obtained an increase or a decrease in velocity (Figure 6c).In particular, the velocities increase if the LL (i.e., Y = 10 MPa) is more plastic than the UL (i.e., Y = 30 MPa).In Model D, free slip is also predicted along the EPT at 90 km from the origin, where the EPT is horizontal.This result represents the relative displacement between the upper and lower layers along the brittle-to-ductile transition, where it coincides with the EPT.The free slip found on the RM suggests a mechanism to facilitate thrust propagation with a low seismic coupling factor.
The predictions of the best-fitting model closely agree with the data.In the first 80 km, both the modelled horizontal velocities and the filtered V h increase to 2.5 mm/yr, indicating that the model is under extension (Figure 6a).The velocities decrease at greater horizontal distances, suggesting an area of compression, also in agreement with the WSM stress data (Figure 7a,b).The modelled vertical displacements show subsidence from 0 to 60 km (Figures 7 and 9a), with the largest amount of subsidence in the Basin.The high-angle faults (unlocked faults A and B) reproduce the irregularities in the filtered V h located at 20 km and 60 km from the origin (Figure 7b).At these two points, the filtered V h is steeper, increasing by 0.5 mm/yr within 10 km.Furthermore, the vertical displacements of faults A and B increase the subsidence of each fault by 25 m (Figure 7c).
In summary, the best model is characterised by a yield stress of 50 MPa for the UL and 20 MPa for the LL, a Young's modulus of 9 GPa for the UL and 10 GPa for the LL, a fully locked ATF, and a 16-km-wide freely slipping fault for the RM.Along the EPT, a 21-km-wide freely slipping segment is located 45 km from the origin.The best model includes two high-angle faults in the UL: an unlocked 9-km-wide freely slipping fault that dips 50° SE (fault A) and an unlocked 20-km-wide fault that dips 40° SE (fault B) (Figure 8).The position of fault A corresponds to the Umbria Valley faults (VUFs) near the surface trace of the ATF (shown in Figure 2b), and fault B corresponds to the Colfiorito Fault (CFF).The geometric and elastic parameters for the best model are listed in Table 1.
The vertical and horizontal displacements, the stress, and the plastic strain (Figure 9) suggest that the faults and the basal traction influence the model results.The vertical displacements (Figure 9a) highlight local subsidence in agreement with the position of the Tiber and Colfiorito Basins (0 to 60 km from the left edge of  the model).The horizontal displacements (Figure 9b) increase at depth, corresponding to the lower tip of the ATF.The stress (Figure 9c) primarily accumulates in the UL and the plastic strain (Figure 9d) increases with depth, with perturbations induced by the slip along the unlocked segments.In the UL and the upper part of the LL, the model mainly behaves elastically, with a neglectable component of plastic strain.

Discussion and conclusions
The shear basal traction appears to be essential factor to modelling the large-scale interseismic pattern in the presence of the low angle fault.This traction allows the reproduction of extension and compression rates for central Italy as also shown by Barba et al. [2008] in a larger-scale model.The lack of a clear discontinuity in the fault-perpendicular horizontal velocities indicates that the ATF is at least partly locked.On the other hand, the high-angle faults induce local variations in the vertical velocities.
The locked state of the ATF and the prior knowledge that it is located within a brittle layer indicate that the ATF is possible seismogenic.In the sense of scaling laws, the locked width corresponds to earthquakes up to magnitude 7.
The main limit of our method is that it is two-dimensional.The rotations of the horizontal velocity vectors and the stress axes that occur far from the modelled sections may be due to local features in the sub-crustal structure or to the presence of minor or unfavourably oriented faults.The two-dimensional geometry also does not allow the discrimination between parallel faults of similar size; in our case, the GF is parallel to the CFF, and in practice, the resulting unlocking width of the CFF may also depend on the GF.Additionally, modelling a simplified stratigraphy induces long-wavelength patterns, especially in the predicted vertical displacements.In our case, the predicted subsidence is wider than the width of the Tiber Basin.Consequently, our model may underestimate the locked width of the VUFs.
In regards to the relationship between the locked width and the potential for earthquakes, we make note of two points: the first is related to historical earthquakes, and the second is related to the seismogenic character of a locked fault.Chiaraluce et al. [2007] noted that if the ATF is seismogenic, historical seismic catalogues should report large-magnitude earthquakes that could be associated with such a large fault although this point is debated.Historical earthquakes have large uncertainties in location and magnitude.Moreover, large earthquakes are not periodic, so comparing the seismic catalogue to the recurrence time of large earthquakes does not provide sufficient evidence to conclude that the ATF is aseismic.For the future behaviour of a locked fault, it is impossible to know if a locked fault will rupture as one large earthquake, as several smaller ones, or through episodes of silent slip, as observed in the central Apennines [Amoruso et al. 2002, Scarpa et al. 2008].Therefore, we suggest considering the steeper segments of the ATF in estimates of seismic hazard, as suggested by other authors as well [e.g., Piccinini et al. 2003, Brozzetti et al. 2009, DISS Working Group 2010].
The slip along the unlocked portions of the VUFs, GF, and CFF as well as the small plastic strain below the ATF may correspond to the frequent microseismicity observed in the study area [Boncio et al. 1998, Carminati et al. 2001, Chiaraluce et al. 2009].The release mechanism of a small amount of accumulated stress into a small amount of plastic strain is an alternative to the one proposed by Collettini and Barchi [2002], essentially differing because of the non-Andersonian state of stress.
The GPS data used to calibrate the model span a time interval of approximately 11 years, from 1998 to 2008.This period follows the Colfiorito earthquake, which occurred in 1997.The unlocked state of the CFF thus corresponds to the behaviour after the 1997 earthquake and is not representative of the times prior to 1997, when the fault may have been locked.
Including tentative faults in the model geometry allows for evaluation of the sensitivity of the model to the fault geometry.The numerical model helps address some characteristics of the faults and detect unmapped buried faults.This approach is indispensable when modelling areas where little is known about the fault geometries.The results of the model are reasonably good in central Italy, which is one of the best-studied areas of normal faulting in the central Mediterranean.Calibrating the method along the ATF allows generalising it for applications to other areas where the faults are poorly known.For the ATF, the model presented here constitutes a starting point to determine the slip rates of the modelled faults.
Our results provide insight to the locked state of faults.The model that best fits the observed data indicates that the ATF is entirely locked due to the lack of a discontinuity in the horizontal and vertical interseismic velocities.The presence of a discontinuity in the horizontal and vertical velocities suggests that the CFF and VUFs behave as unlocked.

Figure 3 .
Figure 3. (a) Section-parallel GPS velocity data (V h ) projected along model section A-A'.Black and grey circles represent the data located within 40 km and 150 km of the section, respectively.The black and grey lines represent the result of Gaussian smoothing (50 km width) applied to the black and grey circles.(b) The azimuths of the WSM data (S Hmax ) projected along model section A-A'.Black circles and grey circles show the data within 20 km and 70 km of the section, respectively.The black and grey lines represent filtered versions of the datasets.

Figure 5 .
Figure 5. Four tests: Models A, B, C, and D. Model A contains two plane-parallel layers, Model B includes the ATF, Model C includes fault A dipping to the NE and Model D includes fault A dipping to the SW (see Figure 4 for acronym definitions and the text for descriptions of the models).

Figure 6 .
Figure 6.(a) The horizontal velocities and (b) vertical displacement for Models A, B, C, and D. (c) The yield stresses for the best model (Model D).The yield stresses for the tests are as follows: for Test Y1, UL: Y=30 MPa and LL: Y=10 MPa; for Test Y2, UL: Y=50 MPa and LL: Y=20 MPa; and for Test Y3, UL: Y=100 MPa and LL: Y=25 MPa.The red line represents the filtered GPS data from within 150 km from the model section (see grey line in Figure3).The RMS values represent the mismatch between the model horizontal velocities at the surface for Models A-D and the red line.

Figure 7 .
Figure 7.The results for the best model.(a) The azimuths of the WSM (S Hmax ) filtered data are in grey for the extensional area and black for the compressional area.The dashed line represents the results for the best model.(b) The horizontal velocities (V h ) for the filtered data and the best model.(c) The vertical displacements and corresponding areas of subsidence in the Tiber Basin, shown within the square.

Figure 8 .
Figure 8.The geometry and grid for the model with the best-fitting boundary conditions, with the free slip segment (grey rectangle).The circles represent nodes that are locked in the vertical direction, horizontal and dipping bars represent the containment bars, and arrows indicate the direction of the shear basal traction (4.3 MPa).UL: Upper Layer; LL: Lower Layer; ATF: Altotiberina Fault; RM: Regional Monocline; GF: Gubbio Fault; CFF: Colfiorito Fault; and VUFs: Umbria Valley Faults.

Figure 9 .
Figure 9.The results for the best model (Model D).(a) Vertical displacement, (b) horizontal displacement, (c) stress, (d) plastic strain (see text for descriptions).Dashed lines represent the unlocked portions of modelled faults.

Table 1 .
The parameters used for the lower layer (LL) and the upper layer (UL) in the tests and the best model.
Z: the approximate depth range of the layer; E: Young's modulus; h: Poisson's ratio; t: density; Y: yield stress; L: model length; l: left side; r: right side; and SH: shear basal traction.