Discrete Element Modeling of Dike-induced Deformation

Introduction: The Tharsis region of Mars is characterized by large volcanic and tectonic centers with distinct sets of graben systems. Many of the radially oriented grabens have been inferred to form in response to intrusion of magmatic dikes. This interpretation is based primarily upon early physical and numerical (boundary element) models that were developed originally to understand surface deformation associated with dike emplacement on Earth. In this study, we constructed and analyzed two-dimensional discrete element models to test the hypothesis of shallow dike emplacement and widening as a primary mechanism for the production of grabens on Mars. In particular, our models are designed to explore the extent to which a widening subsurface dike under varying emplacement conditions will induce near-surface and surface deformation. Methodology: The use of discrete element models allows for permanent deformation and material heterogeneity to be simulated. In contrast, boundary element models are limited by an assumption of homogeneous elastic behavior. We performed discrete element modeling of dike-related deformation using the code PFC2D [1]. PFC2D is based on the discrete element method that was developed initially for analyzing the mechanical interaction of granular materials [2]. With the addition of particle bonding, the method has been extended to a variety of problems in solid mechanics, including rock fracturing and seismicity [3,4,5,6,7]. The discrete element formulation requires selection of a set of micromechanical properties (e.g., friction, normal and shear bond stiffness and strength) that describe the interaction of the elastic particles with each other and the model boundaries (see [4], for details on the theoretical framework of PFC). The result is a set of micromechanical properties that provides the correct bulk (macroscale) material behavior, which includes both the elastic (recoverable) and inelastic (nonrecoverable) components characteristic of natural deformation. For all models, stratigraphic material was modeled by generating particles in an initial model area 10 km wide. The particles have a diameter range of 12.8 – 42.9 m scaled to the kilometer scale model boundaries (smaller particle size ranges generate billions of particles which overwhelm the computations). A half model space was used to conserve computation time. The particles are assigned uniform stiffness and friction values throughout the model area. Blue and red colored layering was added for visualization purposes only and does not affect material behavior. Once the model reached equilibrium (no large out-of-balance forces), the simulated dike was widened. The dike was assumed to have reached a neutral buoyancy level (i.e., stopped vertically ascending) 1.25 km below the model surface and widened in place. This depth to dike tip resulted in the most significant surface deformation while still remaining within the range of depths calculated by [8] for non-eruptive dikes. The dike is widened to 1 km in half-width (2 km overall width); although this exceeds the maximum widths calculated for Martian dikes [8], we found that a large dike width was necessary to produce significant surface deformation characteristics. Internal pressures within the dike are ignored and the dike walls are modeled by forcibly widening the dike into the surrounding model space. During model evolution, however, the stress state was calculated in the area adjacent to the dike plane with a horizontal stress magnitude of approximately 127 MPa, consistent with values reported for internal dike pressures [8]. Results: Previous modeling efforts have examined the role of a widening subsurface dike in the absence of regional extension and under various configurations of mechanical layered stratigraphy [9]. For comparison, the base model is presented here as model A. In model A (Figure 1), no bonds are applied to the particles and their behavior is consistent with unconsolidated rock or very weak regolith material. Dike growth produces compressional forces in the material adjacent to the dike plane and a corresponding topographic uplift at the surface, similar to results predicted by previous boundary element models [10,11,12]. Contractional folds are seen in the subsurface, producing surface deformation away from the dike for approximately six times the width of the dike (6.2 km). The topographic peak occurs at approximately twice the width of the underlying dike (1.9 km), and the magnitude of uplift for the unconsolidated regolith model is 0.6 km, corresponding to the topographic relief of the trough. It should be noted that for the unconsolidated model A, the topographic low of the trough immediately above the dike does not drop below the regional topographic surface.