LitMod2D_2.0: An Improved Integrated Geophysical‐Petrological Modeling Tool for the Physical Interpretation of Upper Mantle Anomalies

LitMod2D integrates geophysical and petrological data sets to produce the thermal, density, and seismic velocity structure of the lithosphere and upper mantle. We present a new LitMod2D_2.0 package with improvements focused on (i) updated anelastic attenuation correction for anharmonic seismic velocities, (ii) chemical composition in the sublithospheric mantle, and (iii) incorporation of sublithospheric mantle anomalies. Sublithospheric mantle anomalies can be defined with different chemical composition, temperature, seismic velocities, and a combination of them, allowing the application of LitMod2D_2.0 to regions affected by mantle upwelling, subduction, delamination, and metasomatism. We demonstrate the potential application of LitMod2D_2.0 to such regions and the sensitivity of thermal and compositional anomalies on density and seismic velocities through synthetic models. Results show nonlinearity between the sign of thermal and seismic velocity anomalies, and that S wave velocities are more sensitive to temperature whereas P wave velocities are to composition. In a synthetic example of subduction, we show the sensitivity of sublithospheric mantle anomalies associated with the slab and the corner flow on surface observables (elevation, geoid height, and gravity anomalies). A new open‐source graphic user interface is incorporated in the new package. The output of the code is simplified by writing only the relevant physical parameters (temperature, pressure, material type, density, and seismic velocities) to allow the user using predefined post‐processing codes from a toolbox (flexure, mineral assemblages, synthetic passive seismological data, and tomography) or designing new ones. We demonstrate a post‐processing example calculating synthetic seismic tomography, Rayleigh surface‐wave dispersion curves, and P wave receiver functions from the output file of LitMod2D_2.0.


Introduction
The characterization of the present day physical state (pressure, temperature, and chemical composition) and architecture of the lithosphere and sublithospheric upper mantle is a fundamental prerequisite for understanding their evolution and underpins large-scale exploration programs (cf. Afonso, Moorkamp, & Fullea, 2016;Hoggard et al., 2020). In this context, information about the density and temperature fields within the lithospheric and sublithospheric mantle typically comes from modeling of geophysical and/or geochemical data (e.g., Afonso et al., 2008;Deen et al., 2006;Fullea et al., 2007;Goes et al., 2000;Griffin et al., 1999;Lachenbruch & Morgan, 1990;Priestley & McKenzie, 2006;Ritzwoller et al., 2004;Zeyen et al., 2005). Although these approaches have provided a great deal of knowledge about the physical state of the lithosphere and shallow upper mantle, significant discrepancies in the predictions from these methods are common in the literature (e.g., Afonso et al., 2008).
Integrated geophysical-petrological modeling of available observables (i.e., elevation, gravity anomalies, geoid height, surface heat flow, seismic velocities, xenoliths, and xenocrysts data) is a particular approach that allows exploring and reconciling observations from different data sets and methods. In particular, LitMod2D_1.0 (Afonso et al., 2008) is a tool developed to study the thermal, compositional, density, and seismological structure of the lithospheric and sublithospheric mantle by combining information from petrology, mineral physics, and geophysical observables within a self-consistent thermodynamic framework. As such, this approach is well suited to (i) handle the intrinsic nonuniqueness associated with the modeling of single observables and (ii) distinguish thermal from compositional effects in geophysical signatures.
LitMod2D_1.0 has been successfully applied in various tectonic settings including continental margins (e.g., Fernàndez et al., 2010;Pedreira et al., 2015) and continental collision regions (e.g., Tunini et al., 2015). It allows considering different chemical compositional domains within the lithospheric mantle, but it assumes a homogeneous sublithospheric mantle with a primitive upper mantle (PUM) chemical composition (McDonough & Sun, 1995). The latter assumption limits the use of LitMod2D_1.0 to regions not affected by mantle upwelling, subduction, and/or delamination. These regions are often characterized by changes in temperature and/or chemical composition in the sublithospheric mantle, as suggested, e.g., by tomography models (e.g., Cammarano et al., 2009;Koulakov et al., 2016).
In this paper, we introduce LitMod2D_2.0, an updated version of the original 2D software by Afonso et al. (2008). Some key novel features of LitMod2D_2.0 include (a) a new open-source graphic user interface (GUI) developed in Python programming language (Hunter, 2007;Rossum, 1995) with improved functionalities to facilitate its updating and cross-platform use, (b) new functionalities for the correction of anharmonic seismic velocities to account for anelastic attenuation using recent laboratory-based measurements, and (c) the option to use different compositions in the sublithospheric mantle. In the following sections, we briefly discuss the general methodological framework and the updated chemical composition in the sublithospheric mantle. Then, we discuss the numerical implementation and effects of considering sublithospheric mantle anomalies of thermal, compositional, and seismic origin (Section 3). The new open-source GUI is discussed in Section 4. Input and output data have been simplified including a post-processing toolbox that users can use and further develop according to specific needs. The general applicability of the new LitMod2D_2.0 package is illustrated in Section 5. We assess the sensitivity of surface observables to sublithospheric mantle anomalies through a synthetic model of a subduction zone setting and compute Rayleigh surface-wave dispersion curves and synthetic P wave receiver functions to show the post-processing possibilities (Section 6).

Methodology and New Improvements
LitMod2D is a software package based on the finite element forward modeling approach of the CAGES code (Zeyen & Fernàndez, 1994). LitMod2D_1.0 incorporated an external GENERATOR module based on Perple_X subroutines (Connolly, 2005) later upgraded (Connolly, 2009) to calculate the mineral assemblages and their physical properties as functions of the bulk composition of mantle domains. It also incorporates a MATLAB-based GUI. For a detailed description of LitMod2D_1.0, the reader is referred to Afonso et al. (2008). Here, we provide a brief summary of the critical aspects of the method and focus on the implementation of the new features/improvements in LitMod2D_2.0.

General Workflow
The general workflow (Figure 1) in both versions is similar, except for the simplified input/output files, new improvements, and the new post-processing toolbox module. The numerical domain extends from the surface to 400 km depth and comprises different crustal and mantle bodies characterized by their individual thermo-physical properties and chemical composition. Crustal bodies are characterized by user-defined thermo-physical properties (e.g., thermal conductivity, volumetric heat production rate, coefficient of thermal expansion, and compressibility), with an option of depth and/or temperature dependence. The composition of upper mantle bodies is assigned within the Na 2 O-CaO-FeO-MgO-Al 2 O 3 -SiO 2 (NCFMAS) system. Stable mineral assemblages are calculated by the external module GENERATOR using a Gibbs free-energy minimization algorithm (Connolly, 2005(Connolly, , 2009 for pressure and temperature ranges in the upper mantle. Here, we use an augmented-modified version of Holland and Powell (1998) thermodynamic database (Afonso & Zlotnik, 2011). Relevant physical properties (density, thermal conductivity, and bulk and shear moduli) are calculated for each mineral and averaged according to Voigt-Reuss-Hill procedure (Hill, 1952) producing look-up tables as function of P and T. In this way, the user can produce its own library of look-up tables for different mantle compositions beforehand.
Temperature distribution is calculated under the assumption of steady-state by solving the heat transport equation using finite elements, with prescribed boundary conditions at the surface (e.g., 0°C) and at the lithosphere-asthenosphere boundary (LAB) (e.g., 1,320°C). Below the LAB, the algorithm considers a 40-km thick thermal buffer with a temperature of 1,400°C at its base to avoid unrealistic discontinuities between a conductive thermal gradient within the lithosphere and an adiabatic thermal gradient within the sublithosphere.
The thermal conductivity within the lithospheric mantle depends on temperature, pressure, and composition and is calculated by GENERATOR according to Grose and Afonso (2013) and read from the look-up tables during modeling. The temperature gradient below the thermal buffer layer is restricted to 0.35 ≤ dT/dz ≤ 0.50°C/km. The density distribution is obtained with an iterative scheme to include the effect of pressure, temperature, and composition. Pressure is calculated at every node as a function of the overburden lithostatic pressure, and it is then used to obtain the P-T-composition-dependent density at the node. In case of mantle bodies with a specific composition, density is read from the look-up tables produced by GENERATOR. This density is then used in the next iteration to update the overburden pressure and to solve for an updated The new version presented in this study may account for sublithospheric mantle anomalies in which case coupled elevation is calculated. It also has an option for post-processing of relevant physical properties (temperature, pressure, density, and seismic velocities) to compare with additional data sets. Flexural isostasy is included by coupling "tao-geo" software (Garcia-Castellanos et al., 2002). "Computer Programs in Seismology" (CPS, Herrmann, 2013) is also coupled to work with passive seismological data sets. Stable mineral assemblages are extracted from the material files produced by GENERATOR subprogram. density. The final density distribution is then used to calculate elevation (under the assumption of local isostatic equilibrium), gravity anomalies (Bouguer and free air), and geoid height. Anharmonic seismic velocities in the lithospheric and sublithospheric mantle are assigned according to stable mineral assemblages at the corresponding P-T conditions from the assigned look-up table.

Anelasticity
Accounting for anelasticity in the modeled seismic velocities is crucial for the physical interpretation of real seismic anomalies in terms of temperature and/or chemical composition (Abers et al., 2014;Afonso et al., 2008;Cammarano & Romanowicz, 2008;Goes et al., 2000;Lau & Faul, 2019;Sobolev et al., 1996;Takei, 2017). In the previous LitMod2D_1.0 version, the anelastic attenuation correction was incorporated a posteriori by using an external code. In LitMod2D_2.0, anelastic effects are incorporated inside the GENERATOR module using Equations 1 to 3 and the empirical parameters proposed by Jackson and Faul (2010) where we assume Q P = (9/4)Q S . V Po (T, P) and V So (T, P) are the anharmonic seismic velocities at a given temperature and pressure for a given chemical composition, and A = 816 s −α μm −α , α = 0.36 is the frequency dependence factor, E = 293 kJ mol −1 is the activation energy, V = 1.2 × 10 −5 m 3 mol −1 is the activation volume, and R the universal gas constant. Here, we use a grain size of d = 10 mm in the mantle and an oscillation period of T o = 75 s, since this combination produces Q P and Q S values in the range of global average attenuation models, particularly close to the LAB where anelastic attenuation is believed to be high ( Figure  2). Both grain size and oscillation period are input parameters that can be modified. Moreover, the user can incorporate her/his own preferred attenuation model into the GENERATOR module. Hereinafter, we use seismic velocities corrected for anelastic attenuation unless specified otherwise.
Anelastic attenuation parameters derived from laboratory measurements on dry polycrystalline olivine often need to be extrapolated to conditions pertaining to the upper mantle, thus introducing unquantifiable uncertainties (Faul & Jackson, 2015;Priestley & McKenzie, 2013). While it is hard to put hard uncertainty values on these parameters, here we assume 10% error and explore how these errors propagate to the calculated seismic velocities. Uncertainties in computed seismic velocities (V P and V S ) are calculated using: where ∂V P=S ∂α is the partial derivative of the velocities (Equations 1 and 2) with respect to α and δα is the assumed error in α; similar terms apply for E and V. Figure 3 shows the resulting errors in V P and V S . We see that these errors have the largest effect on seismic velocities around the LAB and that the dominant source of uncertainty is α affects the seismic velocities the most and V, activation volume, the least. Note that the error in velocity (Equation 4) is directly proportional to the error in anelastic parameters and therefore, increasing uncertainty in an individual parameter will scale the uncertainty in the velocity in a quasi-linear fashion. The effects of melt or water content on seismic velocities are not included here. If melt is thought to be present, an additional correction needs to be applied to the computed seismic velocities (e.g., Afonso, Rawlinson, et al., 2016).

Sublithospheric Mantle Composition and Recalibration of Elevation and Tomography
The chemical composition of the lithospheric mantle may be estimated from available mantle xenoliths and xenocrysts and/or exposed massifs in mobile belts. Although sublithospheric mantle rocks are less represented in these suites, compositional estimates can sometimes be obtained by analyzing the chemistry of primitive basaltic melts (when available; e.g., Shorttle & Maclennan, 2011;Brown et al., 2020). As mentioned in Section 1, LitMod2D_1.0 uses a fixed chemical composition for the sublithospheric mantle, corresponding to the PUM of McDonough and Sun (1995). The upper mantle, however, has been a source for oceanic and continental crust since the onset of plate tectonics and should be less fertile than PUM (Van Keken et al., 2002). Here, we test the most common sublithospheric mantle compositions proposed so far (Table 1) and compare their relevant physical properties with that of PUM. The depleted mid-ocean ridge basalt (MORB) mantle (DMM) is a source reservoir to MORBs and has been computed using a robust geochemical data set (trace elements) on abyssal peridotites (Workman & Hart, 2005). Other commonly used sublithospheric mantle compositions are PUM-3%_N_MORB (Workman & Hart, 2005) and pyrolite (Ringwood, 1977). PUM-3%_N_MORB is computed by extracting 3% of normal MORB from PUM composition, whereas pyrolite is calculated by mixing appropriate fractions of basalts (partial melts from the mantle) and peridotites (the presumed residues from partial melting).

Figure 2.
Depth distribution of anelastic attenuation parameter (Q, quality factor), for range of grain sizes (1-100 mm) and oscillation periods (10-100 s) compared with global average Q P and Q S models. (a) and (b) show depth distribution of Q P and Q S for varied grain size and constant oscillation period of T 0 = 75 s, respectively. (c) and (d) show depth distribution of Q P and Q S for varied oscillation period (10-100 s) and constant grain size of d = 10 mm, respectively. Q P and Q S from ak135 (Kennett et al., 1995) a global average model (continuous red solid line) and Qs from QL6c (continuous gray line, Durek and Ekström (1996)), are also plotted for comparison.
All tested compositions are less dense than PUM at sublithospheric mantle depths (Figure 4a), which is consistent with their higher Mg# (100 × MgO/[MgO + FeO]) and depleted reservoirs. DMM shows a small density change of −0.18% relative to PUM between 120 and 340 km depth and an increase of~0.1% below 340 km, whereas pyrolite shows a maximum change of approximately −0.6%. PUM-3%_N_MORB exhibits variations in between those associated with DMM and pyrolite. In terms of anharmonic P and S wave velocities (Figures 4b and 4c), pyrolite is the slowest, DMM is slightly faster in comparison to PUM, and PUM-3%_N_MORB is intermediate. A water content of 100 ppm (70-160 ppm; Workman & Hart, 2005) in DMM slightly decreases the density compared to its dry counterpart and produces a decrease in P waves velocities higher than in S wave velocities because of its lower bulk modulus (Watanabe, 1993).
Although PUM-3%_N_MORB and pyrolite attest to be depleted with respect to PUM, they still are theoretically computed. Furthermore, the pyrolite composition has been reported not to satisfy trace elements or isotopic data on basalts, violating the chondritic abundances (Anderson, 1989). DMM composition is consistent with trace elements data and with 6% aggregated fractional melting to produce an average 6-km thick oceanic crust at the mid-oceanic ridge (Klein & Langmuir, 1987;Workman & Hart, 2005). In addition, DMM has also been shown to be balanced by continental crust, ocean-island-basalt source, and oceanic crust, which are the products of mantle melting (Workman & Hart, 2005). Although the differences in physical properties  are small, in the new LitMod2D_2.0 version, we fix the sublithospheric mantle chemical composition corresponding to DMM for geochemical consistency.
Changing the chemical composition of the sublithospheric mantle requires a recalibration of the calculated absolute elevations since the total load of the sublithospheric mantle will vary according to the adopted composition. In contrast to LitMod2D_1.0, that considers a reference column at a mid-oceanic ridge (MOR) down to 400 km depth, in our case, we consider a reference column corresponding to a thermally stable oceanic lithosphere and the underlying sublithospheric mantle. Thermally stable oceanic lithosphere is chosen because the depth-dependent T-P-composition is arguably less complicated than at actively spreading MOR, where melt content and short-lived buoyancy sources can complicate its characterization. Four compositional layers have been considered in the lithosphere, and their bulk composition is calculated according to the depth-dependent melt fraction for a standard MOR model (Niu, 1997;Turcotte & Morgan, 2013) as used in Fernàndez et al. (2010). The calibration is done using the formulation described in Afonso et al. (2008) and Fullea et al. (2009) where we solve for a calibration parameter constrained by average elevation, heat flow, and lithospheric structure at stable oceanic lithosphere (for details, see Appendix A).
The calculated seismic velocities in the sublithospheric mantle also depend on chemical composition. Seismic tomography models are usually reported as deviations from global reference velocity models; one of the most widely used reference models is ak135 (Kennett et al., 1995). To compare the seismic velocities obtained from LitMod2D_2.0 with those from ak135, we consider a model with (i) a 35-km thick continental crust and (ii) an 85-km thick lithospheric mantle with Average Garnet Tecton composition (Tc_1 in Griffin et al., 2009), resulting in a LAB depth of 120 km. The composition of the sublithospheric mantle corresponds to DMM. Thermo-physical parameters used in the crust and composition in the lithospheric mantle are listed in Tables 2 and 1, respectively.
The ak135 model was designed to predict arrival times for seismic phases in observational seismology by inverting smoothed empirical travel times (Kennett, 2006;Kennett et al., 1995). P wave velocities are better constrained than S wave velocities since P waves are first arrivals, whereas S waves have low frequency and  (Table 1). (a) Resulting density, (b) P wave velocities, and (c) S wave velocities are compared with that to PUM.

Geochemistry, Geophysics, Geosystems
can be interfered by the CODA of P waves. Considering this and the inherent uncertainties in earthquake hypocenters, we assign uncertainties of 0.5% and 1% to P and S wave velocities, respectively. Different thermodynamic databases produce noticeably different seismic velocities, particularly in the sublithospheric mantle ( Figure 5) resulting in a dispersion of~1%. Formal errors resulting from each thermodynamic database can be calculated from the uncertainties in the elastic moduli of individual minerals. Both Stixrude and Lithgow-Bertelloni (2005) and Holland and Powell (1998) modified by Afonso and Zlotnik (2011) databases reasonably reproduce the ak135 model for depths between 35 and 250 km, whereas Xu et al. (2008) database results in very slow P wave velocities (Figure 5a). Below 250 km depth, only Stixrude and Lithgow-Bertelloni (2005) database shows good agreements with ak135 model, both in P and S velocities, whereas Holland and Powell (1998) modified by Afonso and Zlotnik (2011) and Xu et al. (2008) databases are significantly slower (~1-2%). Slow velocities below 250 km depth are also observed in Cammarano et al. (2005). Fitting the ak135 velocities below 250 km with these thermodynamic databases would require either a lower temperature than that predicted by the assumed adiabatic thermal gradient or a change in the chemical composition at these depths (Cammarano et al., 2009). Another possible contribution to the discrepancies between predicted and reference velocities could be that the temperature and pressure derivatives of elastic moduli in the Afonso and Zlotnik (2011) version of the Holland and Powell (1998) need to be updated. In contrast to the other two databases/formalisms, that of Afonso and Zlotnik (2011) supplements equilibrium phase diagrams with an independent database of velocity derivatives; small corrections to these derivatives can produce changes in the velocities of the same order as the discrepancies in Figure 5. Despite this, the phase equilibria predictions show excellent agreement with experimental data (Afonso & Zlotnik, 2011).
Here, we choose the Holland and Powell (1998) modified by Afonso and Zlotnik (2011) database in LitMod2D_2.0 to calculate the physical properties in the upper mantle given its performance in reproducing experimental data, but we note that a future reassessment of the above derivatives is necessary. Since deviations with respect to the ak135 model are always negative below~250 km (in the case that there are no thermal or compositional perturbations at these depths), we avoid this systematic misfit by using the above reference model (35-km thick crust and 120-km LAB depth; Figure 5) to correct our predictions when calculating synthetic tomography results.  (Kennett et al., 1995), in solid gray circles, with an error of 0.5% in P wave velocities and 1% in S wave velocities are plotted for comparison.

Sublithospheric Mantle Anomalies
A major contribution of LitMod2D_2.0 is the possibility to incorporate sublithospheric mantle anomalies into the modeling workflow, thus opening up the possibility of modeling complex geodynamic processes, such as mantle upwelling, subduction, delamination, and metasomatism (i.e., processes that can easily modify the T and/or chemical composition beneath the LAB). Since seismic tomography has the potential to detect such perturbations in terms of fast/slow velocities, we added in LitMod2D_2.0 a functionality to incorporate (i) thermal anomalies, (ii) compositional anomalies, and (iii) seismic velocity (V P , V S ) anomalies, in the sublithospheric mantle.
In this section, we discuss the numerical implementation of sublithospheric mantle anomalies and its applicability to various geological settings considering the reference model defined for calculating synthetic tomography ( Figure 5). We consider an anomalous region in the depth range of 200-325 km and change systematically its nature to (i) thermal, (ii) chemical composition, and (iii) seismic velocities (V P and V S ).

Thermal Anomalies
Mantle flow can change the temperature in the sublithospheric mantle producing colder and hotter domains, which can be represented as thermal anomalies in LitMod2D_2.0. To account for these thermal perturbations, we modify the constant adiabatic thermal gradient such that the imposed temperature perturbation (ΔT) is added to the steady state temperature distribution. Then, LitMod2D_2.0 recalculates the relevant physical properties (density, seismic velocities, phase changes, and thermal conductivity) at T + ΔT and P conditions, where T and P, are pressure and temperature at a given depth below the LAB.
To illustrate the effect of hot and cold thermal anomalies on the sublithospheric mantle with a DMM composition, we consider a thermal anomaly ranging from −400°C to +400°C in steps of 100°C ( Figure 6). As expected, cold thermal anomalies increase density and P and S wave velocities, while hot anomalies have the opposite effect. The amplitude of the resulting anomalies varies such that density has the least relative change, whereas S waves have the highest change due to the high sensitivity of S wave velocities to temperature. To first order, the absolute density change depends linearly on the sign of the thermal anomaly; thus, similar perturbations (in magnitude) can be seen at both sides of the 0.0 % anomaly in Figure 6b. In contrast, seismic velocities show higher relative changes for positive temperature perturbations than for the equivalent negative ones (Figures 6c and 6d). This is because of the enhanced anelastic attenuation at higher temperatures. Moreover, the amplitude of seismic velocity anomalies

10.1029/2019GC008777
Geochemistry, Geophysics, Geosystems decreases with depth due to the decreasing attenuation of seismic waves (low Q factor, Figures 2) at lithosphere-sublithosphere transition depths. This nonlinear dependence of seismic velocities on temperature has strong implications for qualitative interpretation of seismic tomography models (Cammarano et al., 2003). A second-order variation in density and seismic velocities occurs around 310 km depth (Figure 6, inset) coinciding with the orthopyroxene-clinopyroxene ( Figure S1 in the supporting information) transition and could be the reason for the X-discontinuity imaged in seismic data (Revenaugh & Jordan, 2008).

Chemical Compositional Anomalies
Chemical heterogeneity in the upper mantle can result from a number of processes (e.g., fluid/melt percolation, melt extraction, and replacement of mantle parcels by advection). In turn, delamination, slab break-off, and slab tear processes can induce the sink of pieces of cold lithosphere, with different chemical composition, into the sublithospheric mantle. For the case of chemical compositional anomalies, LitMod2D_2.0 recalculates the relevant physical parameters at the corresponding P-T conditions according to the prescribed composition in a given region beneath the LAB. For the case of combined thermal and compositional anomalies, the relevant parameters are recalculated at T + ΔT and P conditions from the prescribed chemical composition. Figure 7 shows a synthetic example considering compositional and a combination of thermal and compositional anomalies, where the separate effect of temperature is also shown. In the case of a compositional anomaly, we have considered an Average Garnet Tecton chemical composition, which is depleted with respect to DMM (Table 1). Results show that the compositional anomaly decreases density by~0.4% (Figure 7a), which is almost equal to the average decrease in P wave velocities (~0.4%, Figure 7b) and higher than the decrease in S wave velocities (~0.2%, Figure 7c). A combination of compositional and thermal anomalies (of +50°C) results in a maximum decrease in density and seismic velocities (Figure 7) because increasing temperature also tends to decrease density and seismic velocities. Indeed, we can see that P wave velocities are more sensitive to composition than S wave velocities. By studying the effect of melt removal from a fertile composition (range of Mg#) on S wave velocities, Priestley and McKenzie (2006) have also reported low sensitivity of S waves to composition. Note that depending on the chemical composition of the anomaly (degree of depletion) and the temperature perturbation, the sublithospheric anomalies can have competing effects on density and seismic velocities.

Seismic Velocity Anomalies
LitMod2D_2.0 allows for implementing seismic velocity anomalies by providing the absolute velocity values or the velocity variation, in percentage, relative to a reference model. The anomalies are incorporated in two ways: (1) giving the average absolute or relative value of the anomaly over a predefined region beneath the LAB and (2) varying the absolute or relative magnitude of the anomaly along up to five depth levels beneath the LAB. In the first case, a chemical compositional anomaly can also be assigned to the anomalous region. In the second case, a separate input file including horizontal distances and depth levels together with the anomalous values must be specified (Text S2 and Figure S3). LitMod2D_2.0 converts prescribed seismic velocity anomalies (relative or absolute) to temperature using the temperature derivatives of seismic velocities for the assigned chemical composition in the anomalous region. Once seismic velocity anomalies are converted to temperature, they are treated as thermal anomalies (Section 3.1) and in case of assigned chemical composition as a combination of thermal and compositional anomalies (Section 3.2).
In this section, we focus on anomalous seismic velocities (V P and V S ) and their translation into temperature and densities. We consider relative ±1.5% P and S wave anomalies below the LAB and discuss equivalent thermal and density recovery (Figure 8). Results show that for a given magnitude of seismic anomaly, 1.5% in the presented example, P wave translates into a higher temperature anomaly of~200°C than that from S wave (~100°C) (Figure 8a) which is consistence with the higher sensitivity of S waves to temperature. In other words, a given temperature change requires a higher percentage variation in the S wave velocities than in the P wave velocities as shown in Figure 6. Subsequently, density change ( Figure 8b) is higher (~0.6%) in case of P wave velocity input than S wave input (~0.4%). Note that density changes are not symmetric for positive and negative velocity anomalies (Figure 8b) as is the case for temperature change (Figure 8a). Changes in density depend on the P-T conditions, which control stable phase and mineral assemblages.

Open Source GUI and Input/Output Data
To make LitMod2D_2.0 more accessible to users, we have developed an open source GUI in Python, which is not platform depending and can be easily updated/modified by the user. A detailed manual is provided with the new LitMod2D_2.0 version distributed through online GitHub repository (https://github.com/ajay6763/ LitMod2D_2.0_package_dist_users). Here, we briefly discuss and illustrate the main features. The new GUI

10.1029/2019GC008777
Geochemistry, Geophysics, Geosystems allows the organization of different projects in separate folders containing the relevant surface geophysical observables (e.g., elevation, Bouguer and free-air gravity anomaly, geoid height, and surface heat flow) and the material files with petro-physical properties for each project. Mantle material files are lookup tables of all relevant physical properties (density, thermal conductivities, and seismic velocities) as functions of pressure, temperature, and chemical composition. These are computed by GENERATOR (Section 2) using components of the software Perple_X (Connolly, 2005(Connolly, , 2009). An option to plot previously digitized data on the background (e.g., Moho depths, LAB depths, interpretation from active seismic lines) is also provided.
The GUI main window allows defining the geometry of the model made up of different bodies, each one with its physical properties. Different buttons are provided in the top and bottom of the main window to add and delete bodies, edit properties and shape of bodies, and run the model (Figure 9a). After the run finishes, a new window displays the geometry of the model and the fit to surface observables (Figure 9b, left panel), and the calculated temperature, density, and seismic velocities distributions (Figure 9b, right panel). Results are plotted in an interactive Matplotlib environment and can be modified for publication purposes. All the surface observables (elevation, Bouguer and free-air gravity anomalies, geoid height, and surface heat flow) are saved and can be used for customized visualization and further processing of the model output. Likewise, a master output file containing the Cartesian coordinates of the grid nodes and the corresponding temperature, pressure, seismic velocities, density, and material file code is saved to be used in the post-processing toolbox according to the particular needs of the user (Section 6).

Application
To illustrate the applicability of the new LitMod2D_2.0 package, we performed a synthetic model representing a simplified active margin with a subducting slab and the associated mantle wedge. The objective of this section is to show the functionalities and practicalities of LitMod2D_2.0, as well as the possibilities of the post-processing toolbox rather than studying specific aspects of a real subduction setting.

Input Data, Model Geometry, and Physical Properties
The synthetic model is 1,000 km long and comprises three regions: the oceanic domain, the active subduction zone, and the continental domain (Figure 9a). The oceanic domain is characterized by a bathymetry of 5.5 km, and a 6-km thick crust without sediments, and a 100-km thick lithospheric mantle with a composition calculated by 6% fractional melting considering DMM as source. Active subduction zones are characterized by inverted isotherms such that the 1,300°C can be found at three different depths for the same horizontal position, resulting in a triple LAB. This boundary condition cannot be solved by the 2D heat transport equation under steady-state incorporated in the LitMod2D_2.0 code, and hence, the subducted oceanic lithosphere is modeled as a sublithospheric thermo-compositional anomaly, with ΔT = −300°C and with an oceanic lithosphere chemical composition. Similarly, the mantle wedge is modeled as a sublithospheric thermal anomaly with ΔT = +100°C. This simple treatment of the active subduction zone allows calculating the average physical properties of the anomalous bodies and their effects on elevation, seismic velocities, and potential fields. The stable continental domain consists of a 35-km thick crust and a 115-km thick lithospheric mantle with Average Garnet Tecton composition. Near the subduction zone the crust thickens up to 40 km, whereas the lithospheric mantle thins to minimum values of 78 km. For simplicity, we have considered a single crustal layer, but we can add as many layers and bodies as needed through the GUI. Thermo-physical parameters in the crust and chemical composition in the lithospheric mantle are incorporated via GUI and correspond to those listed in Tables 2 and 1, respectively. Figure 9b shows a screenshot of the main results window, which includes the calculated and measured surface observables together with the model geometry and the calculated temperature, density and P and S wave velocities distribution. The anomalous sublithospheric bodies have a clear signature not only on the temperature distribution, as they have been defined as thermal anomalies, but also on the calculated density and seismic velocities. As observed, along the subducting slab, the combined thermal and compositional anomalies increase density and P and S wave seismic velocities, whereas along the mantle wedge, the temperature increase of 100°C results in a decrease of the three observables.
As the modeled profile is a synthetic model, we have considered that the observed elevation, gravity, geoid, and heat Geochemistry, Geophysics, Geosystems anomalies, we made an additional run without sublithospheric anomalies, keeping all other parameters unmodified (Figure 9b).
LitMod2D_2.0 calculates the absolute elevation from the final density distribution, under the assumption of local isostatic equilibrium with a compensation depth at the base of the model (400 km). Lateral density variations beneath the LAB will tend to generate dynamic Stokes flow in the underlying less viscous mantle, which can be enhanced by sublithospheric mantle anomalies. Vertical stresses associated with this dynamic flow can be transmitted to the surface topography depending on viscosity distribution with depth. By default, LitMod2D_2.0 estimates the influence of sublithospheric mantle anomalies on the calculated elevation by considering two end member conditions representing the upper and lower bounds of the dynamic topography contribution Tunini et al., 2016;Jiménez-Munt et al., 2019): (1) coupled elevation, when the vertical stresses induced by the buoyancy of sublithospheric mantle anomalies are completely transmitted to the surface and (2) uncoupled elevation, when sublithospheric mantle anomalies have no effect on topography. A more complete treatment of the problem would require using the density and viscosity distribution, calculated as a function of pressure, temperature, and material type from LitMod2D_2.0, and solve the Stokes flow equation for the whole modeling domain. The resulting normal stresses at the surface or at the base of the lithosphere allow calculating a more realistic dynamic topography contribution (e.g., Afonso, Rawlinson, et al., 2016). In the synthetic model presented here, the net effect of the sublithospheric anomalies on elevation is of ≤1,000 m ( Figure 9b). In the case of coupled anomalies, the cold slab near the subduction zone pushes down the elevation because of its higher density. The relatively lower density of the mantle wedge does not suffice to cancel the negative slab buoyancy. Density perturbations from the cold subducted slab and the hot corner flow primarily affect the geoid height and free-air gravity anomaly across the model, whereas Bouguer anomaly is slightly affected. This is because geoid is sensitive to the density moment which in the modeled synthetic profile is mainly related to the topography and the sublithospheric anomalies. Free air is very sensitive to topography and in a less extent to sublithospheric anomalies due to the inverse square law dependence, whereas Bouguer anomaly is also depending on the inverse square law but corrected for topography effects and then is essentially sensitive to lateral crustal density variations.

Post-processing Toolbox
The main idea of the post-processing toolbox is that from the obtained parameters characterizing the physical state of the crust and upper mantle (pressure, temperature, density, P and S wave seismic velocities, and material type) at each node, the user can further process these data to do a variety of additional calculations, such as strength envelopes, flexural effects on elevation, Stokes flow, distribution of stable mineral assemblages, and synthetic passive seismological data. The post-processing toolbox is a set of independent scripts/codes linking the LitMod2D_2.0 output with other softwares (Figure 1). The user can use, adapt, or create new scripts or codes according to her/his own interest and needs.

Post-processing Example
The following example illustrates the potential of the post-processing toolbox focused, in the presented example, on the calculation of synthetic seismic tomography and synthetic passive seismological data.
We provide shell scripts to compute synthetic tomography from the LitMod2D_2.0 seismic velocities, using the reference model defined in Section 2.3. For the presented synthetic example, maximum positive anomalies are observed along the oceanic lithosphere, whereas the thinned continental lithosphere above the subducted slab results in negative velocity anomaly (Figure 10a). The thick continental lithosphere results in a positive velocity anomaly but with a smaller magnitude relative to the oceanic lithosphere. In the sublithospheric mantle, the subducted slab and corner flow regions result in positive and negative velocity anomalies, respectively, and rest of the sublithospheric mantle shows no deviations from the reference model.
Joint inversion of receiver functions and surface-wave dispersion has been widely used to image major velocities discontinuities and absolute seismic velocity distributions with depth (e.g., Julià et al., 2000;Langston, 1979;Vinnik, 1977). P wave receiver functions are sensitive to the S wave velocity discontinuities where P-to-S converted waves sample the subsurface structure, whereas surface waves are sensitive to the average S wave velocities distribution.

Geochemistry, Geophysics, Geosystems
LitMod2D_2.0 gives depth distribution of absolute P and S wave velocities from a structural (e.g., Moho and LAB geometries) and chemical composition model. Hence, comparing observed surface-wave dispersions and P wave receiver functions with those inferred from LitMod2_2.0 can further constrain the obtained models. This is done by coupling the calculated P and S wave seismic velocities from LitMod2D_2.0 with the "Computer Programs in Seismology" tool by (Herrmann, 2013) through Shell and Python scripts to calculate synthetic Rayleigh surface-wave dispersion curves and P wave receiver functions. Within the crust, P wave velocities are calculated using empirical V P -density relationships from Brocher (2005), whereas S wave velocities are calculated assuming a constant V P /V S = 1.73 ratio (Figure 10b). Indeed, an option including crustal material files is kept for future or if user have a thermodynamic database for the crustal chemical
The depth distributions of the resulting S wave velocity at three locations along the synthetic profile are shown in Figure 10b corresponding to (a) ocean (300 km), (b) subduction zone (500 km), and (c) continent (950 km). Synthetic P wave receiver functions (Figure 10c) show a clear positive converted phase (Ps) at Moho discontinuity, which arrives earlier for oceanic domain because of the lower crustal thickness, and is delayed for thick continental crust. Each of the converted phase at Moho (Ps) has a positive phase (Ppps) and a negative phase (Ppss) multiples which are helpful in determining whether the Ps phase corresponds to a velocity discontinuity. Rayleigh group and phase velocities (Figure 10d) are faster for the oceanic domain than for the continental domain because of the high S wave velocities of the oceanic lithosphere (Figures 9b and 10b). At short periods, velocities are lower than at higher periods because shallow structures (slow velocities) are more sensitive to surface waves at short time periods. Group velocities show an absolute maximum, whereas phase velocities show relative maxima increasing for higher periods. This behavior is because of the differential depth sensitivity of group and phase velocities. For the oceanic domain, the maxima of group velocities is~25 s, whereas for the subduction zone and the continental domains maxima is reached at~75 s. Phase velocities show that local maxima are shifted toward lower periods relative to the maximum of group velocities. The effect of lithosphere thickness is observed in the continental domain showing higher velocities (both group and phase) in regions with thick lithosphere than in the thin lithosphere near the subduction zone. The effect of the cold subducted slab is clearly visible in both group and phase velocities (black dashed line) where maxima for group velocities (local maxima in case of phase velocities) has been shifted toward high periods and velocities are increased in comparison to without slab dispersion curves.

Summary
We have presented an improved version of the LitMod2D package to study the thermal, compositional, density, and seismological structure of the crust and upper mantle (~400 km depth) along cross sections by integrating geophysical and geochemical data set in a self-consistent thermodynamic manner. The main changes and improvements included in the new LitMod2D_2.0 software package are as follows: 1. We use DMM (depleted MORB mantle) instead of PUM as the chemical composition in the sublithospheric mantle for geochemical consistency. Nevertheless, the relative differences with respect to PUM in the calculated density and P and S wave seismic velocities are less than 0.2 %. 2. Anelasticity is now calculated in the GENERATOR thermodynamic module using updated parameters in the power-law attenuation model allowing for the characterization of sublithospheric mantle anomalies. Both oscillation period and grain size remain as input parameters. 3. Sublithospheric mantle anomalies may have a thermal, seismic, or compositional origin or a combination of them, allowing for estimating the response of geodynamic processes such as mantle upwelling, subduction, delamination, and metasomatism on the surface observables (gravity, geoid, and elevation). 4. To compute synthetic V P and V S tomography models, we use an upper mantle model consisting of Average Garnet Tecton composition in the lithospheric mantle, with Z Moho = 35 km and Z LAB = 120 km, and a DMM composition in the sublithospheric mantle down to 400 km depth. This model reproduces the P wave seismic velocities of the global average reference model ak135 in the depth range of 35-250 km and avoids the unsolved mismatch between the thermodynamically calculated seismic velocities and those from ak135 model below 250 km depth. 5. We have developed a new open-source GUI under Python programming language simplifying the input/output data files to obtain across platforms versatility, to gain ease of use, and to allow future user developments. 6. We incorporated a post-processing toolbox allowing the user to apply predefined/customized codes and scripts linking with the LitMod2D_2.0 output file to calculate additional observables according to specific needs (synthetic tomography, regional isostasy, stable mineral assemblages, surface-wave dispersion curves and receiver functions, etc.). The user can customize the provided codes/scripts or develop new ones to be shared with the scientific community.