Modeling Multi‐Material Structural Patterns in Tectonic Flow With a Discontinuous Galerkin Level Set Method

We formulate a numerical framework, in both 2d and 3d, to model the structural patterns emerging from creeping viscous flow typically encountered in long‐term ductile lithospheric deformation by coupling the discontinuous Galerkin level set method with a finite element Stokes‐like flow solver. The level set formulation has the advantage of retaining information on the interface geometry, decreased memory requirement and improved computational efficiency from the two‐way particle‐mesh information transfer compared to particle‐in‐cell methods. Furthermore, our formulation fully exploits the advantages of the finite element method (e.g., the flexibility of mesh geometry and the ease of handling anisotropic materials) by using a unified finite element framework. The novelty of our formulation is the capability to offer a fully dynamic approach for modeling structural patterns resulting from a tectonic flow that is non‐steady and inhomogeneous. The material distribution and the finite deformation patterns predicted from the numerical model can be directly compared with geological map patterns (e.g., lithological distribution at specified depths and on cross‐sections) and field structural analyses (e.g., foliation, lineation and strain patterns), thus offering the possibility of ground‐truthing the modeling results by field evidence. As examples for potential applications of our method, we apply our method to the modeling of a competent inclusion in simple shear flow, as well as a Rayleigh‐Taylor type density overturn. Our models demonstrate good agreement with previous 2d benchmark results and produce 3d lithological and deformation patterns comparable to field observations.

The objectives of our tectonic modeling method are twofold.First, our method implements an efficient and flexible method for capturing material interfaces on an Eulerian mesh, in both 2d and 3d, without the necessity to involve additional Lagrangian mesh or particles.Second, we build the connection between abstract numerical models of crustal dynamics and concrete observations of geological mapping by formulating methodology to extract field-testable lithological distributions and deformation patterns from modeling results.
To achieve the first objective, we use the level set formulation, which is a method that is widely applied to the simulation of multimaterial flows in fluid dynamics (e.g., Kronbichler et al., 2018;Sussman et al., 1994).The level set formulation implicitly "captures" the interface between material domains by attaching a scalar field to the Eulerian mesh (Bourgouin et al., 2007;Hillebrand et al., 2014;Samuel & Evonuk, 2010).The solution of the level set formulation is achieved by the discontinuous Galerkin (DG) discretization, which is a variant of the finite element method optimized for the modeling of transport-type problems.Adaptive mesh refinement localized around the material interface is implemented to improve the accuracy of interface resolution and minimize the mass loss deficiency inherent to the level set method.To build the connection between the material distribution and the tectonic flow, we couple the DG level set solver with a continuous finite element Stokes flow solver.Additionally, we employ an efficient linear system solver technique to obtain fast solution on an adaptively refined mesh with viscosity contrasts.The strength of this combination is the unification of the level set solver and the flow solver under the finite element framework, which then enables us to fully exploit the advantages of the finite element method in terms of flexibility of mesh geometry, simplicity of implementing anisotropic rheology to model mesoscale structures (Jiang, 2016;Ran et al., 2019) as well as large-scale mantle processes (Fraters & Billen, 2021), and accurate treatment of material transport due to the highly scalable DG discretization (e.g., Hesthaven & Warburton, 2007).Our methodology is similar in spirit with Heister et al. (2017), with the difference that the latter computes the advection of the composition fields directly without retaining the geometrical information of the interface.
To achieve our second objective we combine the strengths of the level set description of the interface geometry, the kinematic analysis techniques in structural geology and the mesh-based geodynamic models.We seek a numerical solution of tectonic flow on a discretized mesh, which results in a simulated tectonic flow that is non-steady in time and inhomogeneous in space.Additionally, the distribution of the lithological units are recorded by the evolution of the level set field while the finite deformation field is continuously recorded by the massless tracer particles.A simulated geological map pattern can be constructed from the level set field at any arbitrary level of depth.Regional folding of the lithological units can be visualized by constructing cross-sections WÚ ET AL.

10.1029/2023JB026806
3 of 27 perpendicular to the regional structural grain.Analysis of the deformation field on the stationary mesh as well as on the tracer particles yields structural information (e.g., bedding, foliation and lineation orientations) that can be compared with structural features observable in the field.We note that this approach differs from the traditional PIC in that the markers are entirely passive and therefore not essential to the modeling computations.The modeler has the liberty to choose a suitable number of particles in areas of the model where deformation tracking is needed.
Our modeling method is particularly suited for testing numerical models with observations acquired from field mapping in the form of the distribution of lithological units, regional folding patterns of stratigraphy, the distribution of strain characteristics and intensity, as well as the foliation/lineation patterns.Our present formulation is not confined to any particular scale as long as the viscous rheological law is applicable and the inertial effect of the flow is negligible (the "creeping" flow regime), thereby making it a versatile tool to explore instantaneous and finite strain patterns resulting from both regional crustal-scale deformation (e.g., Archean lithospheric-scale density inversions, Section 3.3) and the evolution of meso-to microscale structural features (e.g., deformation around rheologically distinct inclusions or porphyroclasts, Section 3.1).
In the following sections of the article we begin by briefly introducing the discretization scheme for both the tectonic flow solver and the level set interface solver in Section 2 as well as their mutual coupling.Also discussed in this section is the methodology to extract instantaneous and finite deformation patterns (Section 2.3) from the numerical model.We then present numerical experiments in Section 3 to showcase the validity of our modeling methods and its capability to produce simulation data that is relevant to field structures analysis.Finally, we discuss the correlation between the modeling results with real-world geological observations and our plan for future developments in Section 4.

Numerical Methods
In our method, we couple two solution mechanisms to achieve multi-material simulation of the tectonic flow.The Stokes solver obtains the flow velocity and pressure while the level set solver captures the evolution of regions occupied by materials of different properties.The level set method has been widely applied to the modeling of two-phase flows (e.g., Sussman et al., 1998), free-surface flows (e.g., Groß et al., 2006) and crystal growth (e.g., Merriman et al., 1994;Zhao et al., 1996).It has also been applied to the modeling of tectonic flows under a finite volume (Hillebrand et al., 2014;Samuel & Evonuk, 2010) or a traditional finite element discretization (Bourgouin et al., 2007).
The level set solver is discretized with a DG discretization.To the best knowledge of the authors, our study represents the first attempt to introduce DG level set formulation to modeling viscous tectonic flow.As the discontinuous variant of the finite element method, the DG method relaxes the inter-elemental continuity requirement in the traditional continuous finite element method and reintroduces the numerical flux into the finite element formulation.With this property the DG combines the stability of finite volume discretization for the modeling of hyperbolic problems and the compactness for high-order approximation of the finite element method.Because the approximation basis is local to each mesh cell, the DG method lends itself to massive parallelization and adaptive mesh refinement (Hesthaven & Warburton, 2007).
In the following paragraphs of this section we describe the numerical strategies used to implement both solvers as well as the methodologies to extract information regarding finite deformation from the model results.The numerical procedures herein discussed are implemented in both 2d and 3d based on the finite element infrastructures provided by the deal.II library (Arndt et al., 2021).

The Stokes Flow Solver
The long-term tectonic flow of the Earth's crustal material can be approximated by a highly viscous creeping flow.Furthermore, if the volume change is negligible, the flow can be approximated as incompressible, which leads to the following set of governing equations: where is the rate-of-deformation tensor and g is the constant downward-pointing gravity vector.ρ(x) and η(x) are the spatially-variable material density and the effective viscosity respectively.At the current stage of development the effective viscosity η is a scalar due to the assumption of isotropic material.However, it is possible to model anisotropic material behavior by extending viscosity to a tensorial quantity.A mixed finite element method is used to discretize the equation on a Eulerian quadrilateral mesh with a quadratic approximation for velocity and a linear approximation for pressure (i.e., the  21 elements) (Donéa, 2003).Complemented by the appropriate BC, the velocity v and pressure p can be obtained on the mesh by solving the following linear system: where A is the discrete viscosity matrix and B is the discrete gradient operator. v and  p are arrays of expansion coefficients for velocity and pressure respectively.Equation 2 can either be solved monolithically with a fully coupled (FC) approach (Silvester & Wathen, 1994) or in a two-step manner by decoupling the pressure from velocity by Schur complement reduction (SCR) (Elman & Golub, 1994;L.-N.Moresi & Solomatov, 1995).In both cases, solving the linear system efficiently hinges upon finding a high-quality preconditioner for A. In this study we employ a level-based incomplete LU (ILU) preconditioner as well as an algebraic multigrid (AMG) preconditioner with one V-cycle sweep.An in-depth discussion of these solver techniques as well as performance comparisons are available in the supplementary material "Linear system solver details" (Data Set S1).In summary, we conclude that FC excels relative to SCR regardless of the choice of the preconditioner.In terms of the preconditioner for the FC approach, the ILU with a level-of-fill equal to one is the most cost-effective method for 2d problems but it is no longer as competent as the AMG method for 3d problems due to the large number of induced fills.

The Level Set Solver
While the Eulerian formulation of the Stokes flow solver provides the dynamics of the flow, the passive movement of the regions of different material properties cannot be resolved, thus necessitating the formulation of some material tracking mechanism.In this study, we utilize the level set method as a strategy to capture the material interface rather than representing the material domains with particles (e.g., Gerya & Yuen, 2007;L. Moresi et al., 2003) or tracking the interface with an independent mesh structure (e.g., Braun et al., 2008).
The level set method (Figure 1) captures the interface of a material domain D by associating it with a scalar level set function ϕ D (x, t).The magnitude of ϕ D has a signed distance property, that is at an arbitrary point x, ϕ D (x) is the distance from the point to the interface of D. Whether the point is "inside" or "outside" of the domain is reflected by the sign of ϕ D .The level set field evaluates to zero on the interface ∂D.In this study, we follow the convention that a negative level set value is prescribed to a point inside D and a positive value if otherwise.These properties can be encapsulated in the following equation: where d(x, ∂D) is the distance function from a point x to the interface of a domain D occupied by one single type of material.Due to its property as a signed distance field, the level set field has a unit gradient (‖∇ϕ‖ = 1).
Modeling the evolution of the material domain with the level set method involves the solution of an advection problem followed by a reinitialization problem.Since both problems are of hyperbolic nature, we implement the DG method to discretize the weak form of the equations representing both problems.The temporal integration is treated explicitly with a high-order explicit Runge-Kutta total variation diminishing method.

The Advection Problem
The advection problem describes the passive movement of the level set field ϕ under the velocity field v provided by the flow solver: Because the DG discretization does not enforce continuity across cells, the discretized level set value ϕ h is doubly-valued at the cell face.We use the following shorthand to describe the one-sided value of the interfacial ϕ h from the viewpoint of cell Ω e : where n e is the outward-pointing normal at the cell face of the element Ω e .From a geometric standpoint,   − ℎ, is the value of ϕ h at the cell face approximated by the polynomial basis internal to Ω e , while   + ℎ, is that approximated by the basis local to neighboring cell sharing the same cell face.
By premultiplying with an arbitrary test function ψ and integrating over a finite element cell Ω e , we arrive at the DG discretization of the advection problem: where the subscripts ϕ h and ψ h are the polynomial approximation of the level set function ϕ and test function ψ. f(⋅, ⋅) is the numerical flux located at the interface of the cell Ω e and is upwinded with respect to the transport direction of the velocity field v for stability: Figure 1.Concept of the level set function for the interface of an elliptical material domain on a 2d plane.x 1 is inside the material domain D, x 2 is outside the material domain and x 3 is on the material interface.d i is a shorthand of the distance to interface d(x i , ∂D).In this case ϕ(x 1 ) = −d 1 , ϕ(x 2 ) = d 2 and ϕ(x 3 ) = 0.

The Reinitialization Problem
As the level set field is advected with the velocity of the flow through time, the magnitude of the level set gradient is not guaranteed to maintain unity.As a result, the level set field loses the signed distance property.Because the material parameters in the cells close to the interface are computed as a function of the level set values (see Section 2.2.3 for details), this deviation affects the numerical integration of material parameters in these cells.If left untreated, this deviation can then cause large errors in the simulated results.To correct for this deviation, we reinitialize the level set field by evolving the following reinitialization equation through pseudo-time τ to steady state: where  0() = (, )| =0 is the not-yet-reinitialized value of ϕ and sgn is the signum function.The physical interpretation of this equation is to propagate the information of the unit divergence away from the interface.This is a Hamilton-Jacobi problem with the Hamilton-Jacobi operator H in the form of  (∇) = sgn(0)(‖∇‖ − 1) .A local DG scheme is used to solve the reinitialization equation directly (Karakus et al., 2016;Yan & Osher, 2011).Numerical instabilities that may develop in this formulation are smoothed by applying a weighted essential non-oscillatory limiter (Zhong & Shu, 2013).Detailed derivations of both the discretization and the limiting process are documented in Text S1 and S2 of the supplementary material "Level set details" (Data Set S2).

Flow Solver Coupling
Figure 2 illustrates the coupling between the Stokes flow solver and level set solver as a two-way process.From time t n to t n+1 , the discrete velocity field  v , obtained at each degree of freedom from the Stokes flow solver, is interpolated onto the quadrature points of each element cells via the finite element basis functions of velocity in the Stokes solver.The velocity values at the quadrature points are then employed in the level set update from ϕ n to ϕ n+1 through volume and surface integration in Equation 6.
To inform the Stokes flow solver of the updated material distribution, the level set values ϕ n+1 are first interpolated onto the quadrature points of the Stokes flow solver.The relationship between the updated level set field and the material parameter θ (e.g., density ρ, viscosity η, etc.) in a two-material system is established by: where θ 1 is the material parameter associated with the domain where ϕ < 0 and θ 0 is the "background" material.h ϵ (ϕ) is the smoothed Heaviside function: This method, which is typically used in two-phase flow problems (e.g., Tornberg & Engquist, 2000), creates a transition zone near the material interface of width 2ϵ so as to improve the accuracy of the numerical integration.The material parameters computed from Equation 9 will participate in the assembly of the elemental viscosity matrix in the Stokes flow solver.
In the case where more than two materials are present in the model, we define a background material which underlies all material domains.The overlying materials are each associated with a level set solver and ordered by their stratigraphical relationship.The stratigraphical relationship determines which material takes precedence over others when evaluating the level set functions.For example, if the level set values for material A and B are both negative at point x (meaning in the material domain for both A and B), A overlies B implies the material definition for point x is A instead of B. The overlying relationship comes from the geological model.Geologically, this overlying relationship is equivalent to the stratigraphic order (or the "younging direction") of the modeled lithological units.The details of the efficient implementation of the coupling between the level set and the Stokes solver is discussed in Text S3 of the supplementary material "Level set details" (Data Set S2).
To further reduce the computational cost, we notice that the level set values only enter the flow solver for cells within a narrow band of prescribed width around the zero level set.Therefore instead of solving Equation 8to steady state, we iterate up to τ max = γnh.Here we assume the transition narrow band is of a width of 2n cell and γ is a constant coefficient to compensate for the non-unity propagation speed due to the signum function smoothing (Karakus et al., 2018).Since too frequent reinitialization introduces excessive interface movement, we trigger the reinitialization at a fixed interval N reinit immediately prior to the adaptive mesh refinement of the cells close to the interface.

Deformation Analysis
The modeling approach introduced in this study takes into account the complete dynamics and enables the modeler to analyze both the instantaneous and finite deformation patterns of the simulated flow.This is achieved by applying the methodology of the kinematic structural analysis to the simulated flow field.This process is discussed in detail in this section.The interested reader is also referred to Jiang (2023, chapter 5, and references therein) for a more in-depth and pedagogical account of the flow kinematics analysis.
The instantaneous deformation field is readily available by analyzing the velocity field and the velocity gradient tensor   = ∕ .After decomposing the L into symmetric rate-of-deformation tensor D and anti-symmetric Flowchart showing the coupling between the finite element flow solver and the level set solver.The subscript q implies that the quantity is evaluated at the quadrature points.
vorticity tensor W, the following quantities can be evaluated at any arbitrary time and location in the numerical model: 1.The flow apophyses, or the "fabric attractors" (Passchier, 1998), can be obtained by taking the eigenvectors of L. The eigenvectors represent the direction where the flow is attracted to, assuming L is constant in time; 2. The instantaneous strain axis (ISA) can be obtained from the eigendecomposition of the rate-of-deformation tensor D; 3. The vorticity vector ω can be computed by ∇ × v. ω, which is the local sense of vorticity measured with respect to a fixed reference frame attached to the stationary mesh.We caution that this vorticity is not the same as the vorticity measured with respect to the local instantaneous principal strain rate axes and is therefore not objective (Jiang, 1999;Means et al., 1980).4. The kinematic vorticity number W k can be obtained by  ‖ ‖∕‖‖ .W k measures a nonlinear ratio between the rotation component and the stretching component of the deformation (Tikoff & Fossen, 1995).In a plane strain scenario, W k = 0 and W k = 1 represent two end-member flow types.The former corresponds to pure shear and the latter to simple shear.Differentiating flow types based on W k in 3d deformation is not straightforward, in which case W k = 1 can only be treated as a necessary condition for simple shear flow (Jiang, 2023).
The above-mentioned quantities represent the instantaneous condition of the flow.However, the deformation structures observed in the natural geological materials are commonly the results of accumulated deformation through a prolonged and complicated tectonic history.This process can be modeled by involving a Lagrangian tracking mechanism on top of the stationary mesh.In order to integrate finite deformation through time, tracer particles are injected into and around the regions of the model where deformation tracking is required.The particles are then advected in the modeling domain by the velocity interpolated from the mesh.The particle advection is computed by a forward Euler scheme as soon as the solution of the level set advection problem is completed.
To track finite deformation, each particle is associated with a deformation gradient tensor F(t), which is initialized to identity matrix I upon generation.For time step Δt n = t n+1 − t n , the approximated deformation tensor F n+1 can be integrated with the first-order forward Euler scheme of the well-known relation  Ḟ =  (Malvern, 1969, p. 163): The power of the deformation gradient tensor lies in its potential correlation with the foliation and lineation orientations, thus building a bridge between the abstract numerical model and the concrete field observations.The foliation and lineation orientation can be extracted by computing the eigendecomposition of the left Cauchy-Green tensor C = FF ⊤ , which in 3d space yields three positive eigenvalues (λ 1 ≥ λ 2 ≥ λ 3 ) and three corresponding mutually orthogonal eigenvectors (ν 1 , ν 2 , ν 3 ).The trend and plunge of the eigenvector ν 1 is the maximum stretching orientation of the finite strain and can be interpreted as the orientation of the stretching lineation.In the same way, the orientation of ν 3 corresponds to maximum finite shortening and can be interpreted as the pole to the foliation plane (e.g., Lin et al., 1998).Geometrically, the eigenvectors (ν 1 , ν 2 , ν 3 ) also correspond to the orientations of the three axes of the strain ellipsoid (Elliott, 1972).The square root of the eigenvalues ( ) are the length of the maximum, intermediate and minimum principal axes of the strain ellipsoid.Equipped with this knowledge, we can also compute the slope of any tracer particle on the Flinn diagram, or the Flinn's K-value: Values of K significantly greater than 1 are dominated by prolate strain or "cigar-shaped" strain ellipsoids.These particles mark the locations which favors the development of L-tectonites.On the contrary, particles with K significantly smaller than 1 are dominated by oblate strain or "pancake-shaped" strain ellipsoids, which correspond to the places where the development of S-tectonites is expected.The strain intensity N can also be computed by evaluating the deviation of the strain ellipsoids from a sphere (J.G. Ramsay & Huber, 1983): Journal of Geophysical Research: Solid Earth WÚ ET AL. 10.1029/2023JB026806 9 of 27 These predictions can be compared with the field observations of the deformation style and the relative intensity of deformation for validation.In addition to the spatial distribution, the numerical model provides the continuous temporal evolution of the deformation pattern, e.g. the change of the foliation orientation through time.This information can also be compared to the results of the analysis of strain history typically obtained from the field analysis of overprinting relationships between multiple generations of structures.

Numerical Benchmarks and Applications
Three numerical experiments are presented to test the validity (Sections 3.1 and 3.2) and showcase the real-world application of our technique (Section 3.3).We begin by computing two non-dimensional problems representing flow under simple shearing as well as density overturn.In 2d they serve as benchmarks with established numerical or near-analytical solutions so that we can verify the correctness of our results, evaluate the convergence of our numerical discretization, and establish the robustness of our solution technique for the linear system.We then extend the 2d benchmark cases to 3d to demonstrate the capability of our methods to characterize finite deformation patterns.Despite the non-dimensional nature of these numerical models, the results can be generalized to geological structures on multiple scales and in various tectonic settings as long as the governing equations and the linearly viscous rheology hold.
The reinterpretation of a non-dimensional variable (denoted with an overline, e.g.,   ,   ) to a corresponding variable with a physical unit (e.g., η, v) can be carried out in the following manner: let l c and t c denote the characteristic length and temporal scales such that   = ∕ and   = ∕ .The non-dimensional viscosity   can then be related to the physical viscosity η by   = ∕ .The non-dimensional velocity   and pressure   can then be scaled to their real-world counterparts by   =  ∕ and  =   .We remark that this scaling analysis is only possible due to the linear nature of our problem and is not generally applicable under more complicated material models.
Finally, to demonstrate the geological application of our method to crustal-scale tectonic processes, we model a density overturn process operating under prescribed background shearing in Section 3.3.This model represents a possible scenario for the Neoarchean crustal dynamics when vertical tectonism in the form of sagduction and diapirism as well as horizontal tectonism in the form of modern-day-like plate tectonics operates simultaneously.In terms of the physical process, this experiment demonstrate the 3d interplay between a Rayleigh-Taylor instability discussed in Section 3.2 and the simultaneous deformation of the rheologically inhomogeneous crust under the simple shear flow akin to Section 3.1.
Besides the above-mentioned numerical experiments, we have also included our results for an additional set of benchmark problems in the supplementary materials "Additional benchmark problems" (Data Set S3) and "Linear system solver details" (Data Set S1).

Competent Inclusion in Simple Shear Flow
In this set of experiments we simulate the deformation field produced by a relatively competent spherical inclusion (with viscosity η p ) embedded in a viscous matrix (with viscosity η m ) under a simple shear regime in both 2d and 3d.This model can be viewed as an abstraction for describing the deformation field in and around a competent lithological inhomogeneity undergoing shear deformation on a larger scale.On micro-or mesoscale, such matrix-inclusion system emerges naturally from the study of asymmetric porphyroclasts as shear sense indicators (Q.Zhang & Fossen, 2020), as kinematic vorticity gauges of the bulk flow (Stahr & Law, 2011) and as a proxy for modeling fabric development (Marques, 2016).On a regional scale, the model encapsulates, for example, the deformation field of less competent stratigraphy layers around a mechanically rigid pluton, both undergoing regional tectonic shearing.
Previous modeling of the matrix-inclusion system are typically based on one or more of the following assumptions: 1.The inclusion is considered as a rigid or deformable ellipsoid immersed in a infinitely viscous matrix.In the case of a rigid ellipsoid, the motion of the inclusion can be resolved in both 2d and 3d by applying Jeffery's theory (Jeffery, 1922) either analytically (e.g., Ghosh & Ramberg, 1976) or numerically (e.g., Jeẑek, 1994;Ježek et al., 1999;Jiang, 2007).If the ellipsoidal inclusion is considered deformable, the deformation of the inclusion can be obtained by Eshelby's theory for both linearly (Fletcher, 2009;Jiang, 2007) and power-law viscous rheology (Jiang, 2014).2. The flow is assumed to be steady-state.This assumption is made in several one-step finite element modeling on a bounded domain.The inclusion can be modeled as a separate phase more viscous than the matrix (Bons et al., 1997) or as a domain boundary (Masuda & Mizuno, 1996;Pennacchioni et al., 2000).
Attempts to overcome these limitations are made in 2d by combining a front-tracking algorithm with viscoplastic fast Fourier transform algorithm (Griera et al., 2013) or by a Lagrangian finite element method with remeshing (Schmid & Podladchikov, 2005).By computing two sets of numerical experiments with our formulation, we show that first, the results of our approach is comparable to previous studies in 2d and second, our approach provides a method to relax the aforementioned assumptions in 3d even for large amount of deformation.

2d Benchmark
In the first set of numerical experiments we reproduce the model setup by Bons et al. (1997).A steady-state flow field is computed for a circular inclusion with radius of 1/24 centered at the origin of a [−0.5, 0.5] 2 square box.
The matrix viscosity, denoted as η m , is set at 1 and the inclusion viscosity, denoted as η p , is set at 100.
The experiment is computed under three different BC specified in Bons et al. (1997, Section 2.1), which we briefly summarize as: 1. On all boundaries, v y = 0 and v x varies linearly from 0.5 on the top boundary (y = 1) to −0.5 at the bottom boundary (y = 0) 2. v x = 0.5(−0.5)on the top (bottom) boundary and v y = 0.0 on the left and right boundaries.Normal stresses on all boundaries are set to be equal to the background pressure.3. v x = 0.5(−0.5)and v y = 0 on the top (bottom) boundaries.On the left and right boundaries, v y = 0 and the normal stresses is equal to background pressure.
The computation is carried out on a base grid of 512 × 512.The grid is refined up to an equivalent density of 8,192 × 8,192 at the inclusion-matrix interface for accurate resolution of interface pressure.
To qualitatively compare the stream function ψ, y-velocity and pressure distributions of our results with that of Bons et al. (1997), Figure 3 is constructed to repeat the Figures 4a, 4d and 4g of Bons et al. (1997).The stream function values ψ at mesh nodes throughout the model can be obtained with the following equation, assuming the origin (x 0 , y 0 ) = (0, 0) to be the datum: The integral is computed numerically by a Gauss-Legendre quadrature along the cell faces outward from the datum.The contour patterns of the stream function ψ and y-velocity from the two studies show good agreement.The flow patterns of both studies are also mutually consistent: For BC 1 and 3 a bow-tie-shaped flow pattern is produced, while for boundary condition 2 the flow pattern is eye-shaped.1.Both the geometry of the separatices and the numerical values of ψ, D sp and D s are comparable to the results shown in Figure 4 and Table 1 of Bons et al. (1997).
Finally, the maximum absolute pressure value |p| max of our model is obtained at the inclusion-matrix interface (listed in the last row of Table 1) and can be compared with the analytical solution |p a | max (Schmid & Podladchikov, 2003): where    and   are the boundary pure strain rate and the simple shear strain rate.For this benchmark,    = 0 and   = 1 , thus yielding |p a | max = 1.9604.We remark that the numerical model is computed on a finite domain whereas the analytical solution assumes that the BC are imposed at infinity.Despite this discrepancy, the |p| max for all BC agrees well with |p a | max because the diameter of the inclusion is sufficiently small relative to the domain size.

3d Experiment
In the second numerical experiment we investigate the evolution of the finite deformation through time of a inclusion-matrix system in a 3d box of dimensions [−2, 2] × [−0.5, 0.5] × [−0.5, 0.5].The y-and z-components of the velocity are prescribed to be zero on all boundaries as a 3d non-steady extension to the previous 2d steady-state inclusion problem.The box is subject to dextral shearing by prescribing the x-velocity to be 1.0 on the top boundary and −1.0 on the bottom boundary.The x-velocity varies linearly with respect to the z-coordinate for all other boundary faces analogous to the boundary condition 1 in Section 3.1.1.The matrix viscosity η m is set to 1. Two numerical models are computed, with inclusion viscosity η p is set to a "more competent" value of 100 for one model and a "less competent" value of 10 for the other.The setup of the model is summarized in Figure 4.
Tracking particles are generated within and around the inclusion for tracking the finite deformation.Both models are computed up to a shear strain γ = 8.
The snapshots of our results are shown in Figure 5 and the animations of the model run are given in Movie S1 for "more competent" inclusion and Movie S2 for "less competent" inclusion.The modeling results show contrasting long-term deformation patterns for different inclusion viscosities.The results show that the behavior of both inclusion-matrix systems are consistent with the geological understanding, which we describe as follows.
In the case of a "more competent" inclusion (η p /η m = 100), the inclusion rotates in the flow and shows minimal internal deformation, as shown in the minimal change of the inclusion geometry and the almost spherical strain ellipsoids within the inclusion (Figure 5e).At the incipient stage of the shearing, the matrix material surrounding the inclusion is advected as well as deformed.The combined pattern of the particle clustering and the strain ellipsoid orientations show σ-type asymmetry.As the shearing progresses and more shear strain has accumulated, the bulk rotation of the inclusion exerts drag on the matrix material and the asymmetry gradually evolves into a In the case of a "less competent" inclusion (η p /η m = 10), the inclusion not only rotates, but at the same time undergoes elongation approximately parallel to the instantaneous stretching axis of the simple shear velocity field and shortening perpendicular to the instantaneous shortening axis.Prolonged deformation accentuates the ellipticity of inclusion and rotates the inclusion until its long axis is almost aligned with the shear direction (x-direction) (Figure 5j).Eventually the orientation of the inclusion is stabilized as it continues to undergo deformation by shear flattening.The matrix material around the inclusion shows a σ-type tail at the incipient stage of the model.However, as the rotation of the inclusion gradually halts toward the stable configuration, the matrix material around the inclusion experiences mostly shear-parallel stretching and minimal rotational drag.As a result the tail remains σ-type and does not evolve into a δ-type.
We first compute the standard 2d benchmark problem in Section 3.2.1 to show that our results conform with published results with the proper order of convergence so as to demonstrate our formulation produces correct results for non-steady tectonic flow with multiple materials.In Section 3.2.2we compute a 3d Rayleigh-Taylor instability experiment with tracking particles inserted at the material interface to showcase the capability of our method to visualize the finite strain distribution at different stages of the model evolution.

2d Benchmark
In this experiment we compute a viscous gravitational instability induced by an initially inverted density profile.The modeling domain is set up as in a rectangular box of height 1.0 and width λ = 0.9142.The width λ is chosen to maximize the growth rate of an isoviscous diapir (van Keken et al., 1997).In the initial setup a thin layer of buoyant material with density ρ 0 = 1,000 and thickness 0.2 is placed beneath a layer of denser material with density ρ 1 = 1,010 and viscosity η 1 = 100.A sinusoidal perturbation in the form of  = 0.02 cos( ∕) is applied to the material interface.The top and bottom boundaries of the box are considered to be solid walls (v = 0).A free-slip boundary condition is prescribed for the left and right boundaries.In two separate experiments, the viscosity of the buoyant layer η 0 is prescribed to be 100 (the isoviscous case) and 10 (the non-isoviscous case).The initial configuration of the model is shown in Figure 6.The models are computed on a base grid of 8 × 8 and the mesh is successively refined near the interface up to a mesh density equivalent to N max × N max around the material interface.We also compute the volume-normalized velocity norm (root mean square velocity) to facilitate comparisons between the published results: For the isoviscous case, the temporal evolution of the    is plotted in Figure 6 and the level set scalar fields at t = 500, t = 1,000, and t = 1,500 are shown in Figure 7.The velocity-time graph agrees in general trend with previous calculations (e.g., Hillebrand et al., 2014;Leng & Zhong, 2011).The results on the 128 × 128 and 256 × 256 mesh is almost identical, while the 64 × 64 mesh substantially deviates from the solution on finer meshes at the second peak.The deviation is attributed to the poor resolution of the interface geometry by the coarse mesh, which results in excessive interface smoothing and hinders the time at which the second gravitational instability event takes place.
To show that our modeling framework produces valid physics, we compare the value of the first peak root-mean-square velocity   max and the time t max at which the peak velocity is attained with established numerical results.A synthesis of published results and our data for both isoviscous and non-isoviscous cases are summa rized in Tables 2 and 3.The results show that for both the isoviscous and the non-isoviscous case, the DG level set formulation yields comparable results to previous studies obtained under a different discretization scheme regardless of whether the particle correction is applied or not.The results are also comparable to that obtained from a PIC method where the material information is entirely tracked by the Lagrangian particles.A convergence test is applied to the results by testing the speed at which the relative difference in   max and t max under successive mesh density contracts (Shadpour et al., 2015).The relative difference in   max and t max , denoted respectively as  ̄ and δ t can be computed with the following expression: Tracer-in-cell (Tackley & King, 2003) 64 × 64 209.9 0.003098 Level set (Hillebrand et al., 2014)  Note.The mesh density for this study is reported as N max × N max .The tracer density is 40 per cell for Tackley and King (2003) and 25 per cell for Samuel and Evonuk (2010).

Table 2
Benchmark Results for the Isoviscous Case (η 0 = η 1 = 100) Building upon the relative difference we can then compute the rates of contraction Λ of  ̄ and δ t , as follows: The convergence rate, shown in Table 4, reveals that both δ t and  ̄ contracts at the same rate with respect to mesh density.Hence, our numerical formulation is overall linearly convergent.This is expected from the fact that we use   1 basis for the DG level set solver.This result is also in agreement with the linear rate of convergence observed in the Zalesak's disk benchmark in the supplementary material.The result offers assurance that the application of the limiter does not degrade the overall rate of convergence of the numerical scheme.In summary, our results confirm that the level set formulation under DG discretization is numerically convergent and as effective as other methods for modeling multi-material geodynamic flow.

3d Experiment
We compute a 3d Rayleigh-Taylor instability to showcase the capability of our numerical scheme to track finite deformation by extending the previous 2d Rayleigh-Taylor benchmark to 3d.The model is constructed in a [0, 1] 3 cube with a sinusoidally perturbed material interface described as: Tracer-in-cell (Tackley & King, 2003) 64 × 64 73.5 0.009285 Note.The mesh density for this study is reported as N max × N max .The tracer density is 40 per cell for Tackley and King (2003) and 25 per cell for Samuel and Evonuk (2010).

Table 3
Benchmark Results for the Non-Isoviscous Case (η 0 = 10, η 1 = 100) where  ( ) = √ ( − 0.5) 2 + ( − 0.5) 2 is the distance to the center of the xy-plane.The BCs are all free slip except for the top and bottom, where they are considered as solid walls (v = 0).The material distribution closely resembles that of the isoviscous case of the 2d Rayleigh-Taylor benchmark.A density of 1,010 is prescribed to the material above the material interface given in Equation 19 and 1,000 to that below the interface.The viscosity is prescribed to be 100 throughout the domain.
Three snapshots of the model up to t = 200 are shown in Figure 8.In these snapshots, in addition to the geometry of the rising diapir, the strain ellipsoids are also shown at the material interface colored by their Flinn's K-value.A red color represents a stretching strain while a blue color is used for flattening strain.In Figure 9, the strain ellipsoids are also plotted on the Flinn diagram with respect to their depths (i.e., the z-coordinates).At t = 0 the strain ellipsoids are strain-free and therefore spherical.As the gravitational instability evolves, the areas at the "neck" and the bottom part of the model experiences constrictional strain in the vertical and horizontal directions due to the upward drag by the rising of the buoyant material situated at the lower part of the model.At the same time, the head of the diapir experiences compression at the material interface as the upper part of the diapir expands outwards and push the surrounding denser material sideways.In Figure 9 the same trend is revealed as the particles at deeper levels in the model show predominantly stretching strain and therefore favorable for developing lineation, while the particles at the diapir head preferentially develops strong foliation parallel to the material interface.Both lineation and foliation patterns are radially symmetric.The foliation is shallowly dipping in the center of the diapir and gradually steepens toward the periphery of the diapir.

Sheared 3d Rayleigh-Taylor Instability
The Rayleigh-Taylor experiments in Section 3.2 describe the idealized gravity-propelled density inversion processes where the material motion is predominantly vertical.Indeed, the structural pattern is overall consistent with that observed in the Archean granitoid-greenstone terranes where the partial density overturn is interpreted to be the dominant mechanism shaping the crustal architecture (e.g., Collins, 1989).However, in a realistic situation the vertical tectonic motion takes place under the influence of a background tectonic stress with a horizontal component.For example, in the Neoarchean terranes, evidence for horizontal and vertical tectonics coexist, and it is proposed that crustal scale density inversion may have operated simultaneously with horizontal tectonic strain (Chardon et al., 2002;Lin, 2005).The latter can be the result of far-field tectonic stress transmitted from the oblique convergence in a modern-day-like plate tectonics processes.In this experiment we construct a numerical experiment to simulate the simultaneous operation of the Rayleigh-Taylor instability and a horizontal simple shearing imposed as a boundary condition.Our primary purpose here is not to test interpretations on Archean tectonics, but to demonstrate that our modeling method can be applied to tectonic problems on a prescribed scale and can produce structural predictions readily comparable to field data, which would otherwise be impossible for 2d numerical models and prohibitively complicated for analytical solution.The bedding information and the interface geometry is difficult to obtain straight-forwardly with a PIC modeling approach.A description of the model results under different viscosity and BC, as well as an in-depth discussion of the implications for Archean tectonics is beyond the scope of this paper.
The model setup is illustrated in Figure 10.A box-shaped domain of dimensions 240 × 120 × 24 km is constructed.The geological model consists of a dense and competent supracrustal "greenstone" unit overlying a more buoyant and less viscous granitoid substratum.The upper layer and lower layer has a constant density of 2.8 × 10 3 kg m −3 and 2.7 × 10 3 kg m −3 , respectively.Both units are horizontally layered except for a perturbation of radius 10 km and height 1 km placed at the center of the boundary between two units.Due to positive buoyancy of the granitoid intrusion, a density inversion will be initiated.The viscosity of the supracrustal assemblages η 1 is set to 1 × 10 20 Pa s to represent the typical conditions of the upper crust.The viscosity of the granitoid substratum (η 0 ) is set to be one order of magnitude lower than the supracrustal assemblage (1 × 10 19 Pa s) to simulate a relatively "soft" middle to lower crustal material granitoid rocks of typical tonalite-trondhjemite-granodiorite (TTG) composition.
The BC of the model are set as follows.The top and bottom boundaries are both free-slip boundaries.The northern and southern boundaries have prescribed x-velocity of value +v s and −v s , respectively.The eastern and western boundaries are periodic boundaries with a reflection with respect to xz-plane model to be consistent with a monoclinic asymmetry due to the simple shear BC.In our experiment the magnitude of the shear velocity v s is set to 0.5 cm yr −1 , equivalent to a crustal shear strain rate of 8.3 × 10 −8 yr −1 which is on the same order of typical background strain rate of the continental crust under modern day plate tectonics (Xiong et al., 2021).
As the model evolves, the denser supracrustal assemblage sagducts into the deeper crustal level forming synclinal keels while the granitoid substratum upwells into the higher crustal levels and expands laterally.At the same time, the geometry of the granitoid intrusion is affected by sinistral shear and demonstrate an elongated geometry in the NW-SE direction due to the shear-induced elongation.This is in contrast with the radial symmetric diapir geometry in the 3d Rayleigh-Taylor instability without the simple shear boundary condition (Figure 8).To analyze the strain distribution, at the time (t max = 11.8Ma) when the peak root-mean-square velocity (   max = 1.7585 × 10 −8 cm ⋅ yr −1 ) is attained, we show the distribution of kinematic vorticity and the magnitude-scaled vorticity vector in Figure 11a, the strain ellipsoids colored by Flinn's K-value in Figure 11b and the strain intensity distribution N in Figure 11c.An animation showing the evolution of the model from the initial condition to t max is available in Movie S3.The kinematic vorticity shows that strong non-coaxial rotation is concentrated at the granitoid-supracrustal assemblage boundary, while in the sagducted keels of the supracrustal rocks the strain is mainly coaxial compression or stretching.This pattern is also reflected in the finite strain distribution, where the strain ellipsoids in the greenstone keels mainly show a flattening geometry.The vorticity vectors (Figure 11a) clearly demonstrate granitoid-up/greenstone-down sense of shear at the granitoid-greenstone contact.Because the density inversion is a quick process, large amount of strain is accumulated at the granitoid-greenstone boundary (Figure 11c) as manifested by the high value of N.
Although the results obtained from the numerical model is in 3d, we acknowledge that the data collected from the field survey typically resides on a 2d plane representing the present day ground surface, above which the crustal material is eroded away.Assuming the erosion level is known from the knowledge of the regional geology, we can then construct a slice plane of the 3d model at that depth.The lithological distribution and the structural pattern on that 2d plane is directly comparable to field survey data.To demonstrate this process, a slice surface is constructed at a depth of 8 km of the 3d model.The depth of the slice plane is so chosen because most Archean greenstone terranes display a greenschist facies metamorphism.On top of the lithological pattern, orientation of the primary structures (S 0 , also see Section 4.1), the foliation and the lineation are respectively shown in Figures 12a-12c along with the lithological pattern.A simulated geological cross-section approximately perpendicular to the regional NW-SE structural grain is also constructed in Figure 12d.Local lineation and foliation patterns in a greenstone keel SW of the central granitoid dome are shown in the stereonets in Figures 12e and 12f.The primary structures, which can be interpreted as the bedding in the supracrustal assemblages and pre-tectonic subhorizontal foliation in the granitoid substratum, shows that the greenstone keel corresponds to regional syncline while the granitoid dome forms the core of the regional anticline.
The finite strain patterns in the keels of the supracrustal rocks west-southwest of the granitoid dome is analyzed in detail (with the approximate location shown in Figure 12a).The stereographic projection of the poles to the foliation (Figure 12e) shows that the foliation in the keel is mostly steeply-dipping and strike in a NW-SE direction subparallel to the adjacent granitoid-greenstone boundary.The lineation orientation (Figure 12f) is scattered along the great circle defined by the local average foliation orientation and shows two point maxima that are moderately to steeply plunging.The moderate to steep plunge angle of the lineation is likely the result of the vertical sagduction movement.Meanwhile, the horizontal shearing component contributes to the scattering of  This structural pattern is comparable to that observed in Neoarchean greenstone belts where both horizontal and vertical tectonic processes are postulated to operate synchronously on Earth (Lin, 2005;Lin et al., 2013), detailed discussion of which is beyond the scope of this contribution.

Geological Interpretation of the Level Set Field
Besides being an identifier of material domains, the level set field in our modeling approach contains geometric information that can be helpful for interpreting the modeled deformation patterns.For example, in a crustal-scale model where the initial layers of material domains correspond to the pre-deformation stratigraphic units, the orientation of the level set isosurfaces can be interpreted as the orientation of the primary structures, or S 0 (e.g., sedimentary beddings).The deformation-induced curving and distorting of the isosurface can be interpreted as the regional folding of the beddings of the stratigraphy units.The orientation of the bedding can be extracted by taking the gradient of the level set function (∇ϕ).Assuming the level of erosion is available from the knowledge of the local geology (e.g., mineral assemblages and or thermodynamic modeling), a 2d cut surface at that erosional level can be constructed within the 3d model which shows the lithological distribution and the bedding orientation on that particular 2d surface.Such data can be compared with the present-day geological map pattern and field bedding measurement as a validation for the modeled results.Moreover, the level set field itself represents the distance to the material interface, therefore making our approach suitable for exploring the relationship between variations of structural patterns with respect to distance to the lithological boundaries.
As a side note, one can also extract the curvature of the material interface by taking the Laplacian of the level set field (∇ 2 ϕ).The boundary curvature is related to the boundary mobility and therefore controls the rate of grain boundary migration in microstructure simulations (Chapter 5 Humphreys & Hatherly, 2004).

Geological Interpretation on Deformation Analysis
We remark that compared to the pure kinematic modeling approaches which frequently utilizes the assumption of a steady-state velocity gradient L, our formulation has the ability to incorporate non-steady and inhomogeneous deformation states by allowing L to vary naturally as the strain state changes with the material flow.The finite deformation field is the direct integration of the instantaneous deformation field as the tracer particles are advected through the model.In this way by incorporating the full dynamics, our approach simulates deformation patterns in a more general sense than the traditional kinematic modeling.In fact, the latter can be considered as a steady-state approximation in a small enough interval in both space and time from the viewpoint of a model with complete dynamics.
The foliation and lineation patterns obtained from our modeling formulation are based upon the assumption that the lineation and the pole to foliation always follow the principal strain axes of the local strain ellipsoids.However we acknowledge that this assumption is by itself a simplification because the alignment of the lineationor foliation-defining mineral grains does not necessarily track the principal axes of the local strain ellipsoid in a precise and predictive manner.To probe into the microstructural evolution of the deformation fabric at each tracer particle, our modeling method framework can couple with an efficient solver based on micromechanics to model such processes on a finer scale (Jiang & Bentley, 2012).Furthermore, the information about the strain ellipsoids can not only be compared to the field-observable alignment of minerals, it can also be extended to strain ellipsoid geometry obtained through other techniques e.g. the anisotropy of magnetic susceptibility measurements (Benn, 1994).
Because the deformation field that the numerical model takes into account is not steady-state, the structural patterns recorded on the particles are dependent upon the interval in time through which the deformation is recorded.Taking the quartz crystals deformed in a ductile shear zone as an example, in the temperature regime of subgrain rotation recrystallization, large quartz porphyroclasts are transformed into ribbon grains and meanwhile newly born recrystallized grains will emerge (Stipp et al., 2002).If we spread particles evenly throughout the domain, then the particle corresponding to the recrystallized grain will have to reset its F when the recrystallization takes place, while those corresponding with the ribbon porphyroclast will keep recording F for a longer period of time (Means, 1981).The same situation also applies to geological bodies on outcrop or even regional scale (e.g., boudinage of veins and dykes).Unfortunately, a widely accepted criteria as to the choice of this time interval is problem dependent and the modeler is cautioned to take this into consideration when interpreting the modeled structural pattern.
We also caution that the vorticity-related quantities (e.g., kinematic vorticity W k ) herein obtained are computed on a stationary Eulerian coordinate and therefore not objective.Direct usage of such quantities for determining flow classification and understanding fabric development may be invalid unless corrections are made for external vorticity caused by the spin of the ISA (Astarita, 1979;Means et al., 1980).One algorithm for implementing such corrections using the ISA as a reference is proposed by Jiang (1999).

Summary and Future Extensions
In this study we have developed a modeling scheme by coupling a level set solver with DG discretization for capturing material interface and a finite element solver for resolving flow dynamics.A combination of Eulerian mesh and Lagrangian tracer particles allows us to model large deformation without sacrificing the capability to track deformation history.A complete set of kinematic information from the instantaneous flow field as well as the finite deformation field can be extracted from the numerical model, making it well-suited for analyzing spatially and temporally variant structural responses as a result of the processes on a larger scale and the partitioning of larger-scale tectonic movement into material domains with contrasting mechanical properties.Our model technique also provides the possibility for ground-truthing tectonic models by comparing field observations of the geological map patterns, the strain distribution and the structural patterns with model predictions of material domain geometries and deformation patterns, thus building the bridge between abstract numerical models and concrete field evidence.
Our implementation serves as a proof of concept of integrating DG-discretized level set method into a modeling framework for tectonic flow and is therefore simplified in many aspects.However, this modeling method can be further developed in many ways, making it possible to simulate real world tectonic problems with ever more challenging physical complexities.Future improvements and extensions of the modeling framework will proceed in the following directions: 1. Implementation of more realistic rheological models.This includes, for example, a nonlinear (power-law) viscosity model for dislocation creep processes, evolving anisotropic mechanical properties due to the preferred alignment of minerals in the flow (e.g., Lev & Hager, 2008;Mühlhaus et al., 2003;Ran et al., 2019) and temperature-dependent viscosity model by coupling with a convection-diffusion solver for heat transport.2. Communicating the deformation field from our tectonic modeling approach to an efficient analytical microstructural solver for fabric development and lattice preferred orientations (e.g., Jiang & Bentley, 2012;Jiang, 2016).With this coupling scheme, a direct correlation of deformations across multiple scales can be made possible.

Figure 2 .
Figure 2.Flowchart showing the coupling between the finite element flow solver and the level set solver.The subscript q implies that the quantity is evaluated at the quadrature points.

Figure 3 .
Figure3.The pressure and the velocity field around a spherical inclusion at the center of the modeling domain.Our model reproduces the results ofBons et al. (1997) for linearly viscous flow.The subfigures (a), (b), and (c) correspond to the boundary conditions 1, 2, and 3 (i.e., the Figures4a, 4d, and 4g) inBons et al. (1997).On the left half of each subfigure, the white contours are streamlines at ψ values of 0.0005, 0.001 0.0025, 0.005, and 0.01; the blue line is the flow separatrix and the color shading represents the pressure field.The right half of each subfigure is shaded by the value of y-velocity and the black lines are the contours of y-velocity at −0.015, −0.010, −0.005 and 0.0.See text for discussion.

Figure 4 .
Figure 4. Model setup of a 3d deformable inclusion under simple shear.Upper: the 3d shear box setup at t = 0.The deformable inclusion (shown in dark red) is centered at the origin of the shear box.Blue dots are the tracking particles generated around the inclusion.Lower: A cross-section across the center of the shear box which corresponds to the light blue plane in the upper diagram.

Figure 5 .
Figure 5. Snapshots of the deformation for a 3d deformable inclusion in a shear box up to a shear strain γ = 8 for η p /η m = 100 (subfigures (a) to (d)) and η p /η m = 10 (subfigures (f) to (i)).The tracking particles are colored by the viscosity value.The ellipsoids attached to the tracking particles are strain ellipsoids.The streamlines (in white) and pressure (represented by the color shading) are plotted on an xz slice plane through the origin.Subfigures (e) and (j) show the magnified views of the strain ellipsoids on a 2d cross-section through the origin on the xz-plane for η p /η m = 100 and 10 respectively.The strain ellipsoids are colored by the Flinn's K-value.See text for discussion.

Figure 6 .
Figure 6.(Left) Model setup of the 2d Rayleigh-Taylor benchmark.(Right) The evolution of    up to t = 2,000 for the isoviscous Rayleigh-Taylor instability on a grid with the finest mesh density of 64 × 64, 128 × 128, and 256 × 256.

Figure 7 .
Figure 7. Snapshots of the mesh configuration and the level set field of the 2d Rayleigh-Taylor instability benchmark at (a) t = 500, (b) t = 1,000, and (c) t = 1,500.The solid white line outlines the material interface (the zero level set).

Figure 8 .
Figure 8.The finite strain distribution of a 3d isoviscous Rayleigh-Taylor instability at t = 100 (left), 150 (middle) and 200 (right).The ellipsoidal markers are strain ellipsoids colored by Flinn slope K.The red color range implies constrictional strain, while the blue color range indicates flattening strain.

Figure 9 .
Figure 9.The Flinn diagram of the strain ellipsoids at the interface between the denser material on top and the lighter material below in the 3d isoviscous Rayleigh-Taylor instability test case.The model snapshots are taken at t = 100 (left), t = 150 (middle), and t = 200 (right).The marker points are colored by their z-coordinates: z = 0.0 corresponds to the deepest level (bottom) of the model, while z = 1.0 corresponds to the top surface of the model.The diagonal line in each subfigure is the line of plane strain.The area above the line of plane strain represents constrictional strain while the area below the line has a flattening type of strain.

Figure 10 .
Figure 10.Model setup for the two-layer 3d Rayleigh-Taylor experiment with a sinistral shearing boundary conditions (BC).v s is the magnitude of the prescribed boundary shear velocity.The eastern and western BC are periodic with a reflection such that points A, B, C, and D on the eastern boundary correspond to A′, B′, C′, and D′ on the western boundary respectively.
lineation pattern along the big circle trace and grouping of the lineation orientations into two point maxima.

Figure 11 .
Figure 11.Model results for the sheared 3d Rayleigh-Taylor instability at peak root mean square velocity.An animation for this experiment is provided in Movie S3.In panel (a), the faint white wireframe delineates the interface between the greenstone supracrustal assemblages and the granitoid substratum.The color shading corresponds to the value of the kinematic vorticity number W k value.The conical arrows are the orientations of the vorticity vectors ω scaled by their magnitude.The computation procedures of ω and W k are discussed in detail in the second paragraph of Section 2.3.(b) Shows the geometry of the strain ellipsoid colored by the Flinn's K-value computed by Equation 12; (c) shows the 3d distribution of the strain intensity N computed by Equation 13.

Figure 12 .
Figure 12.Simulated geological map with structural patterns of the sheared 3d Rayleigh-Taylor instability model on a horizontal slice at 8 km depth.The structural patterns shown in (a), (b), and (c) are respectively the primary structures (subhorizontal planar markers before deformation, e.g., bedding), the foliation and the lineation.(d) Is a geological cross-section along the line XX', the location of which is delineated in (a).The shear sense at the granitoid-greenstone boundary is extrapolated from the vorticity vector in Figure 11a.Because the simulated geological maps in (a), (b), and (c) are constructed at 8 km depth, the geological material above z = −8 km (shown in faint colors) is considered to have been removed due to erosion and no longer visible in the field.(e) and (f) are stereographic projections of the poles to foliation (⊥S) and lineation (L) measurements in the synclinal keel with location shown by the gray rectangle in (a).