Petrological Geodynamics of Mantle Melting II. AlphaMELTS + Multiphase Flow: Dynamic Fractional Melting

In this second installment of a series that aims to investigate the dynamic interaction between the composition and abundance of the solid mantle and its melt products, the classic interpretation of fractional melting is extended to account for the dynamic nature of the process. A multiphase numerical flow model is coupled with the program AlphaMELTS, which provides at the moment possibly the most accurate petrological description of melting based on thermodynamic principles. The conceptual idea of this study is based on a description of the melting process taking place along a 1-D vertical ideal column where chemical equilibrium is assumed to apply in two local sub-systems separately on some spatial and temporal scale. The solid mantle belongs to a local sub-system (ss1) that does not interact chemically with the melt reservoir which forms a second sub-system (ss2). The local melt products are transferred in the melt sub-system ss2 where the melt phase eventually can also crystallize into a different solid assemblage and will evolve dynamically. The main difference with the usual interpretation of fractional melting is that melt is not arbitrarily and instantaneously extracted from the mantle, but instead remains a dynamic component of the model, hence the process is named dynamic fractional melting (DFM).\\ Some of the conditions that may affect the DFM model are investigated in this study, in particular the effect of temperature, mantle velocity at the boundary of the mantle column. A comparison is made with the dynamic equilibrium melting (DEM) model discussed in the first installment. The implications of assuming passive flow or active flow are also considered to some extent.\\ Complete data files of most of the DFM simulations, four animations and two new DEM simulations (passive/active flow) are available following the instructions in the supplementary material.\\


INTRODUCTION
Melting in the Earth's and planetary interiors is usually approached either from a physical perspective based on geophysical tools (e.g., seismology) and numerical studies based on dynamic two-phase flow models, or from a petrological perspective that focuses mainly on the chemical evolution of the melt products and minerals. In the first installment (Tirone and Sessing, 2017), hereafter called MF1, a petrological and dynamic numerical model has been introduced with the aim to develop a dynamic and petrological representation of the melting process in the Earth's mantle. The general motivation was that the classical petrological descriptions of melting, namely batch melting and fractional melting are physically flawed. Batch melting implies that melt products stay always in contact with the source, an implausible scenario for mantle melting, except for the cases in which only tiny amount of melt pockets are the end-product. Fractional melting would imply that no melt would exist in the mantle since all melt is immediately extracted. Both are also static atemporal models, the mantle doesn't move and the melt either doesn't move or disappears instantaneously. Time is also irrelevant. To remove these unreal physical limitations the dynamic model, exemplified by a 1-D mantle column, has been formulated considering the presence of two or three dynamic phases (solid, melt or solid, melt and water). The petrological and chemical description is provided by the program AlphaMELTS (Smith and Asimow, 2005) (pMELTS version) which is based on the thermodynamic model developed by Ghiorso over several years (Ghiorso, 1985;Ghiorso and Sacks, 1995;Ghiorso et al., 2002). The ability of the earlier version of MELTS to reproduce observable and experimental data was evaluated in a series of studies (Hirschmann et al., 1998(Hirschmann et al., , 1999aAsimow et al., 2001). A thermodynamic model is not an absolute requirement. A petrological tool based on melting parameterization has been developed recently (Brown and Lesher, 2016), and in principle it could be used in combination with some dynamic model as well. Numerous studies investigated the dynamic evolution of melt and solid mantle using a twophase flow model (McKenzie, 1984;Stevenson, 1984, 1986;Spiegelman, 1993;Richardson, 1998;Ghods and Arkani-Hamed, 2000;Schmeling, 2000;Bercovici et al., 2001;Bercovici and Ricard, 2003;Spiegelman et al., 2007;Hewitt and Fowler, 2008;Katz, 2008;Hewitt, 2010;Rudge et al., 2011;Oliveira et al., 2018). The most recent one (Oliveira et al., 2018) incorporates a complex rheology and a realistic thermodynamic equilibrium model.
The numerical solution and the details of the coupled petrological and geodynamic 1-D model applied here have been presented in the previous study, MF1. The concept of dynamic equilibrium melting (DEM) was also introduced in MF1. Essentially it assumes chemical, thermal and mechanical equilibrium between solid and melt over a certain spatial and temporal interval. Even though some form of hybrid kineticthermodynamic petrological modeling can be conceived , and a dynamic melting formulation based on simplified kinetic principles has been considered already (Rudge et al., 2011), this series of contributions focuses exclusively on the domain defined by thermodynamic equilibrium principles. The reason is that as soon as the petrological batch and fractional melting models are put in a real physical context, the spatial and temporal extension of chemical equilibration can still be interpreted rather freely. Hence it seems reasonable to devote first some effort trying to understand the consequences of different possible equilibrium scenarios. One of these scenarios may be described by a model in which the melt cumulative products are chemically isolated from the local residual solid. It is a common practice in petrology to consider this fractional melting process to explain and model certain geochemical and petrological evidences (e.g., Langmuir et al., 1992;Herzberg and Asimow, 2008;Hole, 2015;Rey, 2015;Carr and Gazel, 2017). However the concept of fractional melting is usually applied without any assumption on the temporal and spatial constraints. Therefore a logical step forward would be to investigate the outcome of a fractional melting process that develops dynamically inside the mantle and evolves over time under certain equilibrium assumptions. This would be also the simplest description of melting departing from the assumption of complete local equilibrium between solid and melt. Such model called dynamic fractional melting (DFM) is the main subject of this study. Asimow (2002) developed the conceptual original idea and first implementation of dynamic fractional melting. The work presented here could be considered an extension of this earlier model in the sense that it removes certain restrictions of the previous implementation. For instance in the current formulation the properties are time dependent, the flow pressure is computed explicitly (relative to a fixed point), the melt velocity is not assumed to be much greater than the solid velocity and no restriction is imposed on the density of melt. This study also presents the complete set of compositional/petrological results and applies the latest thermodynamic melt model for mantle melting, pMELTS (Ghiorso et al., 2002), that at the time of the original development was not available.
The main features and the recipe for a practical implementation of a petrological dynamic melting model have been presented in MF1, therefore mainly new aspects of the DFM model will be discussed in the following sections. The DFM model is presented in the next section "Description of the Dynamic Fractional Melting Model". It is followed in section "Changes in the Thermodynamic Computation" by a discussion of the modifications that need to be implemented to the thermodynamic computation. The section "Changes to the Mass Conservation Equations" begins to outline the changes to the transport model, and includes a brief summary of the transport equations.
The results of the DFM model are first compared with the equilibrium dynamic melting model DEM presented in MF1. Only the effect of certain parameters are then investigated, in particular the temperature at the base of the mantle column, the upwelling velocity of the solid at the boundary and the assumption of passive and active flow approximately represented by a fixed velocity either at the top or bottom boundary. All the results are presented in the section "Results of the Dynamic Fractional Melting Model".

DESCRIPTION OF THE DYNAMIC FRACTIONAL MELTING MODEL
The general idea of the DFM model is illustrated in Figure 1. At every grid point the representative space is occupied by two sub-systems, sub-sys 1 (ss1) and sub-sys 2 (ss2), each of the two developing separate chemical equilibration pathways. The first sub-system is only filled by the solid mantle (and a free water phase, if present in the system), while the second subsystem consists of melt and a different solid that may form by partial crystallization of the local melt. When the thermodynamic computation applied to the solid in ss1 predicts that partial FIGURE 1 | Schematic view of the mantle melting column. In the dynamic fractional melting (DFM) model at each grid point local equilibrium is assumed separately for the residual solid and for melt in two sub-systems. Dynamic and thermodynamic properties represent the average spatial values within each sub-system. Melt formed in the first sub-system (sub-sys1) is transferred in the second sub-system (sub-sys2). Melt in sub-sys2 is allowed to crystallize but the new solid remains in sub-sys2. Distribution of melt and minerals in the figure not in scale. The illustration describes only one possible visual interpretation of the model regarding the melt and solids distribution inside the representative spatial grid points. Details about the numerical discretization and numerical solution can be found in MF1. melting will occur, the melt is transferred to the second subsystem associated to the same depth location and it will not interact chemically any further with the residual solid. In this way the first sub-system is always free of melt. The second sub-system consists essentially of a cumulate melt, however chemical thermodynamics is applied to this sub-system as well to determine whether partial crystallization may take place or whether a free water phase would form at equilibrium. If new solid crystallizes, it remains in ss2 while any amount of free water is transferred in ss1. Both the aggregate melt and the fractionated solid do not interact chemically with the solid in ss1. However from a dynamic point of view, the two solids are assimilated to a single transport phase. Melt is treated as a separated dynamic phase, while the free water phase (if present) is considered the third dynamic phase.
The following sections explain how this process is modeled. The discussion in particular highlights the changes that have been applied to the previous formulation introduced in MF1.

Changes in the Thermodynamic Computation
The main difference in the thermodynamic approach between the DFM and the DEM model is that chemical equilibrium is maintained in two sub-systems separately, that is two separate thermodynamic computations are performed at every spatial grid point. The temperature and the pressure are assumed to be the same in the two sub-systems but the bulk composition is different. As in the previous study, the pressure applied in the thermodynamic computations is the lithostatic pressure and it is assumed to be the same in ss1 and ss2. The energy transport equation is solved for the temperature of the whole system (ss1+ss2) and does not include the viscous dissipation effect (small contribution, Asimow and Stolper, 1999;Tirone, 2016) and the latent heat of melting/crystallization. The latter has been neglected mainly because the numerical implementation presents challenges that cannot be easily addressed yet (Tirone et al., 2012).
The bulk abundance in wt % of a certain component c (oxide) for example in sub-system 1 at grid point i and time t, is computed using the following relation: where the superscripts s, m, w indicate solid, melt and free water. θ sc ss1 i,t is the wt % of the component c defining the bulk composition of the solid in ss1, s ss1 i,t is simply the product (ρ s ss1 i,t φ s ss1 i,t ) where ρ s ss1 i,t is the mass density of the solid and φ s ss1 i,t is the fraction of the control volume occupied by the solid. And last sc ss1 i,t is the short notation for the product (ρ s ss1 i,t φ s ss1 i,t θ sc ss1 i,t ). Since the free water phase (if present) is made only of H 2 O, it follows that for example if c stands for MgO, then in the numerator of Equation (1), w i,t should be set to zero. In sub-system 2 a similar expression can be applied, bulk component c ss2 = 100 ( sc ss2 i,t + mc i,t )/( s ss2 i,t + m i,t ). It can be noted that the superscripts ss1 and ss2 are used only for the solid components, while melt is only located in sub-system 2 and water (if present) only in sub-system 1. The data extracted from the thermodynamic computation are the same that were discussed for the DEM model in the previous study. However in this case they are specific to a particular sub-system: total mass and wt % of oxides of the solid, thermal expansion of the whole sub-system, mass, volume, oxide bulk composition and heat capacity at constant pressure of the mineral components and similar quantities for melt and water. Heat capacity and density are then computed for the total solid, melt and water in each sub-system following a procedure similar to the one described in MF1.
Another significant difference with the previous model is that the computation of the volume fraction of an assemblage (solid in ss1 and ss2, melt in ss2, and water in ss1) with respect to the total volume of the two sub-systems does not simply require the volume of each assemblage, because the thermodynamic computations are applied to sub-systems occupying different portions of the total space associated to a particular grid point. The volume fraction of the solid in ss1 can be retrieved from the following expression: where the various are taken from the transport model and the volumes are computed from the thermodynamic model. Note that the volumes include solid melt, and water in both sub-systems before they are rearranged (e.g., newly formed melt transferred in ss2). Similar relations can be used to determine the volume fraction of solid in ss2, melt in ss2 and water in ss1, given that the sum must be equal to 1, φ s ss1 + φ s ss2 + φ m + φ w = 1. The mass transfer is then computed for each sub-system following the procedure discussed in MF1. For example knowing the density and volume fraction of the solid in ss1 from the thermodynamic model, the mass transfer is: and similarly for the solid in ss2. Since no melt is present in ss1 before the thermodynamic computation is applied (all previous melt was already transferred to ss2), the mass transfer for melt in ss1 is defined as m ss1 =[(density of melt × volume fraction of melt) ss1 ], which is always ≥ 0. Melt in ss2 can also crystallize, therefore the melt transfer in ss2 is m ss2 =[(density of melt × volume fraction of melt) ss2 − m ]. For water in ss1 m ss1 =[(density of water × volume fraction of water) ss1 -m ] and when water is formed in ss2, m ss2 =[(density of water × volume fraction of water) ss2 ] since it was not present before the thermodynamic computation (all water was transferred already in ss1). The chemical mass transfer quantities follow the same principles, for example for the transfer of the chemical component c related to the solid in sub-system 1:

Summary of the Transport Equations
Before describing the changes in the dynamic model, the transport equations presented in MF1 are briefly summarized here. The conservation of mass for the total solid, melt and water separately is given by similar expressions assuming this general form: where v a is the vertical velocity of the assemblage a and the mass transfer term Ŵ a is not explicitly considered in the numerical solution model but instead it is incorporated in the mass at the previous time step. The equation of motion for the solid assemblage (ss1+ss2) is expressed as follow: where S encloses the description of the viscous forces which depend on the velocity of the solid. The M terms describe the dynamic coupling of motion of the assemblages, in general they are proportional to the difference of the velocities of two assemblages (e.g., s − m ⇒ solid − melt). The equation of motion for melt is given by the following expression: after replacing M s−m with an expression that depends on the permeability, melt volume fraction and velocity difference, this equation is solved explicitly for the melt velocity v m . For practical reasons the equation of motion for the solid assemblage (Equation 6) is not used in the numerical solution but a sum of the equations of motion for the three assemblages is applied instead: The three mass conservation equations (similar to Equation 5), the equation of motion for melt and water (Equation 7 for melt and a similar one for water), the total equation of motion (Equation 8) and the volume fraction constraint (φ s +φ m +φ w = 1) are solved simultaneously by iteration following the procedure detailed in MF1 to retrieve the pressure (or pressure gradient), the volume fractions φ s , φ m , φ w , and the respective velocities v s , v m , v w . For every chemical components c in the solid and melt the transport equation is expressed in this form: where the chemical mass transfer term (e.g., Ŵ sc−mc ) is absorbed in the chemical mass component at the previous time step in the same way that Ŵ a was implicitly manipulated in the mass conservation equations. The temperature is described by the following relation: where C p a is the heat capacity at constant pressure of the assemblage a and K sys and α sys are the thermal conductivity and thermal expansion of the whole system and v sys ≈ a v a φ a . The last two terms in the above equation describe the reversible adiabatic effect due to compression or expansion. The effect of viscous dissipation has been ignored because it has a minor effect on the thermal profile over the considered depth interval (less than hundred kilometers). The heat change induced by chemical transformations (mainly melting) has been also excluded from the thermal model even though its overall contribution can be on the order of tens of degrees.

Changes to the Mass Conservation Equations
These transport equations are applied to the DFM model as well but some changes are necessary to account for the different interpretation of the DEM and DFM models. For instance considering the melt phase, it is not allowed to remain in the first sub-system and melt in the second sub-system can also crystallize. As in the DEM model, the melt transfer is added to the related term from the previous time step m i,(t−1) , however in this case it should include two quantities, one to describe the transfer from ss1 to ss2 ( m ss1 ) and one for the change due to crystallization in ss2 ( m ss2 ): where nthermo is the number of time steps between two thermodynamic calculations. The definition of the two m terms have been given in the previous section. The same numerical solution of the mass conservation equations for melt and water (Equation 5) outlined in MF1 is applied here. In the transport model the volume fraction of solid was retrieved from φ s + φ m + φ w = 1. In the DFM model a second solid could be present in sub-system 2, therefore φ s represents the sum of the solid volume fraction in ss1 and ss2, φ s = φ s ss1 + φ s ss2 . To find the volume fraction of the solid in each of the two sub-systems (φ s ss1 and φ s ss2 ) after computing the total φ s from the overall solution of the transport model, an additional step is needed. The ratio: is defined at the previous time step and includes the mass transfer terms. The assumption is that R remains unchanged at the current time. Recalling the definition of = ρφ, then for the solid in ss1, s ss1 i,t = ρ s ss1 i,t φ s ss1 i,t and replacing s ss1 i,t with the ratio expression, R s ss2 i,t = ρ s ss1 i,t φ s ss1 i,t , the volume fraction φ s ss1 i,t can be found using: By replacing φ s ss1 i,t with φ s −φ s ss2 in the above equation, the volume fraction in ss2 can be expressed as: With the above equation the volume fraction φ s ss2 i,t can be computed knowing φ s i,t from the transport model and the volume fraction constraint (1 − φ m i,t − φ w i,t ), and the densities ρ s ss2 i,t , ρ s ss1 i,t from the thermodynamic model. Finally φ s ss1 i,t , s ss1 i,t and s ss2 i,t are then easily retrieved.

Changes to the SIMPLER and SIMPLECR Algorithm
The modifications to the SIMPLER or SIMPLECR algorithm are minimal because the velocity of the solids in the two sub-systems is assumed to be the same. Without the presence of free water the equation of motion for the whole system (Equation 8) and the mass conservation for the total solid (Equation 5) are still solved in four steps and the solution is iterated along with the equation of motion for melt (Equation 7) and mass conservation for melt (Equation 5) to find s i , m i , v s i , v m i and P i (or dP/dz i . The only difference with the DEM model is the new definition of the total solid that is the sum of the solids in the two subsystems s i,(t−1) = s ss1 i,(t−1) + s ss2 i,(t−1) , and similar relations for s i,t and φ s i,t .

Changes to the Temperature Equation
The computation of the temperature using Equation (10) follows closely the procedure already discussed in the previous study. The product C s p s is made of the sum of two terms, C s i,t . The thermal expansion of the whole system is the weighted average of the thermal expansion in the two sub-systems: Thermal conductivity of the two solids is assumed to be the same.

Changes to the Transport Equation for the Composition
The same chemical transport equation used in MF1 (Equation 9) is also applied to the DFM model to determine the composition of the melt and the total solid. Following the procedure for the solution of the mass transport model, the chemical transfer for the total solid is added to the value of sc i,(t−1) at the previous time step and includes two terms, sc i,(t−1) = sc i,(t−1) + ( sc ss1 + sc ss2 )/nthermo. Once sc i,t is computed from the transport equations, the composition of the solid in the two sub-systems can be found considering sc i,t = sc ss1 i,t + sc ss2 i,t and defining a ratio expression similar to the one used to compute the solid volume fraction in ss1 and ss2, that is ℜ = sc ss1 i,(t−1) / sc ss2 i,(t−1) . The ratio is assumed to remain unchanged between two time steps. Combining these two equations, the relation used to compute sc ss2 i,t assumes a simple form: The expression for the ratio is then applied to find sc ss1 i,t . Finally θ sc ss1 i,t and θ sc ss2 i,t are easily retrieved knowing ρ sc ss1 i,t , ρ sc ss2 i,t , φ sc ss1 i,t , and φ sc ss2 i,t .

RESULTS OF THE DYNAMIC FRACTIONAL MELTING MODEL
All the simulations in this and the previous study do not include water which will be considered in the third contribution of the series. Table 1 reports the list of the simulations performed for this study involving the DFM model. Only selected results which subjectively appear to be the most prominent are presented in this section. However the output data may provide additional mean to unfold more findings. Complete data from most of the simulations can be retrieved following the instructions provided in the Supplementary Material. The first case assumes a DFM model with v melt = v solid . Initial and boundary conditions are equivalent to those applied to the dynamic batch melting model discussed in MF1 (simulation PHASE3-P.BRC4-5), additional details of the new simulation can be found in Table 1 (simulation PHASE3B-P.YRC00). The DEM and DFM models can be considered the closest dynamic equivalent to the batch and fractional melting idea commonly referred in classical petrology. Figure 2 summarizes some of the results. In addition the figure also includes the results of the fractional melting model using the pMELTS version of AlphaMELTS (threshold porosity above which melt is extracted set to 0.005).
The difference of the temperature profiles between the DFM model and the fractional melting model with AlphaMELTS (Figure 2A) is mainly due to the contribution of the latent heat of melting, however there is an additional difference between the two models. In the DFM model variable amount of melt is still present in the system at any depth and it affects the thermal profile while with AlphaMELTS only a predefined small mass fraction of melt is always kept with the residual solid. While the net effect of the latent heat on the DFM model in general is difficult to predict, the temperature difference between the DFM model with v melt = v solid and AlphaMELTS should represent the largest possible deviation simply because of the large amount of melt involved compared to DFM simulations with separate melt and solid velocities.
Regardless of the type of dynamic model, the total amount of melt at a particular location is the result of the melt transported (melt moving in/out) and the local melt production from the thermodynamic model. In the DFM model with v melt = v solid the total melt abundance at a given depth is always less than the amount of melt obtained with the dynamic batch melting model ( Figure 2B). It is also not too different from the amount of melt computed by the fractional melting model using AlphaMELTS. The melt volume fraction at a particular depth (or pressure) using AlphaMELTS is given by the sum of the instantaneous volume of melt which is extracted at every point starting from the bottom up to the particular depth.
The larger amount of melt present in the dynamic batch melt model compared with the DFM model with v melt = v solid has only a relatively small effect on the increase of the upwelling velocity ( Figure 2C). The DFM model with v melt = v solid creates melts that are distinctively less rich in silica and with higher TiO 2 content when compared with the melts from the dynamic batch melting model (Figures 2D,E). In the upper portion of the melting column (depth less than ∼40 km) the incompatible component Na 2 O does not show significant differences in the two models ( Figure 2F). The opposite appears to be true for MgO in the same depth range, while at greater depth the MgO content in melt for the two models is very similar (Figure 2G). In the residual solid large differences in the amount of SiO 2 and MgO can be observed at intermediate and low depths (Figures 2D,G). A very small amount of melt crystallizes in sub-system 2 at depths greater than ∼60 km but, the solid appears only in the initial stage of the mantle upwelling and it doesn't form any longer after ∼2 myr (this solid is not reported in the figure, see the data file PHASE3B-P.YRC00.DAT). An animation illustrating the evolution of the DFM model with v melt = v solid over time can be retrieved following the link posted in the Supplementary Material (movie file mf2-movie1.PHASE3B-P.YRC00.avi).     (Figure 4 in MF1, PHASE3-P.YRC4). The conditions imposed to the two models are equivalent, the solid viscosity model is µ s = µ s 0 with µ s 0 = 1e21 Pa s, the permeability constant is C = 1e9 m −2 , the fixed upwelling velocity at the base of the mantle column is −0.03 m/yr (negative upwards), the mantle bulk composition inside the column at the start of the simulation and at the bottom thorough the duration of the simulation are also the same. The only difference between the two models is how melt chemically interacts with the solid. The DFM model has been solved numerically considering the variables defining the pressure gradient dP/dz and the pressure gradient correction dP ′ /dz instead of the pressure P and pressure correction P ′ . The numerical solution obtained using either one of the two solution variables is indistinguishable. The less common application of the pressure gradient may carry some benefits that need to be further investigated (see Appendix 2 in MF1). Additional information about the DFM model can be found in Table 1 (PHASE3B.YRC1). The data file for the DFM model is included in the zip file max-front2-data.zip (see the Supplementary Material for additional information).
Selected properties at two different stages of the melt evolution are shown in Figure 3. The earliest stage refers to the time at which in both models a large amount of melt reaches the top side of the column. For the DEM model it corresponds to the first melt arrival, for the DFM model in general this is not always the case and the large melt occurs at a slightly later time ( Figure 3A). It is more clearly shown in the data file and later in Figure 8. The properties at a later time illustrate the status of the mantle column when the thermal field does not change significantly any further. As discussed in MF1, melt in this particular DEM model never reaches a time invariant state, but for the DFM model instead the melt distribution approaches a condition close to a steady state.
In general the amount of melt at various depths is lower in the DFM model although this is not always the case since the melt distribution observed in the DEM model varies with time.
Velocity of the mantle at the bottom is fixed at −0.03 m/yr (negative upwards) for the DEM and DFM models. The general behavior of the viscous forces is that they decrease with the increase of the melt content (Equation 8), hence within a zeroth order approximation, this observation provides an interpretation for the overall decrease of the upwelling solid velocity in the upper part where more melt is present ( Figure 3B, and also for the next model Figure 4B).
The peculiar feature of the DFM model is the partial crystallization of melt in sub-system 2 (red dotted line in Figure 3A). The amount is relatively small (∼ few % in volume) but most importantly it is a transient feature that disappears with time. By comparing the results from various simulations (see the data files and Table 1), it seems that the duration varies depending on the dynamic and thermal conditions imposed to the mantle column. Figures 3C,D show TiO 2 and MgO in the ss1 and ss2 solid. At the arrival of the large melt volume, the composition of the residual mantle in sub-sys 1 from the two models is quite distinct. In the later stage variation of TiO 2 in the residual solid is not very different in the two models, however MgO shows some significant changes. In the DFM model the composition of the crystallized solid in sub-system 2 varies with depth, in particular it can be noted that TiO 2 is present in the solid only at lower depth. From the data files and following the procedure to determine the bulk composition outlined in the section "Changes in the Thermodynamic Computation", it is possible to apply AlphaMELTS as a post-processing step at any depth and time to re-compute the equilibrium mineralogical assemblage and the minerals composition that were computed but not recorded during the dynamic simulations. Table 2 illustrates the results after 2.35 myr at four depths indicated by horizontal lines in Figure 3A. It is interesting to note that the solid crystallized in sub-sys 2 is made entirely of olivine except at the lowest depth where a small amount of spinel is also present (hence it explains the rise of TiO 2 in ss2). Table 2 also indicates that at this depth the volume fraction of the liquid and the solid in sub-system 2 are comparable and the residual mantle assemblage in ss1 is made of olivine, ortophyroxene and spinel.
The DEM and DFM models create two series of melt that are compositionally quite dissimilar. Figures 3E,F show an example. An interesting point is that the melt composition in the DFM model is very close to the composition observed in the case that assumed v melt = v solid (Figure 2). Inspection of the The first three points include the mineralogy in the two sub-systems. Time = 2.35 myear, depth = 20.2, 35.1, 50.1, 70.2 km, points = 14, 54, 94, 148 (zone 1271). The data file providing the input bulk composition, pressure and temperature is included in the zip file max-front2-data.zip.
Frontiers in Earth Science | www.frontiersin.org data file PHASE3B.YRC1.DAT provides additional insight into the compositional variations of other oxides not specifically discussed here.
An animation of the DFM model shown in Figure 3 can be retrieved following the instruction in the Supplementary Material (mf2-movie2.PHASE3B.YRC1.avi).
The second rheology model considered in the previous study (MF1) assumes that the viscosity of the solid is described by a relation that depends also on the melt fraction µ s = µ s 0 (1 + 1/φ m ). Figure 4 presents the results of the DFM model with this type of viscosity formulation and a comparison with the DEM model results under similar conditions, µ s 0 = 1e20 Pa s and permeability constant C = 1e9 m −2 (DEM model shown also in Figure 5 of MF1, data file PHASE3-P2.YRC13.DAT).
Both DEM and DFM reach a similar steady state in which melt distribution over depth does not significantly change any further over time. In analogy with the observation made for the case that considered the previous viscous model, melt abundance is generally but not always lower in the DFM model compared to the DEM model ( Figure 4A). However with the second rheology model the variation of the melt distribution is related to the presence of only one inversion of the gradient of the melt abundance with depth, while the DEM model develops two inversion points. The nature of the inversion has been discussed to some extent in MF1. It should be only mentioned here that the changed behavior is not the result of turning the petrological condition from melting to crystallization but rather the complex effect of the dynamic transport model. With the DFM model, regardless of the chosen description of the rheology, a second solid crystallizes from the melt in subsystem 2 during the initial stage of the melting process. The melt and solid compositions are similar to those observed with the first viscosity model and they are summarized in Figures 4C-F. The similarity is not entirely unexpected because the net dynamic behavior of the two models is not too different. Even though the value of µ s 0 = 1e20 Pa s is one order of magnitude lower than the value chosen for the previous viscosity model with µ s = µ s 0 , the melt abundance is such that the effective viscosity µ s becomes very similar in the two cases.
The data file for this DFM model (PHASE3B-P2. YRC1.DAT included in the zip file max-front2-data.zip) can be retrieved following the instructions in the Supplementary Material.
In the previous study the temperature dependence was not discussed for the DEM model and only two simulations with fixed temperature at the bottom different than 1,450 • C were included in the zip data files max-front1-data.zip and max-front1-data2.zip (see Table 1 in MF1). In Figure 5 three simulations of the DFM model are presented assuming similar dynamic and chemical conditions but different temperature at the base of the column, 1,450 • C, 1,475 • C, and 1,500 • C. The DFM model with temperature of 1,450 • C is the same shown in Figure 3 (data file PHASE3B.YRC1.DAT). The complete data files of the two simulations with higher temperature are PHASE3B-P.YRC10.DAT (1,475 • C) and PHASE3B-P.YRC12.DAT (1,500 • C).
In general by increasing the temperature at the bottom larger amount of melt is created inside the mantle column ( Figure 5A), but there are additional considerations that can be made. The largest arrival of melt at the top side of the model is not significantly different for the three simulations. However the amount of transient crystallized solid in sub-system 2 increases with the increase of the temperature at the bottom ( Figure 5B). In addition the solid in sub-sys 2 begins to form at greater depth when higher temperature is assumed at the bottom. After a variable amount of time all three models reach a steady melt distribution with only one inversion of the melt gradient with depth. The depth of the inversion point is greater for the model with temperature at the bottom set to 1,500 • C.
The variations of the velocity of the solid reflect indirectly the different temperatures assumed at the base of the model which have a significant control on the overall melt abundance ( Figure 5C).
The compositional difference of the residual mantle (subsystem 1) with temperature is exemplified by TiO 2 and MgO in Figures 5D,E. The behavior does not correlate well, for instance large variations of TiO 2 can be observed as the depth increases (>35 km) while the largest variations of MgO occur between ∼60 and ∼35 km depth. The melt composition does show some dependence on the temperature in particular for the MgO content (Figures 5F,G). The variation is less pronounced for Na 2 O and Cr 2 O 3 and it occurs mostly at depths greater than ∼35-50 km (plots of these oxides is not included in Figure 5). The data files provide a complete overview of the solid and mantle composition.
An animation of the simulation with temperature at the bottom set to 1,500 • C (mf2-movie3.PHASE3B-P.YRC12.avi) is available following the link provided in the Supplementary Material.
A systematic study on the effect of the imposed velocity of the mantle at the boundary (either top or bottom) was also not addressed previously in MF1. Fixing the velocity at the top or bottom boundary may seem a minor detail, after all the total length of the mantle column is only 75 km. From a conceptual point of view the two different conditions describe, at least approximately in 1-D, either passive flow, a regime in which pulling forces (e.g., subducting slab) control the mantle upwelling, or alternatively active flow that could be the result of thermochemical buoyancy (e.g., rising plume).  Table 1.
Perhaps the most interesting observation involves the simulation with the smaller upwelling velocity imposed   at the base (−0.01 m/yr) that somehow stands out compared to the others. The amount of crystallized solid in sub-system 2 is significantly lower (Figure 6A), it also appears at far greater depth (>80 km) and it is persistent for much longer time than for all the other models. The steady state melt distribution includes an inversion point that is much less pronounced than in the other simulations, although it occurs at similar depth (∼35-38 km). The mantle velocities inside the simulation colum are in direct relation with the imposed velocity at the base ( Figure 6B). From Figures 6C-F it is also quite evident that the composition of the melt and residual solid when the largest first melt reaches the top side are remarkably different. The composition of the solid mantle at later time becomes more similar to the composition observed in the other simulations, all being very similar. For the melt instead, the composition remains still noticeably distinct from the other models. The difference also appears for oxides not shown in the figure but accessible through the data file, for example Al 2 O 3 and CaO, in particular at lower depths.
In the previous study (MF1) only one numerical model was carried out assuming a fixed velocity of the mantle on the top side of the mantle column ( Table 1 in MF1) and the results were not discussed. With the DFM model few more simulations with this type of boundary condition have been carried out (Table 1). Figure 7 compares the results of a simulation discussed earlier PHASE3B.YRC1 (Figures 3, 5, 6) that assumes a fixed velocity at the bottom (−0.03 m/yr, µ s = 1e21 Pa s) and two numerical models in which the same constant velocity is imposed on the top side instead. These last two differ in the setup of the solid viscosity parameter µ s that is 1e20 and 1e21 Pa s (data files PHASE3B-P.YRC19.DAT, PHASE3B-P.YRC18.DAT).
The first observation that can be made is that the large first melt arrival on the top side is smaller when the mantle velocity is fixed on this side of the model (Figure 7A, dashed lines). The amount of crystallized solid in sub-system 2 is also smaller with the top velocity fixed, although the depth of the first appearance of the solid in ss2 (∼56 km) is practically the same for the three cases considered (Figure 7B).
The melt distribution at later time is very similar for the model with fixed velocity at the bottom and viscosity= 1e21 Pa s and the simulation with fixed velocity at the top and viscosity = 1e20 Pa s (Figure 7A, green and red solid lines). The other simulation with fixed top velocity but viscosity = 1e21 Pa s shows instead larger melt variation at lower depths ( Figure 7A, blue solid line). Perhaps the most important feature of this simulation is that the melt distribution, in particular in the upper portion of the mantle column is extremely variable over time. This characteristic behavior is not quite evident from Figure 7 but it is clear by inspecting the data file and from the animation that can be retrieved following the link provided in the Supplementary Material (mf2-movie4.PHASE3B-P.YRC18.avi). Figure 7C shows the significant increase of the mantle velocity with depth for the models with fixed velocity on the top side. It should be mentioned that for the model with a transient melt distribution, the velocity of the solid also varies over time inside the mantle column.
When the velocity at the top is kept fixed, the composition of the residual mantle (sub-system 1) appears to be quite different at lower depths ( 40 km) (Figures 7D,E). The model that does not reach a steady state develops temporal and spatial variations of the MgO content at low depths.
In the earlier stage the melt composition shows some variability among the three models (Figures 7F,G). Moving forward in time these differences tend to disappear for two of them, however significant variations of the melt composition are observed for the unsteady model.
The model with fixed velocity on the top side that exhibits large temporal and spatial variations of the dynamic and petrological properties has been tested quite extensively by running several simulations varying the number of spatial grid points, time step, maximum velocity allowed in the model and time interval between thermodynamic computations ( Table 1). All these test simulations seem to confirm the conclusion that, at the conditions imposed on the model, the unsteady behavior is not the result of a spurious numerical artifact.
Because in MF1 the effect of the location of the fixed boundary velocity was not properly discussed for the DEM model, the Appendix presents the results of two new DEM simulations that assume fixed velocity of the solid, one at the top and the other at the bottom side of the mantle column.  (Figures 8C,D). Other simulations instead seem to approach the composition of the DFM model with v m = v s , but only when the model assumes the same fixed temperature at the bottom (1,450 • C). Once the formation of the second solid in sub -system 2 ceased, starting at the point indicated by an arrow in Figure 8D, the composition of the melt products afterwards seems to be affected to some extent. The complete data sets related to these models are PHASE3-P.BRC4-5.DAT (zip file max-front1-data.zip in MF1), PHASE3B-P.YRC00. DAT (DFM with v m = v s ), PHASE3B.YRC1.DAT (green Fewer simulations involving the DFM model have been carried out with the second viscosity model µ s = µ s 0 (1 + 1/φ m ) ( Table 1). Figure 9 includes the results at the point of extraction for the model shown in Figure 4 PHASE3B-P2.YRC1.DAT (purple lines) along with two other simulations not discussed previously, one with the same reference viscosity µ s 0 = 1e20 Pa s but higher velocity at the bottom (−0.09 m/yr) PHASE3-P2.YRC7.DAT (red lines), and the second one with lower viscosity µ s 0 = 1e20 Pa s and fixed velocity on the top side (-0.09 m/yr) PHASE3B-P2.YRC8.DAT (green lines) (all data included in the zip file max-front2-data.zip). Batch melting from MF1 and the DFM with v m = v s are also included for comparison. As expected the models with higher solid velocity, either at the top or the bottom, generate larger melt flux ( Figure 9A). The one with the imposed solid velocity on the top side develops the fastest melt transport (Figure 9B). From a compositional point of view, once the models reach approximately a condition of steady state and the solid in sub-system 2 does not crystallize anymore, they all seem to approach the composition of the melt found by the DFM model with v m = v s (Figures 9C,D).

DISCUSSION AND CONCLUSIONS
The simplest and possibly the most general description of dynamic melting under the condition of chemical disequilibrium between the residual solid and the local melt has been developed for this study. The previous section presented some of the results of this model that has been referred as dynamic fractional melting (DFM). Many of these results will require further analysis, particularly for a comparison with real petrological data. Only a simple observation is made here regarding the crystallization of a second solid which consists essentially of olivine. The results from the DFM model suggest a possible mechanism for the formation of dunite melt channels found in several ophiolites (e.g., Kelemen et al., 1995;Suhr, 1999;Godard et al., 2000;Müntener and Piccardo, 2003;Morgan et al., 2008;Abily and Ceuleneer, 2013).
A comparison with the dynamic equilibrium melting (DEM) introduced in the previous study (MF1) highlighted the main differences between the two models. The effect on the DFM process of the fixed bottom temperature, upwelling mantle velocity at the top or bottom side of the mantle column have been investigated in some detail.
Summarizing some of the findings, the first observation that can be made is that with the DFM model the amount of melt distributed inside the model is often but not always less than the amount obtained with the DEM model. The residual solid and the melt compositions are quite distinctive in the two models. In most cases the unique presence of a second solid in sub-system 2 is a transient feature of the DFM model during the initial stage of the melt evolution. The differences between the DFM and DEM model appear to be persistent regardless of the chosen description of the solid mantle rheology.
In general there is a positive correlation between the temperature at the base of the mantle column and the amount of melt inside the DFM model. The inversion point of the melt gradient, a feature that has been found also in the previous study, appears at greater depth. The composition of the melt and solid are affected by the temperature change but there is no clear predictive pattern to the variations with depth and time.
Considering different values for the fixed upwelling velocity of the solid mantle at the base of the model (active flow), the results obtained with the lowest velocity that was assumed in this study (−0.01 m/yr, negative upwards) somehow deviate from the general behavior observed in DFM models with higher imposed upwelling velocities. In particular this is quite evident for the melt abundance, melt and solid composition in the two sub-systems when the first melt reaches the top side of the model and melt composition at later times.
Although the depth range of the mantle column is rather limited (75 km), the definition of the boundary fixed velocity of the mantle at the top or the bottom side (passive/active flow) has an effect on the DFM process. In particular, under certain conditions, by imposing a fixed velocity on the top side (passive flow), the melt shows large temporal and spatial variations which are not observed when the same velocity is assumed at the bottom side (active flow).
The travel time, flux and composition of the melt at the extraction point (top side of the model) depend on various mantle conditions such as the bottom temperature, imposed upwelling velocity of the solid (either on the top or bottom side), and rheological properties.
Even though the DEM and DFM models allow us to explore the evolution of mantle melting in time and space, they can still be considered simplified representations of the melting process. Some of the limitations have been outlined earlier. The thermodynamic model carries its own limitations which were extensively investigated in a series of studies (Hirschmann et al., 1998(Hirschmann et al., , 1999aAsimow et al., 2001) addressing the earlier version of the melt model (Ghiorso and Sacks, 1995). The newer version used in this study and in MF1 was better calibrated for melting of mantle rocks (pMELTS, Ghiorso et al., 2002), nevertheless some discrepancies with experimental data could not have been avoided (e.g., error on the MgO content in melt is ∼1-4%).
One point that perhaps deserve a further consideration is the composition of the solid mantle source. There is a lot of debate regarding the presence of an eclogitic or component contributing to the mantle melt products (e.g., Ito and Mahoney, 2005a,b;Strake and Bourdon, 2009;Brown and Lesher, 2014;Shorttle et al., 2014;Lambart et al., 2016). It has clear repercussions on our understanding of the petrology of the deep Earth Interior (e.g., Allègre and Turcotte, 1986;Hauri, 1996;Sobolev et al., 2005;Barker et al., 2014;Elkins et al., 2016). From a modeling perspective it won't be difficult to assume a local system at the bottom of the melt column that is made of a weighted average of a peridotitic and eclogitic component. More interesting would be to assume instead that two distinct lithologies coexist in a predefined volumetric proportion, each forming a separate sub-system in chemical equilibrium. This last scheme would be still approachable but one problem that would clearly emerge is how to treat the melt product(s). Once one of the two sub-systems creates some melt or both of them separately create some melt, if some kind of fractional melting similar to the one described in this work is imposed, then some assumptions need to be made on how the melt product(s) will be transferred into a third or fourth sub-system where it/they will mix with all the previous melt products. Perhaps even more challenging is the case of equilibrium melting because, for example, when two melts are present within the same local volume, they may or may not mix and they may or may not equilibrate (together or separately) with one or both lithologies.
Despite the fact that these conceptual ideas can be modeled numerically, it is not known yet whether they can be associated to some realistic melting process or they are nothing more than numerical exercises. There is one additional complication that discourage (at least the author of this study) from pursuing this type of melt modeling for complex petrological assemblages. The characterization of the eclogitic component in the mantle at the present time is based on indirect evidences. There is no understanding of the level of chemical equilibration between eclogite and peridotite in the mantle nor the physical distribution of the two assemblages (Tirone et al., 2015). Only a combination of specific experimental studies to determine the kinetics of the chemical equilibration in combination with large scale geodynamic modeling could provide at least some information on this issue. Unfortunately for various reasons it seems that this type of study will not materialize anytime soon.
Regardless of the above considerations, the DEM and DFM models can be used to try to interpret real petrological data. However for a petrologist using geochemical data a dynamicthermodynamic model may not always be an absolute necessity. A thermodynamic study using for example alphaMELTS should be sufficient to gain some general idea about the temperature/pressure and composition of the source at the beginning of a melting process when it is driven by mantle upwelling. This is quite reasonable in particular if the "true" melting process can be related to the DEM model. In this case often (but not always) the dynamic evolution does not affect significantly the composition of the melt at the top of the mantle column which depends essentially on the initial conditions of the mantle. Any further interpretation that involves temporal and spatial variations of melt abundance and composition would require some model similar to the one in this study provided that the petrological/thermodynamic description is sufficiently accurate. In terms of the procedure that should be followed, there is no easy shortcut or simple formula/tool that can be applied for the task. A possible course of action could be the following. If the petrological data consist only of melt composition on lavas or intrusive components, then plots similar to Figures 8, 9 could be used with the available simulation data to provide a first assessment. Since the exit point is located at 15 km depth (5 kbar), it may be necessary to revise the numerical results or correct the observed melt composition for partial crystallization that may have occurred at lower depths. If samples from the mantle are also available then other plots showing variations of the certain properties with depth may turn to be useful for a direct comparison. When few of the chemical or mineralogical components do not match the real data then it may be necessary to revise the starting mantle composition chosen in this and the previous study, and new numerical simulations would be needed. Petrological data alone may not be sufficient to put a firm constraint on the temporal and dynamic nature of the melting process and more information should be included in the study. Geophysical data, for example seismic, gravimetric, electrical conductivity, could be helpful to constrain the variation of the melt distribution and/or the variation of the thermal state with depth. For this reason it would be easier to apply the DEM or DFM models to present day or recent magmatic events. Geochronological data could also provide support to bracket the temporal window of the melting process and perhaps the volumetric distribution over time.
One important contribution has been left out of the discussion describing the melting process so far, that is the presence of water. It will be the main topic of the third and last installment on modeling dynamic melting in the mantle.

AUTHOR CONTRIBUTIONS
The author confirms being the sole contributor of this work and approved it for publication.