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, 1999a,b; Asimow 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 two-phase flow model (McKenzie, 1984; Scott and 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 kinetic-thermodynamic petrological modeling can be conceived (Tirone et al., 2016), 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”.

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 sub-system 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 melting will occur, the melt is transferred to the second sub-system 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.

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₂ 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₂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₂ 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 3 presents a comparison between two 2-phase flow dynamic melting models. One is a DFM model and the other one is the DEM model discussed in the previous study (Figure 4 in MF1, PHASE3-P.YRC4). The conditions imposed to the two models are equivalent, the solid viscosity model is μs=μ0s with μ0s=1e21 Pa s, the permeability constant is C = 1e9 m⁻², 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₂ 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₂ 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₂ 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₂ 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 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=μ0s(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, μ0s=1e20 Pa s and permeability constant C = 1e9 m⁻² (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 sub-system 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 μ0s=1e20 Pa s is one order of magnitude lower than the value chosen for the previous viscosity model with μs=μ0s, 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.DATincluded 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 (sub-system 1) with temperature is exemplified by TiO₂ and MgO in Figures 5D,E. The behavior does not correlate well, for instance large variations of TiO₂ 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₂O and Cr₂O₃ 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).

Figure 6 summarizes some of the results of four DFM simulations that assume different velocities at the bottom ranging from −0.01 to −0.09 m/yr (negative upwards) while all the other parameters and conditions are kept unchanged. The complete data files are available (see the Supplementary Material), PHASE3B-P.YRC15.DAT (−0.09 m/yr) PHASE3B-P.YRC13.DAT (−0.06 m/yr), PHASE3B.YRC1.DAT (−0.03 m/yr) included in max-front2-data.zip. The data file PHASE3B-P.YRC22.DAT (−0.01 m/yr) is included separately in max-front2-PHASE3B-P.YRC22.zip. Additional information is provided in 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₂O₃ 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.

As in the previous study, the variations of some of the results over time at the point of extraction (exit point or top side of the model) are briefly illustrated here. Selected DFM simulations from those presented earlier assuming the viscosity model μs=μ0s are considered in Figure 8. In addition results of the dynamic batch melting from the previous study and DFM with v^m = v^s are also plotted. Figure 8A illustrates the amount of melt extracted per unit of time. The largest melt flux is observed for the model with the largest velocity of the solid imposed at the bottom (blue line). The earliest arrival of the melt at the exit point is consistent with the large upwelling velocity of the solid and the melt (Figure 8B). The second largest flux (purple line) is the result of the high temperature at the bottom (1,500°C) that generates large amount of melt. Periodic oscillations over time are created by the model with solid velocity fixed at the top side and solid viscosity μ^s = 1e21 Pa s. These temporal variations are not only evident in the melt flux and melt velocity but also in the melt composition (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 line), PHASE3B-P.YRC15.DAT (blue line), PHASE3B-P.YRC12.DAT (purple line) PHASE3B-P.YRC19.DAT (magenta line), PHASE3B-P.YRC18.DAT (black dotted line) (zip file max-front2-data.zip this study).

Fewer simulations involving the DFM model have been carried out with the second viscosity model μs=μ0s(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 μ0s=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 μ0s=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).

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, 1999a,b; Asimow 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 dynamic-thermodynamic 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.

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