Diagnosing open-system magmatic processes using the Magma Chamber Simulator (MCS): part II—trace elements and isotopes

The Magma Chamber Simulator (MCS) is a thermodynamic model that computes the phase, thermal, and compositional evolution of a multiphase–multicomponent system of a Fractionally Crystallizing resident body of magma (i.e., melt ± solids ± fluid), linked wallrock that may either be assimilated as Anatectic melts or wholesale as Stoped blocks, and multiple Recharge reservoirs (RnASnFC system, where n is the number of user-selected recharge events). MCS calculations occur in two stages; the first utilizes mass and energy balance to produce thermodynamically constrained major element and phase equilibria information for an RnASnFC system; this tool is informally called MCS-PhaseEQ, and is described in a companion paper (Bohrson et al. 2020). The second stage of modeling, called MCS-Traces, calculates the RASFC evolution of up to 48 trace elements and seven radiogenic and one stable isotopic system (Sr, Nd, Hf, 3xPb, Os, and O) for the resident melt. In addition, trace element concentrations are calculated for bulk residual wallrock and each solid (± fluid) phase in the cumulate reservoir and residual wallrock. Input consists of (1) initial trace element concentrations and isotope ratios for the parental melt, wallrock, and recharge magmas/stoped wallrock blocks and (2) solid-melt and solid–fluid partition coefficients (optional temperature-dependence) for stable phases in the resident magma and residual wallrock. Output can be easily read and processed from tabulated worksheets. We provide trace element and isotopic results for the same example cases (FC, R2FC, AFC, S2FC, and R2AFC) presented in the companion paper. These simulations show that recharge processes can be difficult to recognize based on trace element data alone unless there is an independent reference frame of successive recharge events or if serial recharge magmas are sufficiently distinct in composition relative to the parental magma or magmas on the fractionation trend. In contrast, assimilation of wallrock is likely to have a notable effect on incompatible trace element and isotopic compositions of the contaminated resident melt. The magnitude of these effects depends on several factors incorporated into both stages of MCS calculations (e.g., phase equilibria, trace element partitioning, style of assimilation, and geochemistry of the starting materials). Significantly, the effects of assimilation can be counterintuitive and very different from simple scenarios (e.g., bulk mixing of magma and wallrock) that do not take account phase equilibria. Considerable caution should be practiced in ruling out potential assimilation scenarios in natural systems based upon simple geochemical “rules of thumb”. The lack of simplistic responses to open-system processes underscores the need for thermodynamical RASFC models that take into account mass and energy conservation. MCS-Traces provides an unprecedented and detailed framework for utilizing thermodynamic constraints and element partitioning to document trace element and isotopic evolution of igneous systems. Continued development of the Magma Chamber Simulator will focus on easier accessibility and additional capabilities that will allow the tool to better reproduce the documented natural complexities of open-system magmatic processes.


Introduction
Phase equilibria in magmatic systems at given P-T-conditions are controlled by the major element composition and thermodynamic properties of the magma. Major (± minor) elements are usually reported as oxides in geochemical datasets and are primary stoichiometric constituents of mineral phases that are stable in the magma. Modeling of phase equilibria and major element evolution of an igneous system requires knowledge of the thermodynamic properties of the liquid + solid + fluid phases present in the system. These include the standard state properties of all end-member components in all phases as well as the activity-composition relations of all phases.
The Magma Chamber Simulator (MCS; Bohrson et al. 2014) is a computational tool that uses the family of MELTS engines (Ghiorso and Sack 1995;Ghiorso et al. 2002;Gualda et al. 2012;Ghiorso and Gualda 2015) to numerically quantify these parameters in an isobaric open-system resident magma body subsystem (M: melt ± solids ± fluid) that evolves by fractional crystallization (FC) and can simultaneously experience magma recharge/mixing (R) and assimilation of either anatectic melts (A) and/or stoped blocks (S) of wallrock (WR). The utilization and application of MCS-PhaseEQ and its bearing on the evolution of phase equilibria and major elements have been described in a companion paper (Bohrson et al. 2020).
Trace elements are found in dilute quantities (usually ≤ 0.1 wt%) and are not primary stoichiometric constituents of common phases in magmatic systems. They generally do not influence phase equilibria, but their relative concentrations and isotopic compositions can nevertheless record important information about the sources and evolution of magmatic systems. The activities of trace elements vary in direct relation to their concentrations, and their behavior is controlled primarily by their relative compatibility among melt, solids (i.e., minerals) and fluid. This behavior is defined by a partition coefficient ( K ), which in igneous petrology is a ratio of the concentration of an element in a mineral versus the concentration of the element in the melt ( K sm ) or fluid ( K sf ) in equilibrium with the mineral. Because information on the composition and relative amounts of stable phases in a magmatic system is required to constrain bulk partition coefficients ( D sm and D sf ), detailed modeling of trace elements is not possible without a record of phase equilibria of the system.
In this study, we describe the MCS-Traces tool that uses primary MCS-PhaseEQ output as input and adds equations for modeling up to 48 trace elements in an R n AS n FC system. The program input consists of MCS-PhaseEQ output, initial trace element concentrations for the M parental melt, R n + S n (n tot ≤ 30), and WR, and phase-specific K sm (± K s∕H 2 O ± K s∕CO 2 for possible fluid phases) values for M and WR. An option to take into account the T-dependence of K values is available. In addition to trace elements, MCS-Traces computes radiogenic isotope ( 87 Sr/ 86 Sr, 143 Nd/ 144 Nd, 176 Hf/ 177 Hf,206 Pb/ 204 Pb,207 Pb/ 204 Pb,208 Pb/ 204 Pb,and 187 Os/ 188 Os as defaults) and stable O isotope (δ 18 O) compositional evolution of the open-system resident melt based on user-input isotope ratios for the parental melt, WR, and each R magma. A description of the tool is provided along with results and comparisons of case studies that correspond to the ones introduced in the companion paper. An example of MCS-Traces applied to a natural magmatic system (continental flood basalts) is given in Heinonen et al. (2019).

How does MCS-Traces work?
MCS-Traces utilizes the primary output ("RunSummary") of the MCS-PhaseEQ modeling tool. This output is read as the manual for what is going on in the R n AS n FC system (in terms of phase equilibria and mass balance) and MCS-Traces executes the trace element and isotope calculations based on it and user-input initial concentrations, isotope ratios, and K values. Primary output of fluid composition is used by the MCS-Traces tool to calculate K sf values and temperature information is required for runs with T-dependent K values. Because of these utilizations, MCS-Traces cannot be used without output from a successfully executed MCS-PhaseEQ model. The MCS-derived trace element solution has the advantage of being self-consistent and more complete in the sense that a proper reckoning is made between trace elements and isotopes and major element phase compositions, modal abundances, pressure, and temperature.
The tool itself is composed of two main parts: (1) the user interface built in Visual Basic environment, and (2) the trace element and isotope engine that is a separate Excel worksheet, in which all the computations are performed. The current versions are available for both Mac and PC; the reader is referred to the website for the most up-to-date information on the different versions. In short, the user interface feeds input into the engine and, after the calculations have been performed, collects the output in an easily accessible and digestible form that is then amalgamated as a set of additional worksheets with the original MCS-PhaseEQ output. Consequently, the entire petrological-geochemical solution to a complex RASFC scenario is given in one workbook with multiple tabs. More detailed descriptions are given in the following sections and on the MCS website (https ://mcs.geol.ucsb.edu/).

Trace element and isotope equations
The basic framework of modeling trace elements and isotopes in MCS-Traces is described in Bohrson et al. (2014); we will complement and update that description here. For M, each temperature step is considered to be an equilibrium crystallization step. At the end of the step and before the model proceeds, however, the newly formed solid and fluid phases are fractionated from M in separate cumulate reservoirs. The model thus never represents a pure equilibrium or fractional crystallization case and thus better mimics natural systems-the larger the user-defined temperature step is, the larger is the role for equilibrium crystallization. Wallrock is always in internal equilibrium. The trace element equations for equilibrium crystallization in a melt ± solid ± fluid system developed by Spera et al. (2007) form the backbone of the trace element engine: In Eqs. 1, 2, 3, C m , C s , and C f stand for trace element concentration in the melt, solid, and fluid, respectively, whereas C 0 is the initial trace element concentration in the system in internal equilibrium. Because of the structure of the model code, C 0 in Eqs. 1, 2, 3 represents the composition of M melt initially or at the end of the previous fractionation step. In Eqs. 1, 2, 3 for WR, which is always in internal equilibrium, C 0 represents the bulk composition of WR initially or after the previous anatectic melt removal event. D sm and D sf are solid/melt and solid/fluid bulk partition coefficients for a given element, and f m and f f mass fractions of melt and fluid, respectively. D values for M ( D M sm and D M sf ) and WR ( D WR sm and D WR sf ) are calculated on the basis input element-specific K values and the mass fractions of stable equilibrium phases ( x ) for each step, as shown below: Respective K sf values have to be input only if a separate fluid phase consisting of either H 2 O or CO 2 or both (possible in MCS-PhaseEQ runs using rhyolite-MELTS 1.1.0 or 1.2.0) saturates in M and/or WR. When the fluid phase contains both components, the effective K sf is computed according to: in which w H 2 O is the mass fraction of H 2 O in the fluid phase and K s∕H 2 O and K s∕CO 2 are the solid-H 2 O and solid-CO 2 partition coefficients, respectively. (1) Note that although phase proportions are calculated in MCS-PhaseEQ for recharge magma and stoped WR blocks, they are homogenized before mixing with M, and thus, the partitioning of elements into solid, melt and fluid is not relevant to their mass balance. Therefore, constraints on K sm ± K sf are only needed for M and, in the case of AFC, for WR. Note that any possible negative solid and fluid masses produced by MCS-PhaseEQ (these are related to how MELTS algorithm performs the calculations; see MCS website) have been approximated as zero for the trace element calculations. This ensures that these very small negative masses of mineral phases do not cause negative trace element concentrations for those elements which have high K sm or K sf values. The effect of this approximation on the trace element model results has been determined to be negligible (< 0.1% difference in regards to total masses, i.e., less than the analytical error of any major and trace element geochemical data).
Also, T-dependence of K sm and K sf can be invoked in MCS-Traces. The engine calculates the K at a given temperature [°C] logarithmically based on A and B values, which, in turn, are calculated based on two user-input K values at two different temperatures [°C]: The temperature-dependence of K is identical in form to that of any reaction with the isobaric heat capacity change of the reaction set to zero. Note that for many elements and solid phases, the effects of composition (of melt and the solid phases themselves) and temperature (and other physical and chemical factors) on K values are difficult to distinguish from each other (e.g., Blundy and Wood 2003). These effects, however, are also effectively intertwined in magmatic systems, and, for example, T-dependency can thus be used to approximate the compositional dependency with some caution (e.g., Blundy and Wood 1991).
When a mass of WR anatectic melt above the user-input percolation threshold value (FmZero) or R or S is added to M melt, simple mixing equations are utilized to calculate the resulting hybridized composition of M melt: In the case of assimilation of WR partial melt, the assimilated mass of the trace element is then removed from the WR reservoir.
In the case of isotopes, isotopic equilibrium within the subsystems (M, WR, R) is assumed and high-T isotope fractionation during melting and crystallization is not taken into account. Minor high-T fractionation of relatively light O isotopes during melting and crystallization is a well-known process, but the geochemical effect of this fractionation is generally quite small (usually < 1% in a magmatic series) compared to that of adding isotopically distinct components into resident melt (see, e.g., Eiler 2001). Therefore, and because additional input of such fractionation factors would be needed for each stable phase, high-T fractionation of O isotopes is not considered in MCS-Traces. Simple mixing between the subsystems is assumed in the case of all isotope compositions ( ): For O isotopes, O mass fractions are calculated for the M melt, R, and WR melt (or S) on the basis of MCS major element output at each step.

General overview of using MCS-Traces
For detailed instructions on the operation and input and output of both MCS tools, the reader is referred to the MCS website (http://mcs.geol.ucsb.edu); a generalized overview of the MCS-Traces is presented here.
MCS-Traces should be run with the same platform that was used to run MCS-PhaseEQ. MCS-Traces should never be updated or configured separately from MCS-PhaseEQ and, in the case of possible updates, we always recommend downloading the most up-to-date versions of both programs from the website. An internal compatibility check performed by the software will prevent the use of outdated or otherwise incompatible MCS-PhaseEQ output in the MCS-Traces tool. If one has used MCS-PhaseEQ to produce output required for MCS-Traces, the installation of the whole MCS package has been successful.
The general layout of MCS-Traces interface is presented in Fig. 1. The command buttons are set so that the default work flow proceeds primarily downwards and secondarily from left to right. Input consists of initial concentrations and phase-specific K sm values (+ K s∕H 2 O and/or K s∕CO 2 if fluid is present in the simulation) for up to 48 trace elements and initial isotopic compositions for up to eight isotopic systems (Sr,Nd,Hf,3xPb,Os,and O). Input is given separately for the different subsystems (M, WR, and R). Note that Ti, Cr, Ni, Mn, and Co are forced disabled as trace elements if these (10) were modeled as respective oxides in MCS-PhaseEQ. The concentration input is not bound to a unit so the user can input the concentrations in ppm, ppb, or any other unit as long as the unit is uniform across the different subsystems. Input can be saved and previous input loaded in the form of PAR files. Some recommended sources for K sm values for igneous systems are Rollinson (1993), the Geochemical Earth Reference Model (GERM) Partition Coefficient Database (https ://earth ref.org/KDD/), and EarthChem traceDs (http://earth chem.org/trace ds), which is in preparation, but already includes up-to-date compilations of K sm values for some common igneous phases. A compilation of K sf values is available in Spera et al. (2007), and the reader is also directed to the original publications that support this compilation. PAR files which store input for all the MCS-Traces simulations of this study are included in Online Resource 1. Although the associated partition coefficient input data can be loaded in MCS-Traces to be utilized for any simulation, we recommend users to search for partition coefficients that are the most suitable for the geological environment to be modeled.
Once the computation is complete, the MCS-Traces output can be accessed by opening the original MCS-PhaseEQ output file (see Online Resource 2). Two to three new tabs in the MCS-PhaseEQ output worksheet have been produced. Information regarding the trace element and isotope budget in the magma chamber and wallrock (if present in the simulation) are compiled in the "Traces_X_Magma" and "Traces_X_Wallrock" tabs, respectively. Both are appended with echoed input and relevant data from the primary MCS-PhaseEQ output. Finally, tab "Traces_X_OUTPUT" tabulates the full output from the MCS-Traces computation. Note that MCS-Traces can be run multiple times against the same MCS-PhaseEQ output file, the double/triplet tab sets are then distinguished by numbers (X = 1, 2, 3…) as they are introduced.

Comparison of closed (fractional crystallization) and open-system magma evolution illustrated by MCS
We present and compare trace element (Ni, Sr, Zr, Ba, La, Nd, and Yb) and isotopic ( 87 Sr/ 86 Sr, 144 Nd/ 143 Nd, and δ 18 O) results of the same five MCS simulations introduced in the companion paper (Bohrson et al. 2020). These models are: (1) fractional crystallization (FC), (2) fractional crystallization + recharge (n = 2) (R 2 FC), (3) fractional crystallization + assimilation of wallrock anatectic melts (AFC), (4) fractional crystallization + assimilation of stoped wallrock blocks (n = 2) (S 2 FC), and (5) fractional crystallization + recharge (n = 2) + assimilation of wallrock anatectic melts (R 2 AFC) ( Table 1). The input and results of MCS-PhaseEQ modeling are only briefly outlined here to place the trace element and isotope results in context. The outcomes of the simulations are illustrated as block models built with the help of MCS Visualizer (see MCS website) in Fig. 2. The reader should refer to the companion paper for details of initial conditions and thermal, mass, phase equilibria, and major element output. Note that the presented models have been built dominantly for the purpose of displaying the features of MCS, and are not directly applicable to a particular natural case. They broadly represent a setting of anorogenic magmatism, where magmas from a hot asthenospheric source intrude shallow levels of cratonic continental crust. We encourage MCS users to carefully select the parameters for the specific system to be simulated (see Bohrson et al. 2014).
The five example simulations use a depleted mantlederived tholeiite from the ~ 180 Ma Karoo large igneous province (sample P27-AVL; Luttinen and Furnes 2000) as the initial melt and recharge magma major and trace element composition with ~ 2 wt% of H 2 O added (Table 2). Initial redox conditions were set at quartz-fayalite-magnetite (QFM) buffer, but the run itself was closed for oxygen, i.e., oxygen fugacity was determined by the relevant internal chemical equilibria of the system. The Sr and Nd isotope compositions were set to those of mean MORB (Gale et al. Fig. 1 The user interface of the MCS-Traces tool. The number of buttons on rows 3-5 vary according to the uploaded MCS-PhaseEQ output. In this run, which is the R 2 AFC case of this study, both M and WR contain a fluid phase composed of a mixture of H 2 O and CO 2 , and two R events took place Table 1 MCS-Traces models discussed in this study The K values for the different MCS-Traces models are presented in Table 3 The underlying MCS-PhaseEQ models are presented and discussed in detail in Bohrson et al. (2020) (Eiler et al. 2000). The major element, trace element, and isotope compositions of the wallrock and the stoped blocks represent the granodioritic average upper continental crust of Rudnick and Gao (2003) with ~ 2 wt% of H 2 O and 1 wt% of CO 2 added and initial Fe oxidation state set at QFM ( Table 2). The Sr and Nd isotope composition has been constrained on the basis of global river waters (Goldstein and Jacobsen 1988) and the O isotope composition is within the range of felsic crustal rock types (Eiler 2001). The mass of WR was twice that of the initial M, and the percolation threshold (FmZero) of the WR was set to a melt mass fraction of 0.1. In the models involving the assimilation of wallrock anatectic melts (AFC and R 2 AFC), the initial temperature of the wallrock was set at 700 °C (close to its solidus) to simulate preheating by hypothetical earlier magma pulses and maximize the assimilation potential for comparison with the other cases.

MCS-Traces
The degree of contamination of the resident magma can be presented in two sensible ways in the models: (1) relative to the mass of the parental melt or (2) relative to the mass of the resident magma chamber (the whole M subsystem). The former may be relevant for cases without a recharge, but here we use the latter approach that takes account recharge pulses Fig. 2 Results of MCS-PhaseEQ and MCS-Traces simulations for the default five cases (FC, AFC, R 2 FC, S 2 FC, R 2 AFC) shown in annotated snapshots that depict the situation after the final magma crystallization step (AFC and R 2 AFC include one additional step of wallrock equilibration before the simulation ends). Completions of R and S events and beginning of A are indicated in the cumulate pile where applicable. Note that the phase proportions are based on mass fractions not volume fractions and that the wt% of the subsystems are relative to the whole magma-wallrock system; M melt, M fluid, and M cumulate (M cumul.) comprise the total magma chamber mass. For the subsystems, SiO 2 is given in wt% and trace elements Ni, Sr, and Ba in ppm. See Bohrson et al. (2020) for more information. Mineral abbreviations: ol olivine, opx orthopyroxene, cpx clinopyroxene (FC and S 2 FC include two separately output cpx solid-solution phases), plag plagioclase, qtz quartz, spl spinel, rhm rhombohedral oxide and is in line with the companion paper. Note that the models that do not involve assimilation of WR anatectic melts (FC, R 2 FC, and S 2 FC) still included WR in the modeling for the purpose of studying its thermal evolution due to heat released by the cooling and crystallizing magma. Therefore, initial composition and K WR sm values were part of the input in MCS-Traces modeling also in these three cases. Because the temperature of the wallrock does not reach its solidus, and hence no anatectic melt forms or is assimilated by the resident melt, the input values do not affect the outcome of these models, and can thus be ignored.
Models were run at 0.1 GPa with an M temperature decrement of 5 °C. The relatively low pressure was set to accommodate separation of fluid.
K sm values used in the models have been compiled from Mahood and Hildreth (1983); Bacon and Druitt (1988); Rollinson (1993); Bea et al. (1994); Ewart and Griffin (1994), a recent version of PELE software (see Boudreau 1999), and Kessel et al. (2005). and are listed in Table 3. We selected consistent sets of K sm values widely used and appropriate for the modeled compositions. Because MCS is thermodynamically controlled modeling tool, we utilized T-dependent K sm values for all phases that are stable in a wide range of temperatures and for which applicable element-specific data has been published (Table 3). Two additional trace element simulations were run for the AFC case to illustrate the effects of using hypothetical maximum and minimum values for K ol∕m (Ni), K cpx∕m (Ni), and K plag∕m (Sr). K sf values (for H 2 O and CO 2 ) are 1000 for all general runs to minimize incompatible element concentrations in the fluid phase. We also ran a FC trace element simulation that utilized more realistic D M sf values (Kessel et al. 2005), and illustrates the partitioning of mobile elements ( D M sf for Sr, Ba, La, and Nd < 1) into a separate fluid phase in the resident magma. The total number of MCS-Traces runs illustrated and discussed below is thus 8 ( Table 1).
The results of the simulations are discussed below and illustrated in Figs. 2, 3, 4, 5, 6, 7. Online Resources 1-3 include complete input and output for both MCS tools in case-specific files. We recommend viewing the outputs (Online Resource 2) simultaneously when reading the following sections. For a more detailed discussion on the phase equilibria and major element evolution in the MCS-PhaseEQ models, the reader is referred to the companion paper (Bohrson et al. 2020); concise summaries of those results are presented herein.

Case 1: fractional crystallization (FC)
In the primary MCS FC simulation, the melt composition evolves from basaltic (SiO 2 ≈ 49 wt% at liquidus temperature of ~ 1129 °C) to dacitic (SiO 2 ≈ 69 wt% at ~ 899 °C) by ~ 76 wt% of crystallization of the magma chamber ( Fig. 2a). Because of the high H 2 O content in the initial magma, a separate H 2 O fluid phase forms after ~ 46 wt% crystallization and accumulates up to ~ 1 wt% of the magma chamber in total by the end of the run. The fractionated phases and their approximate relative proportions in the bulk solid cumulate are in order of appearance: olivine (15 wt%), clinopyroxene (40 wt%), plagioclase (39 wt%), spinel (6 wt%), and rhombohedral oxide (0.4 wt%).
In terms of trace elements, Ni is compatible in the incremental solids throughout the default FC run (Fig. 3a), because all of the associated solid phases except plagioclase have K M sm (Ni) values of ≥ 1 (Table 3;   Major elements discussed and modeled in Bohrson et al. (2020). Sources for the trace element and isotope values are given in the text. PM is the recharge magma composition in the R 2 FC and R 2 AFC simulations and WR is the "recharge magma" (stoped WR block) composition in the S 2 FC simulation PM WR   Rollinson (1993;original sources: Philpotts and Schnetzler, 1970;Gill 1981;Nash and Crecraft 1985), c Ewart and Griffin (1994; average of minimum and maximum values for low-Si rhyolite; if not available, for leucitite, dacite, or peralkaline rhyolite), d Bacon and Druitt (1988); e Bea et al. (1994); f Mahood and Hildreth (1983) Note that K s∕H 2 O and K s∕CO 2 (and thus K sf ) are 1000 for all elements in all phases, except for the case with fluid mobility taken into account (Table 1), where bulk D sf values (uniform element-specific K sf value input for all phases) are as follows (after Kessel et al. 2005): a negligible effect on rare earth element (REE) ratios, and isotope compositions obviously remain unchanged (Fig. 4). The effect of adjusted D M sf values is negligible for elements with D M sf > 0.5 in the alternative FC trace element simulation (Fig. 5a, b). However, the effect becomes more important for more soluble elements, such as Sr and Ba with D M sf < 0.1, even when the total mass of fractionated fluid in the FC run does not exceed 1 wt % of the magma chamber mass (Fig. 5c, d). The Sr and Ba contents at the end of the adjusted D M sf run are about 70% and 30% lower, respectively, than at the end of the default run. Compared to Ba, the effect on Sr is proportionally larger because Sr has a higher D M sm after the onset of plagioclase fractionation for most of the run, whereas D M sm (Ba) stays below 0.3 throughout.  (e), and Yb (f) (ppm) diagrams for the resident melt in the example case studies (FC, R 2 FC, AFC, S 2 FC, and R 2 AFC). In (a), beginning of assimilation and R and S events are indicated. Percentages shown in SiO 2 vs. La for the cases involving assimilation (e) reflect the mass of assimilated material relative to the mass of the M subsystem. All of these are applicable for all diagrams at a given SiO 2 but are only shown in (a, e) to preserve clarity

Case 2: recharge-fractional crystallization (R 2 FC)
The R 2 FC model simulates replenishment of the magma chamber with two batches of mantle-derived magma that are equivalent to the parental melt composition ( Table 2). The two R events of 75 mass units (m.u.) are introduced at ~ 1049 °C (resident melt SiO 2 ≈ 52 wt%) and at ~ 998 °C (resident melt SiO 2 ≈ 58 wt%). Because of the large R events, the total magma chamber mass at the conclusion of the simulation is 250 m.u.. The second R magma is below liquidus temperature (at 1080 °C; T liq ≈ 1129 °C) and thus contains ~ 21 wt% of crystal cargo (but no fluid); this does not affect the trace element run, however, as bulk R magmas (i.e., all crystals + melt ± fluid) are always fully hybridized with the resident melt in MCS.
The R 2 FC model terminates after ~ 77 wt% of crystallization of the magma chamber with ~ 22 wt% (~ 55 m.u.) of dacitic (SiO 2 ≈ 69 wt% at 900 °C) melt left and ~ 1 wt% of fluid separated (Fig. 2b). The composition of the bulk solid cumulate is approximately 13 wt% olivine, 40 wt% clinopyroxene, 39 wt% plagioclase, 5 wt% spinel, 2 wt% orthopyroxene, and < 1 wt% rhombohedral oxide (listed in the order of appearance), which is very similar to that of FC with the exception of orthopyroxene being stable in the R 2 FC simulation.
At almost any given SiO 2 content, the concentrations and ratios of the modeled trace elements in the R 2 FC run are very similar to those produced by the FC runs (Figs. 3,  4a). The only notable differences are in regard to largely compatible elements Ni and Sr (Fig. 3a, b). The resident melt is enriched in Ni due to the added primitive R pulses (that have higher Ni contents than the resident melt at the instance of R), especially the second R event (R2) that triggers considerable precipitation of clinopyroxene (~ 15 m.u.), which has consistently lower K M sm (Ni) than olivine ( Table 3) that is the sole crystallizing phase before R2. In the case of Sr, R1 lowers the Sr content of the resident melt slightly, whereas R2 increases it; this is because the recharge events take place when resident melt Sr concentrations are higher (~ 234 ppm before R1) and lower (160 ppm before R2) than in the recharge magma (226 ppm; Table 2). In terms of fully incompatible elements (Fig. 3c-f), the recharge events do not cause notable deviation from the rather linear geochemical trends. The implications of these observations will be discussed later.

Case 3: assimilation of wallrock anatectic melts-fractional crystallization (AFC)
The AFC model is the first one involving the assimilation of wallrock that has a considerably higher initial temperature Fig. 4 Nd/Yb vs. La/Nd diagram (a) and Sr, Nd, and O isotope diagrams (b-d) for the resident melt in the example case studies (FC, R 2 FC, AFC, S 2 FC, and R 2 AFC). Note that FC and R 2 FC are not shown in the isotopic diagrams because these models do not result in isotopic changes in the resident melt. Percentages given for the model curves involving assimilation reflect the mass of assimilated material relative to the total mass of the M subsystem (magma chamber). Simple binary mixing models are also shown in the isotopic diagrams; percentages indicate the relative amount of WR in the mixture and are broadly comparable to the percentages in the MCS runs (there is no fractional crystallization in binary mixing, however) (700 °C) than in the FC and R 2 FC (and S 2 FC) cases (100 °C). The assimilation of WR partial melts begins when the temperature of the resident melt has dropped to ~ 1069 °C (~ 51 wt. SiO 2 ; ~ 26 wt% crystallization at that point). The heat released from the crystallizing M before this has heated the relatively large mass of wallrock (200 m.u.) from its initial temperature above the temperature (~ 740 °C) at which its percolation threshold (10 wt% of melting) is exceeded.
The AFC model terminates after ~ 43 wt% of crystallization at a total magma chamber mass of ~ 171 m.u. (~ 71 m.u. of WR melt added which is ~ 41 wt% of the magma chamber) with ~ 57 wt% of dacitic (nearly rhyolitic; SiO 2 ≈ 70 wt% at 856 °C) melt left and < 1 wt% of fluid separated (Fig. 2c). The degree of partial melting of the WR is ~ 42 wt%. The notably large degree of WR partial melting and the amount of assimilation are related to the high initial T of the WR and because the simulation proceeds until thermal equilibrium between M and WR has been reached; it takes a lot of crystallization and release of latent heat and thus assimilation of the large WR mass to reach this thermal equilibrium state at ~ 840 °C. The comparatively small amount of total fractionated solid cumulate consists of approximately 8 wt% olivine, 35 wt% clinopyroxene, 38 wt% plagioclase, 5 wt% spinel, 12 wt% orthopyroxene, and 1 wt% rhombohedral oxide (listed in the order of appearance). Note that addition of felsic WR melts increases the stability of orthopyroxene at the expense of olivine relative to the FC run. The residual WR during the run consists of varying amounts of orthopyroxene, plagioclase, alkali feldspar, quartz, spinel, rhombohedral oxide, and apatite. The WR phases that contribute most to the partial melting are alkali feldspar (disappears completely from the WR at ~ 782 °C) and quartz.
The default AFC run shows drastic differences relative to that of the default FC run. Rather surprisingly, the compatible element Ni is higher in the resident melt in the AFC model than in the FC model after assimilation begins (Fig. 3a). This is because the crystallization of olivine with high K M sm (Ni) is subdued in the AFC simulation and because the WR anatectic melt batches have higher Ni (~ 30 ppm) than the evolved resident melt in the FC model (< 30 ppm below 1000 °C). Strontium is also enriched in the AFC resident melt (Fig. 3b) at any given SiO 2 content after assimilation begins even when D WR sm (Sr) remains relatively high (mostly above 1, up to 4; Fig. 6d). This enrichment in the resident melt is primarily caused by subdued plagioclase crystallization in AFC relative to FC (~ 16 wt% vs. ~ 28 wt% of the magma chamber, respectively). In addition, the WR anatectic melt batches introduced into the resident melt become increasingly Sr-rich (from ~ 108 to ~ 139 ppm) as the simulation proceeds.
Incompatible trace elements ( D sm < 1 in both M and WR) that are enriched in the WR (i.e., Zr, Ba, La) are considerably enriched in the AFC resident melt compared with the FC resident melt after the onset of assimilation in the former (Fig. 3c-e). The contrary is true for Yb that is found in relatively lower quantities in the WR (Fig. 3f). Assimilation has drastic effects on REE ratios and isotope compositions of the resident melt (Fig. 4). Note that compatibility of Sr in the WR (Fig. 6d) results in a delayed increase in 87 Sr/ 86 Sr at a given 143 Nd/ 144 Nd in comparison to a simple binary mixing model (Fig. 4b). Similarly, the incompatibility of Nd and its enrichment in the WR anatectic melt decreases 143 Nd/ 144 Nd more efficiently than δ 18 O is increased in the resident melt relative to a binary mixing model (Fig. 4c). The opposite is true when comparing 87 Sr/ 86 Sr and δ 18 O, the latter showing a stronger increase in the resident melt (Fig. 4d).
The effects of using constant "maximum" or "minimum" K ol∕m (Ni), K cpx∕m (Ni), and K plag∕m (Sr) values compiled in Rollinson (1993) are compared to those of the default AFC model in Fig. 6. Varying these values in the model input Fig. 6 Comparison of the resident melt Ni and Sr concentrations (a-c) and D sm (Sr) (d; for both M and WR) of the AFC trace element simulations using default (same as in Fig. 3), constant maximum (AFC-hK), or constant minimum (AFC-lK) K ol∕m (Ni), K cpx∕m (Ni), and K plag∕m (Sr) values. In (d), D sm (Sr) of 1 is highlighted with a dashed line, and the thermal evolution of the resident melt and WR are indicated. The oscillating pattern for D M sm (Sr) in the AFC runs reflects the working principles of the primary MCS tool after assimilation begins: low points in the pattern are first steps after assimilation of WR anatectic melt, whereas each high point in the diagram represents a normal 5 °C decrement step. See text for more details Fig. 7 Ni (ppm) versus Sr (ppm) in the resident melt, bulk cumulate, WR melt, and WR bulk residue (solid + melt + fluid) of the default AFC model. Initial parental melt (PM) and WR compositions are labeled. Arrows indicate the direction of evolution in the MCS simulation. The causes for the major inflection points are marked for the resident melt and bulk cumulate composition causes drastic differences to the outcomes of the otherwise uniform simulations. For example, the final Sr concentration in the resident melt of the "maximum K sm " simulation is ~ 9 ppm, considerably lower than in the default run (~ 102 ppm) or in the "minimum K sm " run (~ 314 ppm).
The jigsaw pattern for D M sm (Sr) (Fig. 6d) reflects the working principles of the primary MCS tool: low points in the jigsaw pattern are the result of crystallization in relation to assimilation of the WR melt that causes decreased plagioclase stability and thus low D M sm (Sr), whereas each high point in the diagram represents a normal 5 °C FC step with increased plagioclase stability and thus high D M sm (Sr). The reason for the continuous depletion of Sr in the resident melt in the default AFC run (Fig. 6b)-in spite of D M sm (Sr) being less than 1 at some temperatures (between ~ 1060-1010 °C; Fig. 6d)-is due to the dilution effect of adding relatively Sr-poor WR anatectic melt into the resident melt (Fig. 7). The inflection points in the AFC-T D WR sm (Sr) trends at around WR temperature of 780 °C (Fig. 6d) are caused by the disappearance of alkali feldspar and increased consumption of plagioclase both of which possess relatively high K sm (Sr) in the residual wallrock (plagioclase higher than alkali feldspar; Table 3). Again, all these are examples of complex relationships between the different processes and variables included in MCS, and are natural consequences of open system behavior that are counterintuitive to the reasoning that either implicitly or explicitly assumes closed system behavior.
Finally, Fig. 7 illustrates the compositional evolution of not only the resident melt, but also the bulk cumulate, WR melt, and WR residual (solid + melt + fluid) in the default AFC run. Nickel concentration in the bulk cumulate is highest in the beginning when olivine ± clinopyroxene with high K M sm (Ni) are the sole fractionating phases. Strontium content in the bulk cumulate begins to notably increase after the addition of plagioclase in the fractionating assemblage. Because Ni and Sr are compatible in WR solids, the WR partial melt shows little change in its composition, whereas the WR residual gets progressively enriched in these elements. Note that all the data presented in Fig. 7 can be collected from MCS-Traces output with minimal effort. Correlating specific trends in trace element evolution is simple, as a major element and phase abundances for each step of the scenario are output alongside the trace element concentrations.

Case 4: assimilation of stoped wallrock blocks-fractional crystallization (S 2 FC)
The S 2 FC model simulates complete assimilation of two stoped wallrock blocks-or two sets of multiple small blocks each with a combined mass equal to the respective S event-by the resident melt. The stoped blocks are of WR composition, and are introduced into the magma using the recharge function at resident M temperatures of ~ 1014 °C and ~ 907 °C. The masses of the blocks are 17 and 38 m.u.. These masses roughly correspond to the total mass of assimilated WR partial melts in the AFC model at the given magma temperatures to enable comparison. The temperatures of the blocks at homogenization are 760 °C (composed of ~ 78 wt% solid, ~ 20 wt% melt, and ~ 2 wt% of fluid) and 795 °C (composed of ~ 60 wt% solid, ~ 38 wt% melt, and ~ 2 wt% of fluid), respectively. These temperatures are above WR liquidus to enable MELTS engine to calculate their phase state at the time of stoping. Note that the total mass of the stoped blocks cannot represent the total mass of WR partial melts assimilated in the AFC scenario (~ 71 m.u.), because wholesale melting of the rather cold stoped blocks requires relatively lot of heat. Since the WR blocks are completely melted and homogenized in M, the decrease in temperature and mass of crystallization in M after each stoping event is considerable (~ 47 °C and ~ 15 m.u. and ~ 45 °C and ~ 25 m.u., respectively).
The S 2 FC model terminates after ~ 68 wt% of crystallization at a total magma chamber mass of 155 m.u. (55 m.u. of stoped WR assimilated which is ~ 35 wt% of the magma system) with ~ 31 wt% of rhyolitic (SiO 2 ≈ 72 wt% at 840 °C) melt left and ~ 1 wt% of fluid separated (Fig. 2d). The bulk solid cumulate consists of approximately 7 wt% olivine, 28 wt% clinopyroxene, 48 wt% plagioclase, 6 wt% spinel, 11 wt% orthopyroxene, and < 1 wt% rhombohedral oxide (listed in the order of appearance). Although WR blocks are composed of solids, melt, and fluid at hybridization, their phase identity is irrelevant for the outcome of the MCS-Traces modeling as the R function in MCS produces a new, completely homogenized and equilibrated state. Similar to the AFC case, the addition of Si-rich WR increases orthopyroxene stability at the expense of olivine relative to the FC run.
Relative to the R 2 FC and FC runs at a given SiO 2 content, S 2 FC resident melt shows varying degrees of enrichment of Sr, Ba, and La (elements notably enriched in the WR), and depletion of Yb (element notably depleted in the WR; Fig. 3). The complete homogenization of a stoped WR block in the resident melt results in a significant compositional change for most elements and isotope ratios immediately after each stoping event (Figs. 3, 4). Interestingly, the SiO 2 vs. Zr trend of the S 2 FC run does not considerably differ from those of R 2 FC or FC (Fig. 3c), because Zr and SiO 2 are equally enriched in the stoped WR blocks relative to the resident melt. The rather small differences in the respective SiO 2 vs. Ni trends (Fig. 3a) are caused by the compatibility of Ni in the crystallizing minerals lowering the Ni concentrations in the resident melt close to that of the WR blocks (47 ppm; Table 2).
Relative to the AFC run, S 2 FC shows some differences that are related to the contamination process. In resident melt that has experienced wholesale assimilation of stoped blocks, Sr, which is compatible during WR melting, gets relatively enriched in the resident melt following assimilation by stoping (Fig. 3a, b), whereas elements that are enriched in the WR and behave incompatibly (e.g., Zr, Ba, and La; Fig. 3c-e) get less enriched; the comparisons are again made at a given SiO 2 content after the beginning of assimilation. On the other hand, incompatible element Yb is depleted in the S 2 FC resident melt relative to the AFC run right after the stoping events because of the low concentration of Yb in the bulk WR (1.96 ppm; Table 2) compared to that of resident melt before the stoping events (3.0 and 4.4 ppm, respectively). The differences between AFC and S 2 FC arise because contamination by stoping represents the incorporation of the bulk WR, rather than the assimilation of anatectic melt, in which compatible elements are relatively depleted and incompatible elements enriched in relation to bulk WR. The degree of enrichment in the assimilated WR anatectic melt batches is dictated by the thermodynamically derived degree of melting and the userinput percolation threshold (FmZero) and K WR sm and K WR sf values. For isotopes, 87 Sr/ 86 Sr and δ 18 O are higher at a given 143 Nd/ 144 Nd in the S 2 FC resident melt than in the AFC run or in a binary mixing scenario (Fig. 4b, c). In the case of 87 Sr/ 86 Sr, this is both because Sr is progressively depleted from the resident melt by fractional crystallization of plagioclase, and because assimilation of bulk stoped WR adds more Sr to the resident melt (per unit addition of added mass) than relatively Sr-poor anatectic melt does. Strontium released from the assimilated WR blocks thus changes the Sr isotopic composition of the melt more effectively than in the case of binary mixing or AFC. The explanation for the higher δ 18 O at a given 143 Nd/ 144 Nd relative to the binary mixing scenario is less intuitive but stems from the fact that Nd is enriched in the resident melt at the time of the stoping events relative to the parental melt used in the binary mixing model. Therefore, 143 Nd/ 144 Nd of the S 2 FC resident melt is not as easily affected by assimilation of stoped WR. The higher δ 18 O of S 2 FC resident melt relative to that of the AFC model at similar degrees of contamination (e.g., around 35 wt% of the total M mass, Fig. 4c, d) can be associated with relatively higher resident melt mass in the AFC simulation; δ 18 O is more efficiently buffered by resident melt composition in a larger melt pool.

Case 5: recharge-assimilation of wallrock anatectic melts-fractional crystallization (R 2 AFC)
In the exemplary R 2 AFC model that combines magma replenishment with the assimilation of WR partial melts, the R events are set as in the R 2 FC model and WR parameters are set as in the AFC model.
Assimilation begins in the R 2 AFC simulation as in the AFC model when the temperature of the resident melt has dropped to ~ 1069 °C (~ 51 wt. SiO 2 ; ~ 26 wt% crystallization of the magma chamber at that point). This precedes recharge events at M temperatures of ~ 1041 °C (~ 54 wt. SiO 2 ) and 998 °C (~ 60 wt. SiO 2 ). The R 2 AFC model terminates after ~ 38 wt% of crystallization at a total magma chamber mass of ~ 352 m.u. with ~ 60 wt% of dacitic (SiO 2 ≈ 63 wt % at 967 °C) melt left and less than 0.02 wt% of fluid separated (Fig. 2e). The large magma chamber mass includes ~ 102 m.u. (~ 29 wt%) of added WR anatectic melt and 150 m.u. (~ 43 wt%) of added R magma. The fractionated bulk solid cumulate consists of approximately 6 wt% olivine, 50 wt% clinopyroxene, 26 wt% plagioclase, 6 wt% spinel, 12 wt% orthopyroxene, and < 1 wt% rhombohedral oxide (listed in the order of appearance). The residual WR mineral composition and behavior is similar to the AFC model, although additional heat provided by the R magmas has melted the WR to a greater degree and has exhausted it of quartz by the end of the simulation.
When assimilation is simulated together with recharge events, the resident melt begins to show complex trace element trends (Fig. 3). These effects are the result of mass exchange between the resident melt and either "depleted" primitive R magmas or "enriched" WR partial melt batches and the resulting changes in phase equilibria and D values that can be reviewed in detail in the related model output. The effects of R events for contamination-sensitive trace element and radiogenic isotope ratios are not drastic, however, and the respective trends of the R 2 AFC run closely correspond to those of the AFC simulation (Fig. 4). What is very different in these models, however, are the masses of the components (see Fig. 2). The final resident melt composition of the R 2 AFC run is the result of assimilation of ~ 102 m.u. of WR anatectic melt by 250 m.u. of initial and R magmas. The resulting compositions correspond to those of the AFC run after ~ 30-35 m.u. of WR anatectic melts has been assimilated by 100 m.u. of initial melt. Note that at the end of the R 2 AFC run, the amount of assimilation relative the resident magma chamber mass is nevertheless similar (~ 30 wt %) to that of the AFC run at a given isotopic composition (Fig. 4).

Discussion on selecting partition coefficients
For any reasonable trace element and radiogenic isotope model of a magmatic system, selection of feasible partition coefficients is crucial. Given the number of studies on the partitioning of elements in different phases in different magma compositions and at different P-T-conditions, it is easy to get overwhelmed in the quest for such information.
For our example models presented here, we used a set of partition coefficients compiled from other modeling software or from individual publications of systems similar to the resident melt and WR partial melt compositions in the MCS models. In addition, we applied T-dependence whenever it was feasible and enough relevant data was available ( Table 3).
The example of using hypothetical "maximum" or "minimum" K sm reported in the compilation of Rollinson (1993) demonstrates the issue of selecting partition coefficients (Fig. 6). Partition coefficients are often considered constant in the still widely used AFC model of DePaolo (1981) and in the EC-AFC class of models (Spera and Bohrson 2001, 2004Bohrson and Spera 2007) that only model trace elements and do not consider phase equilibria. In the MCS models presented here, the resident melt evolves from basaltic to dacitic/rhyolitic and crystallizes and melts a wide range of silicate, oxide, and phosphate phases, which propel considerable difficulties for selecting constant K sm values for the whole run. Therefore, we strongly recommend running MCS-Traces models with T-dependent K input, especially for cases that span large compositional and temperature ranges. If such data are not available, running with different sets of constant partition coefficients (e.g., using relevant minimum and maximum reported K sm values for each phase) could be a viable option. A multifaceted approach will provide a framework for the relative importance of reported K sm variance and strengthen the presented models. The benefit of relying on highly compatible or incompatible elements is that possible magnitude-scale variations reported for K sm values (e.g., from 100 to 10 or from 0.01 to 0.001) will have a minute effect on the resulting geochemical trends. Whatever is the selected approach, it is important to clearly state which K sm values were used, what was the reason they were selected, and what the referred sources are.
Finally, the additional complexity in element partitioning in MCS models is caused by the appearance of a fluid phase, either in the resident magma or in the wallrock (Fig. 5). Unfortunately, K sf values are highly underrepresented relative to K sm values in experimental petrology and petrologic literature in general. This still is an important aspect of the modeling, because if an element is highly fluid-compatible, even small amounts of fluid separation may have drastic effects on its concentrations in the resident melt (e.g., Cs, Rb, Th; for Sr and Ba, see Fig. 5). Geochemical models based upon the assumption that trace elements are insoluble can lead to inaccurate source composition estimates if, in fact, appreciable amounts of trace elements are removed from the system by volatiles. Furthermore, assimilation of the fluid phase from the wallrock is not presently handled by MCS (see Bohrson et al. 2020), so highly fluid-compatible elements may be primarily concentrated in the residual wallrock and not assimilated with the wallrock partial melts in AFC scenarios. Therefore, we recommend using caution in selecting K sf values for the resident magma and especially for the wallrock in modeling fluid-saturated systems.

Trace elements and isotopes as indicators of open-system processes
Trace elements and isotopes have traditionally been regarded as useful indicators of open-system processes in igneous systems. This is because major element and mineral compositional evolution is controlled by phase equilibria, and effects of recharge and assimilation processes on these may be minute or difficult to identify at least from bulk compositions (as illustrated in the companion paper). In the sections below, we review and, on the basis of MCS-Traces modeling, re-evaluate the usefulness of trace elements and isotopes as such indicators.

Identifying fingerprints of magma recharge
Identifying recharge (and complete homogenization) of magma pulses having the same composition with the parental melt based on major element evidence is challenging (Bohrson et al. 2020), and the same is true for incompatible trace elements: Zr, Ba, La, and Yb trends for both the FC and R 2 FC scenarios are practically uniform (Fig. 3c-f). Some variation is evident in the case of Ni, which behaves mostly compatibly throughout the run, as well as Sr, which is incompatible before the onset of plagioclase fractionation (Fig. 3a, b). In both cases, slight offsets from the FC trend can be seen the following recharge. Such offsets may nevertheless be difficult to distinguish from natural sample variation due to other factors (such as sample inhomogeneity) and analytical uncertainty.
We conclude that magma recharge processes are very difficult to identify in differentiating magmatic systems on the basis of major or trace element trends alone if the recharge melts are of uniform composition. Instead, one or more of the following three premises should preferably be fulfilled: (1) the recharge magma has a distinct composition relative to earlier initial or parental melts of the system (e.g., Pietruszka and Garcia 1999); (2) stratigraphic or other temporal sample control enables detection of "reverse" geochemical 105 Page 16 of 21 and/or mineral compositional trends (e.g., increase in MgO or Ni contents) that are best explained by the addition of more primitive magma into the system (e.g., Cox 1988;Blight et al. 2010;Yuan et al. 2017);and/or (3) there is a representative crystal cargo available to study crystalline microtextures and textural relationships, and major and trace element compositional variations. In the case of zoned crystals, this can be regarded as an intracrystalline stratigraphic control, and its importance in indicating magma recharge has been recognized in different magmatic environments, especially those related to subduction settings (e.g., Ginibre and Wörner 2007;Streck 2008;Ginibre and Davidson 2014;Bouvet de Maisonneuve et al. 2016). However, evidence from MCS-PhaseEQ modeling suggests that recharge is not necessarily recorded similarly in all phases (e.g., in those that show varying stability in a replenishing system; Bohrson et al. 2020). In addition, post-crystallization reequilibration of crystals by diffusion may readily erase some elemental zoning in certain minerals (see, e.g., Costa et al. 2008). Note that for studying the effects of recharge on individual phase compositions, mineral major and trace element compositions are also recorded in MCS output.

Identifying fingerprints of crustal assimilation
Crustal assimilation was also surprisingly difficult to recognize based only on major element oxides from MCS model output in standard AFC and SFC scenarios (Bohrson et al. 2020), even though the amount of assimilation in the presented models is very large. Notable differences in the modeled AFC scenario relative to FC or R 2 FC scenarios were only identified in terms of K 2 O (5 × enrichment relative to FC or R 2 FC model due to melting of alkali feldspar in the wallrock) with subtler differences also present in Al 2 O 3 , CaO, and Na 2 O. These differences were even less obvious if the degree of contamination was lower or happened by wholesale incorporation of stoped blocks (S 2 FC case).
In contrast, AFC, S 2 FC, and R 2 AFC cases are readily distinguished from FC and R 2 FC cases in terms of incompatible elements and isotopic compositions (Figs. 3, 4). These will be reviewed in case-specific sections below.

Assimilation of wallrock anatectic melts
Assimilation of wallrock anatectic melts (AFC in MCS jargon) will dominate in settings where wallrock xenoliths are not efficiently introduced into the resident melt, but assimilation takes place within a contact zone of the wallrock and the resident melt. Such conditions largely prevail in middle to lower crustal magma chambers, where the wallrock is not highly fractured, behaves in ductile manner, and boundary temperatures are high (see Bohrson et al. 2014).
MCS modeling shows that partial melting of the wallrock considerably enriches the wallrock partial melt and thus the contaminated resident melt in incompatible elements that are found in high quantities in the wallrock (e.g., Zr, Ba, and La) already during early stages of assimilation (Fig. 3c-e). In fact, the first batches of WR partial melts are most important in enriching the resident melt, since they are loaded with such elements. The melting WR, and thus WR partial melts, will get progressively depleted in incompatible elements after each assimilation step, which then dilutes this effect. The presented MCS models suggest that the relative amounts of assimilation in natural systems on the basis of incompatible trace element concentrations or their isotopes, without some kind of control on the thermodynamics, phase equilibria, and style of assimilation, may often be significantly overestimated (see Heinonen et al. 2016). The contrary is true if elements are incompatible but present in low quantities within the wallrock (e.g., Yb; Fig. 3f). Ytterbium is initially slightly enriched in the contaminated resident melt relative to the FC case, because Yb effectively enters the assimilated wallrock partial melt ( D WR sm (Yb) is 0.1-0.3 throughout the run). Because of the initially low concentration in the wallrock (1.96 ppm relative to 1.73 ppm in PM; Table 2), however, the wallrock partial melts become depleted in Yb quite rapidly relative to, for example, SiO 2 , Zr, Ba, or La.
Compatible element Ni shows unexpected behavior (Fig. 3a): its concentrations in M melt in the AFC case are relatively higher than in the FC case, although assimilation of relatively Ni-poor WR (Table 2) would be expected to dilute its concentrations in the contaminated melt. As stated previously, this surprising effect is because assimilation suppresses the crystallization of olivine with high K sm (Ni) in the AFC case relative to the FC case, and because Ni is not as compatible in WR ( D WR sm (Ni) ≈ 1.7-1.9) as it is in the magma chamber ( D M sm (Ni) ≈ 2-10). The different degrees of enrichment/depletion shown by variably compatible elements in AFC resident melt underline the importance of defining phase equilibria and evaluating the role of assimilation in natural systems based on several and not just one or two different elements.
The isotopic effects of assimilation of wallrock partial melts are clear and most drastic for elements that show incompatible behavior in the wallrock and are relatively depleted in the resident melt (Fig. 6). If the partially melting wallrock contains residual plagioclase in significant proportions such that D WR sm (Sr) > 1, 144 Nd/ 143 Nd is expected to show more drastic changes than 87 Sr/ 86 Sr. Oxygen isotopes have been considered as useful indicators of crustal contributions to mantle-derived magmatic systems because O concentrations in magmas and wallrock are usually very similar whereas their O isotopic compositions are different. Therefore, binary mixing or traditional AFC curves (with O isotopes often considered together with radiogenic isotopes) have been used to estimate the relative amounts of mantle and crust end-members and importance of assimilation vs. source heterogeneity in such systems (e.g., Taylor 1980;James 1981;Martinez et al. 1996;Baker et al. 2000;France et al. 2016). The results of MCS modeling illustrated in Fig. 4 indicate that the differences in the compatibility of the trace elements (both in the magma and in wallrock), O budget of the crystallizing resident melt vs. assimilant, and style of assimilation have considerable effects on such isotopic comparisons in different contamination scenarios. It is evident from Fig. 4 that relative contributions from crustal sources may be considerably overestimated using binary mixing models based on O isotopes.
If magma recharge is concurrent with anatectic wallrock assimilation, the latter process controls enrichment in incompatible elements. We emphasize, however, that this observation is highly dependent on the composition of the recharge magma and the wallrock, and the mass balances among the three subsystems. Finally, in addition to being able to follow the geochemical evolution of the resident melt, MCS records the geochemical evolution of the residual wallrock (Fig. 7). In cases where such information is available (e.g., contact zones of intrusions, inclusions of anatectic melt), MCS can be used in accordance with other information on the magmatic system to confirm or refine scenarios involving wallrock partial melting.

Wholesale assimilation of stoped blocks
Wholesale assimilation of stoped blocks (SFC in MCS jargon) will dominate in settings where wallrock xenoliths are efficiently introduced into the resident melt. Such conditions may prevail in upper crustal magma chambers or feeding channels, where the wallrock is highly fractured and behaves in a brittle manner, and boundary temperatures are low (see Bohrson et al. 2014). The smaller the stoped blocks are, the higher the chance of rapid thermalization and partial or wholesale melting and, eventual equilibration with resident melt because small blocks heat up internally faster than larger ones (Carslaw and Jaeger 1959). Note that a single S event in MCS can be considered as stoping of multiple small xenoliths, the combined mass of which equals the total mass of the S event.
When wallrock is completely homogenized and equilibrated with the resident melt (i.e. bulk assimilation), partition coefficients for the wallrock phases do not play a role: i.e., the effect of assimilation on trace element contents of the resident melt is fully dictated by their concentration in the wallrock and the mass of the assimilated wallrock block. Therefore, in assimilation by stoping, any element that is relatively enriched in the assimilated wallrock block will also be enriched in the resident melt. Such effects may subsequently be diluted by fractional crystallization, however.
The isotopic model results very clearly illustrate the differences between assimilation of stoped blocks and anatectic melts (Fig. 4). Whereas Nd is incompatible and Sr compatible in wallrock in the case of wallrock partial melting (AFC; see previous section for details), in S 2 FC, the isotopic composition of the resident melt is only controlled by bulk assimilation of the wallrock blocks into the resident melt. Depending on the compatibility of the associated elements and the process of contamination, the curves in isotope diagrams can be convex upwards or downwards (cf. Sr and Nd isotopes relative to δ 18 O in AFC and S 2 FC in Fig. 4b, c), making the distinction between wallrock contamination vs. source heterogeneity rather difficult (cf. Taylor 1980;James 1981).
If AFC and S 2 FC are thought of as end-members of the assimilation process, anything in between could take place in natural systems, where there is always a middle ground between wholesale melting of stoped blocks and partial melting of the wallrock. In most cases, MCS modeling would benefit from considering both options. Note that it is also possible to model SFC and AFC together in a single run to simulate such a system.

Model assumptions and limitations and future development of the MCS-Traces
MCS is a versatile, mass-and energy-balanced thermodynamic tool that can model complex scenarios that involve R n AS n FC. One of its shortcomings is that it is "only" an equilibrium thermodynamic model; it, therefore, does not account for time-dependent kinetic, diffusion, or transport phenomena, or other processes that may result in chemical disequilibrium-a well know feature of most igneous systems. Many of these limitations have been discussed in the companion paper (Bohrson et al. 2020) and the reader is referred to it for more information, but we will further concentrate on those that are particularly relevant to modeling trace elements and their respective isotope systems.
A potentially important limitation in terms of trace element modeling is related to the MELTS engines. Current versions do not incorporate some minor phases such as zircon or monazite. Such phases are not usually important in terms of major element evolution of igneous systems, but they serve as important hosts for many otherwise incompatible trace elements (e.g., in the case of zircon and monazite, Zr and Th + REE, respectively). For example, Zr is modeled here as a generally incompatible trace element in the wallrock, but in an average crustal granitoid, zircon is almost always present, and it will thus influence the behavior of Zr during wallrock partial melting. If zircon remains in the wallrock residue, wallrock partial melts will be relatively depleted in Zr, which will, in turn, affect the Zr evolution of the evolving resident melt. We recommend that the possibility of minor phases retaining some trace elements in the wallrock residue should be evaluated in MCS scenarios. In addition, MELTS engines provide limited thermodynamic data for water-bearing phases (e.g., biotite and amphiboles; see Gualda et al. 2012;Bohrson et al. 2020) which are often not stabilized by the engines in systems where they are expected. This could have a negative impact on comparisons of MCS-Traces models to fluid-rich natural systems. If fluidbearing phases are common in the studied natural system, trace elements incompatible to them should be preferred in the modeling until more robust thermodynamic data is available for the MELTS engines.
A designed future extension of MCS-Traces is to allow for isotopic disequilibrium during the partial melting of wallrock as demonstrated by Iles et al. (2018). In its current form, wallrock has an isotopic ratio that is the same as all phases in the wallrock. In the case of radiogenic isotopes, for example, this implies that differential ingrowth has not occurred in any of the minerals, and hence the isotopic ratio in anatectic melt does not change during assimilation. However, pre-magmatic ingrowth of daughter isotopes in wallrock phases, dictated by phase-specific parent-daughter ratios and time leads to variations in the respective isotopic ratios among mineral phases. In addition, although radiogenic ingrowth does not concern O isotopes, O isotopic heterogeneity between mineral phases in the same rock is also a common feature. It is easy to envision, for example, an Archean granitic wallrock where coexisting minerals (e.g., biotite, alkali feldspar, and plagioclase with distinct Rb/Sr) have distinct 87 Sr/ 86 Sr. As wallrock undergoes partial melting, partial melt will have 87 Sr/ 86 Sr that is a function of the contribution each crystalline phase makes to the melt. Hence the isotopic composition of wallrock melt will change during the partial melt assimilation process. The potential effects of isotope fractionation between an isotopically heterogeneous wallrock and its partial melt batches can currently be approximated by, for example, running two models using representative minimum and maximum radiogenic isotope ratios for the wallrock. Extension of MCS to include disequilibrium isotopes would allow a more quantitative and high-fidelity approach to this issue.
Similar to MCS-PhaseEQ, there is an urge to translate MCS-Traces to a programming language that would combine the two and be platform-independent. While the advantages to the two-stage calculation presented here are many, merging the two tools would enable the Monte Carlo approach to trace element and isotope calculations, as has been proposed for the major element and phase equilibria calculations (see Bohrson et al. 2020). This modeling environment would allow users to produce a large number of models relatively rapidly, and with less human bias in choosing parameters. Future improvements will also likely include statistical tools for enhanced handling of the model results. The users are encouraged to check the MCS website for software updates.
Finally, one of the most important things to realize as a petrologist and a potential user of MCS is the scale and the relative importance of the involved processes. If some crystals in the studied natural system are zoned, full of solid inclusions, or exhibit any disequilibrium textures, this does not instantly mean that MCS could not be used to constrain the general magmatic evolution of the system. Although current versions of MCS assume thermodynamic equilibration of a magma during open-system processes, the processes which produce such disequilibrium textures, such as recharge and magma mixing, are those which MCS was designed to model, and thus can reasonably approximate. In addition, with the ability to record major and trace element concentrations of melt, fluid, and each mineral in every step, MCS can provide important constraints on the thermal histories and geochemical conditions during the generation of crystal zoning; of course, interpretations of model results must be placed in context and compared with those derived from other data sources. In other words, when attempting to understand a natural system, one can search for the best MCS model by perturbing initial conditions and parameters iteratively, while remaining within the realm of reasonability. Such a 'best' model will, in all likelihood remain imperfect, perhaps in significant ways. Such imperfection should not be seen as a failure. Instead, by struggling with the root cause of the imperfections, valuable clues to the systems petrological and geochemical evolution may be elucidated. That is, the nature of the transport, kinetic and dynamic phenomena, precisely those ignored in thermodynamic models, may rise to the top to be addressed by 'beyond thermodynamic' analysis. This confrontation with reality and attempt to model complex systems is the ultimate goal of petrology. We close the MCS companion papers citing the final few words in N.L. Bowens opus on the Evolution of the Igneous Rocks (Bowen 1928): "No implication is intended that our knowledge is all that could be desired. That can never be."

Summary
The Magma Chamber Simulator MCS-Traces computer program calculates the simultaneous evolution of 48 trace elements, seven radiogenic isotope systems (Sr, Nd, Hf, 3xPb, and Os), and O isotopes for a magma (melt + crystals ± fluid phase) in an igneous system influenced by magma recharge, assimilation (either of wallrock partial melts or stoped blocks), and fractional crystallization (R n AS n FC). In addition, the trace element evolution of the cumulate pile and residual wallrock is output. MCS-Traces requires output from the MCS-PhaseEQ computational tool, which is described in a companion paper (Bohrson et al. 2020), to perform the calculations. Our analysis of example case studies of depleted mantle magma crystallizing ± experiencing recharge pulses ± assimilating average granodioritic crust at 0.1 GPa reveals key characteristics of the different processes. For the case studies presented here, the effects of magma recharge are difficult to recognize from trace element or radiogenic isotope data alone; a temporal (e.g., stratigraphic or intracrystalline) reference frame is required. We also note that recharge may be more evident when there is greater contrast between resident melt and intruding magma(s). In contrast, assimilation (especially in the case of assimilation of wallrock anatectic melts) may have drastic effects on incompatible trace element and isotopic compositions of the resident melt. The magnitude of these effects depends predominantly on phase equilibria, element partitioning, style of assimilation, and the geochemical differences between the resident melt and the wallrock. In many cases, the effects of assimilation on trace elements and isotopes are counterintuitive and considerable caution should be practiced in ruling out potential assimilation scenarios. The case studies presented here highlight the challenges that petrologists and geochemists face when diagnosing open-system processes. Discrepancies between observations and the 'best' open system model(s) have the potential to reveal details of the role of kinetic and transport phenomena in petrogenesis. There are plethora of ways to use MCS to solve petrogenetic questions related to magmatism in a wide range of igneous environments.