Journal of Computational and Applied Mathematics Stability and finite element approximation of phase change models for natural convection in porous media

In this paper we study a phase change problem for non-isothermal incompressible viscous flows. The underlying continuum is modelled as a viscous Newtonian fluid where the change of phase is either encoded in the viscosity itself, or in the Brinkman– Boussinesq approximation where the solidification process influences the drag directly. We address these and other modelling assumptions and their consequences in the simulation of differentially heated cavity flows of diverse type. A second order finite element method for the primal formulation of the problem in terms of velocity, temperature, and pressure is constructed, and we provide conditions for its stability. We finally present several numerical tests in 2D and 3D, corroborating the accuracy of the numerical scheme as well as illustrating key properties of the model. 0 . 01) T on the left wall, no-slip velocities on the particles, and Neumann conditions on the remaining walls. We observe a drag force exerted on the fluid by the array of particles. The velocity is reduced on a very small interval, but flow is allowed through the permeable field. These observations suggest why a viscosity-based formulation would require a much smaller mushy region in order to produce the same melting front as in porosity-based models, and why slow flow is in the solid regime in the porosity For (2.12)


Introduction
The phenomenon of natural convection driven by variations in temperature distribution has been extensively studied from the viewpoint of physical properties and also using computational methods. Common applications include ocean and atmosphere dynamics, design of double glass windows and ventilation devices. If the density of the fluid is approximately constant and the buoyancy contribution depends on temperature, the model equations consist of the so-called Boussinesq approximation [1]. Phenomena that involve phase change in addition to these elements have also a great relevance in many industrial and natural processes, as in e.g. the melting and solidification in the refining of metals.
As the nature of the physical scenario abruptly changes, modelling and computing formalisms usually have difficulty in reproducing the behaviour of the system especially near the liquid-solid interface. Phase changes have been incorporated into the Boussinesq approximation mainly using two different approaches. One is based on enthalpy-porosity models (as in e.g. [2]), where a jump function arising from the so-called Carman-Kozeny equations (see e.g. [3,4]) enforces a large drag force in the solid regime. In other approaches, phase change has been modelled by embedding a jump function into the viscosity, as in e.g. [5]. One objective in the present paper is to give a numerical comparison between these two models. Difficulties in incorporating phase change models are related to the choice of regularisation and jump size parameters. We address this issue with a new viscosity-based model that highlights an appropriate choice of parameters. This model considers the presence of microscopic particles in the solid, which resembles porosity-based models. We choose a transition from fluid to solid having a large gradient, which creates additional numerical challenges.
Recent numerical methods dedicated for phase change Boussinesq models include a class of stabilised discontinuous Galerkin and finite volume methods proposed for porosity-based models in [2] and [6], respectively; and the primal finite element scheme for viscosity-based models, introduced in [5]. However these contributions do not address the stability of the discrete or continuous problems. Theoretical studies are available for other (typically simpler) related stationary models, as natural convection [7,8] including stabilisation analysis and errors estimates, and also timedependent Boussinesq-type problems under different contexts in [9][10][11][12][13][14]. The discretisation of these problems has been associated with Taylor-Hood finite elements for mass and momentum equations and piecewise quadratic approximations for the temperature (as in [7,8,12,15]), or Taylor-Hood and piecewise linear elements as in [5], the MINI-element and Lagrange elements [10] and also piecewise quadratic, piecewise linear and piecewise quadratic for velocity-pressuretemperature as in [9]. Exactly divergence-free methods are available for the stationary Boussinesq equations [16], and other related mixed formulations including a posteriori error estimates can be found in [17][18][19][20][21]. Finite volume, finite difference and Lagrangian schemes have also been used to simulate solidification problems [22][23][24][25]. In summary, a large variety of methods could be employed to solve numerically the equations we look at here. However the stability of the numerical methods applied to this specific case has not yet been addressed, and this is precisely another objective of this paper.
The remainder of this paper is structured as follows. Section 2 recalls the model problems of flow and temperature with and without phase change. In Section 3 we derive a weak formulation of the governing equations and outline a solvability and stability analysis. In Section 4 we introduce two finite element methods based on the primal velocity-pressuretemperature formulation and on the mixed-primal stress-velocity-temperature formulation of the generalised Boussinesq equations. We specify the fully discrete implicit scheme and write down the corresponding Newton linearisation. Next, in Section 5 we present several tests serving as numerical validation for the enthalpy-free case, and we then present a set of comparisons and concluding insights drawn from the simulations of the melting case, collected in Section 6. These tests also include a qualitative analysis on the micro-structure and its relationship with our modelling assumptions, and we close with some remarks and discussions on alternative models.

Main assumptions and model equations
Let t ∈ (0, t f ] denote time and let us consider a homogeneous isotropic porous structure, occupying a spatial domain Ω ⊂ R d , d = 2, 3 and saturated with an incompressible viscous fluid. This fluid has kinematic viscosity ν, thermal expansion coefficient α, and nondimensional specific heat C . The model problem arises from the description of flow using Navier-Stokes and Stefan problems. Applying the so-called Oberbeck-Boussinesq approximation, one ends up with the following set of governing equations written in terms of the velocity u(t ) : Ω → R d , the pressure p(t) : Ω → R, and the temperature θ(t) : Ω → R: ] + ∇p + η(θ)u = f (θ )k,

2)
∂ t θ + u · ∇θ − 1 C Pr div (κ∇θ ) + ∂ t s + u · ∇s = 0, (2.3) and which state the conservation of momentum, mass, and energy with enthalpy, respectively. In (2.1)-(2.3), ε(u) = 1 2 (∇u + ∇u T ) is the strain rate tensor, the function s is the enthalpy, the symbol k stands for the unit vector pointing in the opposite direction to gravity, η, µ are nonlinear functions of temperature that encode the permeability of the porous material and the viscosity of the fluid, respectively. These functions will assume different specifications depending on the phase change model, to be discussed later on. We also define the Reynolds number In order to analyse the coupled system (2.1)-(2.3), we will suppose that the functions µ, η are uniformly bounded and Lipschitz continuous: there exist positive constants µ 1 , µ 2 , η 1 , η 2 , L µ and L η , such that (2.5) Similar assumptions will be placed on the source function f : we suppose that there exists positive constants C f and L f such that On the other hand, we will suppose that for every ψ ∈ H 1 (Ω), we have s(ψ ) ∈ H 1 (Ω), and that there exist positive constants s 1 , s 2 , L s 1 and L s 2 such that for all ψ, ϕ ∈ R. Finally, we suppose that κ is a uniform bounded and uniformly positive definite tensor, meaning that there exist positive constants κ 0 and κ 1 such that conditions are prescribed on the velocity over the whole ∂Ω, and therefore an additional condition is required for pressure uniqueness; as usual we impose a zero-mean property. Regarding the energy equation, we assume that the domain boundary admits a splitting between two disjoint sets Γ θ D and Γ θ N , where temperature and normal heat fluxes are prescribed, respectively. The system is supposed to be initially at rest and isothermal, and so we set u(0) = 0, p(0) = 0 and θ(0) = θ 0 with θ 0 constant.

Enthalpy-porosity models for phase change
The permeability function η appearing in the drag term is usually defined in such a way that (2.1) behaves as the well-known Carman-Kozeny equations (see e.g. the review [26]). That is, one uses a phase change (or liquid fraction) where m > 0 is a small parameter that prevents division by zero. The term δθ represents the temperature range corresponding to the width of the mushy region, and θ f is a constant the jump function is regularised about and it corresponds to the melting point subject to appropriate scaling (in the sense that in the fluid one has φ = 1 by setting η = 0, and in the solid φ = 0 corresponds to η = ξ [m] −1 ). This also implies that in the solid region one imposes a low permeability field that generates a large drag force. The parameter ξ is a large constant that represents the morphology of the melt front. It is related to the imposed permeability through ξ = 180 where d m is the particle diameter, and the constant 180 depends on the material under consideration [27]. Alternatively to (2.9), in our examples we will include a regularised jump function that takes the form where η s corresponds to the relative size of the imposed force and M η is the size of the mushy region. These constants determine the degree of regularisation of the jump. As above, in the liquid phase we have η = 0, and in the solid η = η s , with η s a large constant accounting for the morphology of the melt front.
In many models from the literature, the phase change is often described by combining the permeability (or porosity) regularised jumps with an enthalpy formulation [27]. The non-dimensional enthalpy function s should ideally be a Heaviside function assuming the values s s in the solid, and s l in the liquid. After regularisation using the phase change field φ we employ (2.11) We will adopt this form in all of our models, so that the mushy region for the enthalpy will be predetermined by the temperature range, δθ . The corresponding constants are s l = 0 and s s = [Ste] −1 , where the Stephan number is inversely proportional to the latent heat of fusion L a , and proportional to the temperature range δθ , and specific heat α. More succinctly Ste = αδθ La . For example a material with a larger latent heat of fusion, will have a larger Stephan number embed a larger jump in the enthalpy jump function. Another observation regarding this jump, is that its height depends on the mushy region size. For instance, if δθ decreases then Ste decreases and s s increases. In Section 5 we refer to mushy regions but always in relation with the jump functions imposed on η, µ. In order to avoid cross-effects associated with having effectively two mushy regions (one for viscosity/drag terms and another for enthalpy), we fix in advance a relation that specifies a fixed latent heat.

Enthalpy-viscosity models for phase change
The incorporation of phase change can be alternatively embedded in the form of a temperature-dependent viscosity combined with an enthalpy formulation (as in e.g. [5]). The information about phase variation from solid to fluid is then carried by two different scaled viscosities. This family of models can also be linked to the principle of packing spheres (more naturally associated with the porosity model above, as a dense packing of spheres in the solid translates into a drag force that slows down the flow), and its influence on viscosity variations.
Defining Φ as the ratio of volume occupied by solid particles, we notice that it is related to the phase function φ by Φ = Φ m (1 − φ), where Φ m is a constant depending on the maximum packing of the imposed particles, and so the latter corresponds to the ratio of volume not corresponding to solid particles. As a given solid melts, then there are less solid particles suspended in the liquid, implying an adequate decrease in the viscosity µ = φ −BΦm . Defining n = BΦ m and using a security constant m (mimicking the regularisation in the Carman-Kozeny equation) we obtain µ = 1 φ n + m . (2.12) In analogy to (2.10) we will also employ (2.13) so that in the solid we have µ = µ s , and in the liquid µ = µ l . Here the constant M µ encodes the width of the mushy region.

Relationship with the rheology of suspended particles
Historic models for the rheology of suspensions go back to [28], where the relation µ = 1 + BΦ, with B = 2.5 is postulated. This model can be derived by the consideration of slow flow past a sphere, and it corresponds to linear Newton rheology and is only valid for very dilute suspensions (Φ ≲ 0.01), so it would not be suitable for phase change models.
Extensions to the case of particles with varying size were derived from first principles in [29] and [30]. These models also capture the property that the viscosity tends to infinity as the ratio of solid volume to liquid volume tends to 1, and they deal much better with non Newtonian viscosities observed at higher values of Φ. This form of the viscosity model is still used for nano particles [31], or for the sedimentation-consolidation in macroscopic models [32], and it relates to our model of phase change. Traditionally, these models incorporate differences in density and enthalpy by considering nanoparticles made of different materials. See for instance [33], where thermal and density properties of copper nano-particles are taken into account. In this context, the models in (2.12) and (2.13) consider nano-particles with the same density and thermal properties as the liquid, another important distinction is that the concentration of nano-particles is dependent on temperature in such a way that the particles are not present in the liquid phase, and they exhibit maximum concentration in the solid phase. This is linked to the concept of critical fraction introduced in [30]. There, a packing density becomes high enough that the fluid behaves as a solid, µ = (1 − Φ Φm ) −2.5 where Φ m is a constant depending on the maximum particle packing, and the model was later extended to µ = (1 − Φ Φm ) −BΦm , see [34].

Weak formulation
Firstly let us recall some recurrent notation. For instance, we will write L 2 (Ω) to denote the space of square integrable functions, and will use H 1 (Ω), H 1 (Ω) to refer to the scalar and vector-valued Sobolev spaces W 1,2 (Ω) and W 1,2 (Ω), respectively; whose norms will be denoted as ∥ · ∥ 1,Ω . The inner product in L 2 (Ω) (or in its vectorial and tensorial counterparts) will be simply denoted as (·, ·) and its associated norm as ∥ · ∥. In addition, the space L 2 0 (Ω) denotes the restriction of L 2 (Ω) to functions with zero mean value over Ω. In view of incorporating the boundary conditions for velocity and temperature, we also introduce the space H 1 0 (Ω) of vector functions in H 1 (Ω) whose trace vanishes on ∂Ω, and the space H 1 D (Ω) of scalar functions in H 1 (Ω) whose trace vanishes on the sub-boundary Γ θ D . Associated with the spaces introduced above, the following nonlinear, bilinear and trilinear forms are defined for all u, v, w ∈ H 1 (Ω), p, q ∈ L 2 (Ω), and θ, ψ ∈ H 1 (Ω)  On account of these definitions, we proceed to test (2.1)-(2.2) against adequate functions and integrate by parts conveniently in order to arrive at the following problem in weak form.
such that The forms defined in (3.1) enjoy the following properties, established in e.g. [35] |a θ We can now define the kernel of the bilinear form b(·, ·), characterised by the space of divergence-free velocities V, as Thus, based on (3.3) and the definition of V, the incompressibility condition is included in the functional space and the pressure can be removed from the formulation. That is, problem (3.2) is equivalent to the following problem: For all (a proof can be carried out following e.g. [12]). Moreover, one can see that for all u ∈ V, v ∈ H 1 (Ω), and ϑ, ψ ∈ H 1 (Ω), we have

Stability analysis
Note also that the convective and advective terms can be rewritten using skew-symmetric forms as follows Single phase flows. Let us consider an enthalpy-free counterpart of (3.2), and proceed to derive energy estimates. Testing the energy equation against the temperature solution, and using (3.5) we obtain and then one can use Gronwall's Lemma to assert that Testing now the momentum equation against the velocity solution, exploiting again (3.5), and applying Young's inequality we obtain where we have also used that |k| = 1. We can then we invoke Gronwall's Lemma once again to get Enthalpy-based flows. The following two lemmas will be employed to show the stability and uniqueness analysis of (3.2).
We now establish the stability analysis of problem (3.4).

Proof. It follows from an argument based on Galerkin's method and assumptions (2.4)-(2.7). For more details see [9, Theorem 2.3]. □
The following result establishes the uniqueness of problem (3.2).

Two families of finite element schemes
Let {T h } h>0 be a shape-regular family of partitions of the regionΩ, by triangles (or tetrahedrons in 3D) K of diameter h K , with overall meshsize h := max{h K : K ∈ T h }. In what follows, given an integer k ≥ 1 and a subset S of R d , P k (S) will denote the space of polynomial functions defined locally in S and being of total degree ≤ k.

A conforming method in primal formulation
The spatial discretisation will be based on the finite element method. Accordingly, we define the following finitedimensional spaces for the approximation of velocity, pressure, and temperature respectively: Then the semi-discrete Galerkin method associated with (3.2) reads: A fully discrete method will be obtained after applying the method of lines. Regarding the time discretisation of (4.3), and in view of the overall second order space discretisation expected when choosing k = 1, we here employ a fully implicit second-order backward differentiation formula (BDF2). This choice provides unconditional stability and permits to take sufficiently large timesteps to reach approximate steady state solutions, should they exist. Let 0 = t 0 < t 1 < · · · t N = t f be a uniform partition of the time interval into equi-spaced subintervals of size ∆t, then the method reads: starting from the initial values u 0 h taken as interpolates of the initial data onto V h and Z h , solve for n = 1, . . . the nonlinear and for n = 0 one applies a first order backward Euler method. An advantage of introducing these additional terms is that we have the following relations (for a proof, see e.g. [36]) ), the properties (4.5) immediately hold and the corresponding analysis is substantially simpler.
In much the same way as the continuous case, we define V * the inf-sup condition (4.2), consider the following equivalent problem (see [15], Lemma 3.1): In what follows, we establish the stability result for (4.4), for which the following algebraic identity will be essential: for any real numbers a n+1 , a n and a n−1 , we have where Λa n = a n+1 − 2a n + a n−1 . (4.4). Assume that κ 0 > 2k 1 s 2 . Then, where C 1 and C 2 are constants independent on h and ∆t.  Summing over n and using the estimate (4.10), we finally deduce that Finally, from bounds (4.10) and (4.11) we obtain the results (4.8) and (4.9). □ We now establish the existence of the solution of (4.4).

Theorem 4.2.
Assume that the data satisfy Then, problem (4.6) admits at least a solution (u n+1 Proof. We proceed analogously as in [12,Thm. 3.2]. For notational convenience, we introduce the following constants ) .

A mixed-primal finite element method
Numerical methods based on mixed-primal variational formulations have the advantage that important physical variables (as pseudostress and vorticity) can be approximated directly. Before stating the corresponding Galerkin scheme, we proceed to derive a mixed-primal weak formulation for the model problem. In order to simplify the exposition of this section we will restrict to the stationary counterpart of (2.1)-(2.3), in this case written as follows Also the total stress (or pseudostress tensor, including diffusive, convective, and pressure contributions) is regarded as a new unknown which implies that the second equation in (4.14) together with (4.15) are equivalent to the relations Consequently, we arrive at the following coupled system without pressure t + γ = ∇u in Ω, (4.17) η(θ) u − div σ = f (θ )k in Ω, (4.18) −CPr −1 div(κ∇θ ) + u · ∇θ + u · ∇s(θ ) = 0 in Ω,  ∫ Ω tr(σ + u ⊗ u) = 0. (4.21) Note that the incompressibility constraint is implicitly present in (4.16), and the pressure being in L 2 0 (Ω) is guaranteed by the equivalent statement in (4.21).
A weak form for this problem is obtained after testing (4.16)-(4.17) against suitable functions and imposing the symmetry of σ; multiplying (4.18) by τ ∈ H(div; Ω), integrating by parts and using the velocity boundary condition in (4.20); and taking a similar weak formulation for the energy equation as in (3.2), appropriately modified to incorporate the temperature boundary data in (4.20). That is, the temperature trial and test space will be Moreover, the specific structure of the problem and the orthogonal decomposition H(div; Ω) = H 0 (div; Ω) ⊕ RI, allows us to look for stresses in the space The weak formulation then reads: Find (t, σ, u, γ, θ) ∈ L 2 tr (Ω) × H 0 (div; Ω) × H 1 (Ω) × L 2 skew (Ω) × H 1 0 (Ω) such that A further inspection of this formulation eventually reveals the lack of sufficient regularity (in particular for velocity and stress), which suggests the use of augmentation techniques in the spirit of e.g. [38]. We therefore incorporate residual, Galerkin-type terms in weak form where κ i , i ∈ {1, 2, 3, 4} are positive parameters. Denoting H := L 2 tr (Ω)×H 0 (div; Ω)×H 1 (Ω)×L 2 skew (Ω), ⃗ t := (t, σ, u, γ), and ⃗ s := (s, τ, v, η), we arrive at the following augmented mixed-primal formulation of the initial coupled problem (4.14): where, given an arbitrary (w, φ) ∈ H 1 (Ω) ×H 1 0 (Ω), the forms A φ , B w , a, and the functionals F φ , F D , and G w,φ are defined for all ⃗ t, ⃗ s ∈ H, θ , ψ ∈ H 1 0 (Ω). The solvability of (4.22) can be established by invoking a fixed-point approach between a mixed formulation for the momentum and mass equations and a primal formulation for the energy equation; which can be carried out following [37][38][39]. In contrast with those works, the solvability analysis of the Navier-Stokes equations can be done in our case by applying a further fixed-point iteration that relies on existence results for Brinkman equations, and we require four residual terms and the use of Korn's inequality. Since θ lives in H 1 0 (Ω), we can then apply a similar treatment as in [17,40]. These steps go beyond the scope of the paper and will be addressed in a forthcoming contribution.
The Galerkin scheme and specific finite element subspaces. For each K ∈ T h let us recall the definition of the local Raviart-Thomas space of order k as RT k (K ) := P k (K ) d ⊕ P k (K ) x , where x is a generic vector in R. Then, given k ≥ 0, we The well-posedness of (4.23) will be also based on a suitable adaptation of the continuous analysis of (4.22) to the present context (see, e.g. [37,39]).

Consistent linearisation
Problem (4.4) entails solving a set of nonlinear equations at each timestep. For this we will employ Newton-Raphson's method, which features quadratic convergence provided the initial guess is sufficiently close to the zone of attraction. For a generic nonlinear problem F(w) = 0, one produces a sequence {w k } k converging quadratically to w, through the iterates where DF(v) [δv] denotes the Gâteaux derivative of the functional F along the direction δv. Then at each Newton step k we solve the linear problem where the terms R i h stand for the Newton residuals associated to the momentum, mass, and energy-enthalpy equations. One readily notes that as we increase the Rayleigh number, the coupling between the Navier-Stokes and the temperature equation becomes stronger. This makes the radius of convergence for the Newton method smaller (see e.g. [41]). As the new initial guess in the time-dependent method is the solution at the previous timestep, this condition reflects on a restriction on the timestep, independently of the CFL condition.

Numerical verification
We stress that the zero-mean condition enforcing the uniqueness of the pressure (for the schemes based on the primal finite element formulation from Section 4.1) is implemented using a pressure penalisation approach. All nonlinear systems undergo a Newton linearisation with fixed residual tolerance of 1E−6. In addition, the resulting linear solves are performed with the direct method SuperLU.

Experimental convergence for the semidiscrete and fully discrete methods
For our first example we produce the error history associated with the finite element approximation. Let us consider the following closed-form solutions to the stationary Boussinesq equations with enthalpy, defined on the unit square domain Ω = (0, 1) 2 : u(x, y) = ( sin(π x) 2 sin(π y) 2 cos(π y) , p(x, y) = 10(x 4 − y 4 ), θ(x, y) = 1 + sin(π x) cos(π y). respectively, and the remaining parameters specifying this steady-state version of (2.1)-(2.3) are Re = 10, Ra = 100, Pr = 0.71, γ =1E−6. We then construct a sequence of successively refined meshes for Ω and proceed to compute errors between approximate and exact solutions, together with local convergence rates. The outcome of this error study Table 1 Error history (errors on a sequence of successively refined grids, experimental convergence rates, and Newton iteration count at each refinement level) associated to the spatial discretisation using the finite element spaces (4.1) with k = 1 and k = 2. is depicted in Table 1, which shows optimal convergence O(h k+1 ) for velocity, pressure, and temperature in their natural norms.
The error associated to the time discretisation is assessed by considering the original transient problem, with enthalpy, and employing the following exact solutions (proposed in [12] for the study of the Boussinesq approximation without enthalpy), defined on the three-dimensional domain Ω = (0, 1) 3 : θ(x, y, z, t) = 2 + [x 2 + y 2 + z 2 + 1] sin(t). The regularity of these solutions implies that the spatial finite element discretisation using a method with k = 1 will be machine-precision accurate, and so the total error will practically coincide with the time discretisation error. The temperature-dependent functions are taken as in (5.2). We proceed to discretise the time interval into a sequence of successively refined grids and compute accumulative errors, defined for a generic field scalar or vector field v as The error history is displayed in Table 2, showing a second order convergence, consistent with the BDF2 algorithm employed. We also generate the error history associated with the mixed-primal discretisation discussed in Section 4.2. The closedform solutions are those in (5.1) and the model parameters and temperature-dependent functions are as in the first example above; while the stabilisation constants needed for the augmented formulation take the values κ 1 = κ 2 = µ 1 µ −2 2 , κ 3 = µ 1 /2, κ 4 = µ 1 /4 (and where the constants µ 1 , µ 2 are the bounds for the viscosity introduced in (2.5)). The problem in this case is solved in terms of strain rate, total stress (including viscous and convective contributions), velocity, vorticity tensor, and temperature; and each individual error is measured in its natural norm. The error decay and the corresponding convergence rates are reported in Table 3, which confirms an optimal convergence rate. The Newton iteration count is similar to the one observed in Table 1.

Benchmark test: natural convection of air
We further validate the numerical method against a well-documented benchmark, the natural convection of air in a differentially heated square cavity Ω = (0, 1) 2 . In this problem we do not have the enthalpy terms, we do not consider Table 3 Error history (errors on a sequence of successively refined grids, experimental convergence rates, and Newton iteration count at each refinement level) associated to the spatial discretisation using the mixed-primal formulation from Section 4.2, using a lowest-order scheme with k = 0.  The simulation is run until the final, dimensionless time t f = 0.5, using a Taylor-Hood finite element family for the approximation of velocity and pressure (i.e., k = 1), and the pressure penalty parameter takes the value γ =1E−7. The flow is driven by the difference of temperature and examples of velocity, pressure and temperature distribution at the final time are depicted in Fig. 5.1(a-c). We observe well-defined temperature profiles and the expected recirculation velocity patterns. A more quantitative study is done by extracting the approximate solutions for temperature and vertical velocity on the horizontal midlines at y = 0.5 and plotting them against properly rotated published benchmark values from [42] (which were generated using the method of discrete singular convolution). A reasonably close match is confirmed by looking at Fig. 5.1(d,e), where we emphasise that our results come from coarse mesh computations.
We carry out further comparisons based on the average Nusselt number on the hot wall of the cavity, that is, at x = 0.
The value is here defined as where M denotes the hot wall, and it encodes the rate of heat transfer along M (including the total flux, even the part coming from advection). We also record the maximum and minimum velocities and temperatures attained on the symmetry lines x = 0.5 and y = 0.5. The computed values are collected in Table 4, where we also include reference values from the literature (see also [25,44]).

Simulating the melting of N-octadecane
We now consider the melting of a solid phase change material (N-octadecane) contained in a square box and subject to heating on the left side of the domain. The problem set up is taken from [5], where the boundary conditions are as above (no-slip velocities on the whole boundary, the temperature has zero flux on the top and bottom walls, and high and low temperatures are maintained on the left and right walls, respectively). The low temperature imposed on the right wall θ C = −0.01 is lower than the phase change temperature θ = 0, in order to allow the phase change to occur. We here employ a structured mesh of 100 000 elements.
We first consider an enthalpy-viscosity model. The model parameters are as follows Ra = 3.27E5, Pr= 56.2, Re= 1, κ = 1, Ste= 0.045, high temperature on the left wall θ H = 1, and size of the mushy region δθ = 0.1. We then use (2.11) with s l = 0, s s = [Ste] −1 , and we also employ the regularised viscosity specification (2.13) with the constants µ l = 1, µ s = 10 8 , M µ = 50, and θ f = 0. These values give a constitutive relationship for the phase change that can also be recovered from (2.9) using the phase change function φ together with the viscosity µ(φ) = µ s + (µ l − s s )φ. In that case, the jump is sufficiently large so that the mushy region acts as a solid near the melting point. Notice that in [5] the convection of enthalpy (the last term in the LHS of (2.3)) is neglected. This term is zero almost everywhere, except within the mushy region. So if this region is relatively large, then the term will have an important effect in the generated flow patterns.
In Fig. 6.1 we present snapshots of the numerical solutions obtained using an enthalpy-viscosity model (and setting constant drag), and an enthalpy-porosity model (using constant viscosity). After an initial stage where the dynamics of the system are dictated predominantly by heat conduction, the convective effects in the mixture start to dominate and the layer that determines the phase change moves on to the right on the top of the cavity. We can observe that even if the jump definitions are qualitatively similar, substantial differences appear in terms of the pressure profiles in the solid region. The velocity patterns are also different (the porosity-based model permits recirculation even in the solid), but this can be straightforwardly remediated by taking a larger value for η s . Secondly, we also observe that the phase change boundary is more advanced using the enthalpy-porosity model, and this is more pronounced at the top of the container. In the next set of tests we will address these differences in further detail.

Changing the size of the mushy region and the jump nonlinearity
In view to investigate the differences between porosity and viscosity based models, we recall that enthalpy-viscosity models can be very sensitive to dynamic viscosity effects, thermal conductivity variations, and the size of the mushy region (see [45]). On the other hand, the enthalpy-porosity model for the melting of gallium, and studied in [2] is sufficiently accurate for a specific mushy region size. We then proceed to vary the size of the mushy region in an adequate manner to conciliate enthalpy-viscosity and enthalpy-porosity models. It is of note, we do not vary the enthalpy from that of the last experiment. This is because the enthalpy-porosity formulation allows the regularisation of the jump for the drag force η(θ) independent to the regularisation of the jump function for enthalpy, i.e. by changing the value for n. So for consistency we vary M µ without varying s(θ ).
We start with different values for M µ in the enthalpy-viscosity model, and we observe that with large mushy regions one can mimic the results obtained using enthalpy-porosity models. However, smaller mushy regions in enthalpyviscosity models do not necessarily imply a higher irregularity in the model coefficients, as their counterpart in enthalpy-porosity models need to impose a larger jump in order to prevent unwanted flow in the solid region. The results of mushy region variations in enthalpy-viscosity models are collected in Fig. 6.2(a).
A second investigation is done by modifying the degree of nonlinearity in the Brinkman term. Setting the parameters ξ = 10 7.4 , m = 10 −2.6 in the specification of the enthalpy-porosity model (2.9) imply that the permeability η is regularised over the desired mushy region, with a constant larger than 10 8 reducing the flow in the solid region. We recall that ξ is related to the material morphology, and in the context of the study of melting N-octadecane, one could choose a more appropriate value. Also, the phase change function φ could be regularised around some artificial melting temperature with a different choice of mushy region δθ . A closer inspection of the model coefficients reveals that the present form is equivalent to (2.10) with a mushy region of size M η = 150 (which is a large value, especially considering that enthalpy-porosity models require less regularisation than enthalpy-viscosity models). Making the link with variable viscosity we recall that in the model by Brinkman one has (2.12). We use the value n = 2.5, and choose m = 10 −8 to produce a jump of order 10 8 . Plotting then η it is clear that the phase change function φ requires a different parametrisation in order to get the jump over the required mushy region. We then choose φ = 1 2 [tanh(50(x − 0.074)) + 1], and stress that this model is based on the theory of suspended particles and how the packing density affects viscosity. Again, using (2.12) and varying n, one can mimic the effects of (2.9). Comparisons of the melting fronts produced with these two last family of models are displayed in Fig. 6.2(b), whereas a sample of the numerical solutions choosing n = 5 (and having the following regularisation of the phase change function φ = 1 2 [tanh(50(θ − 0.0366)) + 1]), and a locally refined mesh are portrayed in Fig. 6.3.

Flow patterns in a local element
We next turn to the simulation of buoyancy-driven flow within a mesoscopic subdomain (or local element) lying on the boundary layer between liquid and solid. In line with the observations made above, the solid phase can be regarded as a porous medium filled with solid particles. It is expected that these obstacles will generate a drag as the one encoded in the macroscopic form (2.9). In a local element of diameter 0.01, we create a so-called blockage arrangement of randomly distributed particles of varying sizes with mean d m = 0.0002 and a 1% variance. The boundary conditions for the flow are different than in the macroscopic case. We impose a slip velocity u = (0, 0.01) T on the left wall, no-slip velocities on the particles, and Neumann conditions on the remaining walls. We observe a drag force exerted on the fluid by the array of particles. The velocity is reduced on a very small interval, but flow is allowed through the permeable field. These observations suggest why a viscosity-based formulation would require a much smaller mushy region in order to produce the same melting front as in porosity-based models, and also why slow flow is allowed in the solid regime in the porosity formulation. For instance, for the form (2.12) one should consider nano particles that melt at the melting temperature of the material, and have the same properties as the liquid, therefore their only effect on the fluid would be an increased viscosity. This behaviour is achieved simply by not being fluid and by colliding with each other. This phenomenon amounts to impose a moving permeable field, eventually leading to the formation of a larger mushy region. In summary, in enthalpy-porosity models the particles are fixed and the drag force hinders the flow, whereas for enthalpy-viscosity models, the holes are allowed to move and the amount of fluid present is reduced. An example is given in Fig. 6.4, and a somewhat similar study can be found in [46].

Concluding remarks.
We have addressed the modelling of phase change in Boussinesq models within porous media. A finite element method has been proposed for its numerical approximation, and we have established stability of the continuous and discrete equations. We have tested the performance of the method using a classical benchmark for air convection, where the scaled viscosity is one, there is no porosity, and no enthalpy terms. Secondly, we have simulated the melting of a material, where the phase change is incorporated in two alternative ways: either using viscosity or porosity as main effects producing the interface movement. Extensions of this work include the derivation of error estimates, the  generalisation of the enthalpy-viscosity and enthalpy-porosity models to include nonlocal contributions and performing further comparisons with alternative models, such as those based on the error function and reviewed in [47,48]. We would also like to extend the applicability of the methods proposed here, in the study of large scale models including the thermal evolution of magma/ocean interfaces [49], or ice-shelf melting [50]. A preliminary test is given in Fig. 6.5 where we have implemented an enthalpy-viscosity model to simulate the ice-shelf melting around Antarctica, using an unstructured mesh consisting of 20500 elements and the transient version of the mixed-primal finite element scheme described in Section 4.2. For a relatively large mushy region we can observe a stabilisation of the recirculation patterns and a concentration of the temperature gradients towards the phase change layer.