Modelling binary alloy solidiﬁcation with Adaptive Mesh Reﬁnement

The solidiﬁcation of a binary alloy results in the formation of a porous mushy layer, within which spontaneous localisation of ﬂuid ﬂow can lead to the emer-gence of features over a range of spatial scales. We describe a ﬁnite volume method for simulating binary alloy solidiﬁcation in two dimensions with local mesh reﬁnement in space and time. The coupled heat, solute, and mass transport is described using an enthalpy method with ﬂow described by a Darcy-Brinkman equation for ﬂow across porous and liquid regions. The resulting equations are solved on a hierarchy of block-structured adaptive grids. A projection method is used to compute the ﬂuid velocity, whilst the viscous and nonlinear di ﬀ usive terms are calculated using a semi-implicit scheme. A series of synchronization steps ensure that the scheme is ﬂux-conservative and correct for errors that arise at the boundaries between di ﬀ erent levels of reﬁnement. We also develop a corresponding method using Darcy’s law for ﬂow in a porous medium / narrow Hele-Shaw cell. We demonstrate the accuracy and ef-ﬁciency of our method using established benchmarks for solidiﬁcation without ﬂow and convection in a ﬁxed porous medium, along with convergence tests for the fully coupled code. Finally, we demonstrate the ability of our method to simulate transient mushy layer growth with narrow liquid channels which evolve over time.


A B S T R A C T
The solidification of a binary alloy results in the formation of a porous mushy layer, within which spontaneous localisation of fluid flow can lead to the emergence of features over a range of spatial scales. We describe a finite volume method for simulating binary alloy solidification in two dimensions with local mesh refinement in space and time. The coupled heat, solute, and mass transport is described using an enthalpy method with flow described by a Darcy-Brinkman equation for flow across porous and liquid regions. The resulting equations are solved on a hierarchy of block-structured adaptive grids. A projection method is used to compute the fluid velocity, whilst the viscous and nonlinear diffusive terms are calculated using a semi-implicit scheme. A series of synchronization steps ensure that the scheme is flux-conservative and correct for errors that arise at the boundaries between different levels of refinement. We also develop a corresponding method using Darcy's law for flow in a porous medium/narrow Hele-Shaw cell. We demonstrate the accuracy and efficiency of our method using established benchmarks for solidification without flow and convection in a fixed porous medium, along with convergence tests for the fully coupled code. Finally, we demonstrate the ability of our method to simulate transient mushy layer growth with narrow liquid channels which evolve over time.
c 2019 Elsevier Inc. All rights reserved.

Introduction
A wide variety of physical processes involve binary alloy solidification. Within industrial settings, understanding this phenomena is relevant to metal casting [1] whilst, in a geophysical context, notable applications include the study of the Earth's core [2] and sea ice formation [3]. Solidification of a multi-component melt from a cold boundary generally leads to the formation of a solid phase with a different composition from the liquid. Heat diffuses away from the solid-liquid interface faster than the solute, causing liquid adjacent to the freezing interface to be constitutionally supercooled [4]. Under these conditions, the solid-liquid interface is unstable to dendritic growth [5] and a mushy layer will form: a porous solid matrix bathed in its melt.
Variations in solute concentration within a mushy layer control variations in density, which can initiate natural convection [6]. Flow through the porous matrix transports heat and solute, promoting local dissolution and internal solidification which modify the porosity. In certain circumstances, dissolution leads to the formation of entirely liquid regions through which the flow is focused. Within sea ice, a binary alloy of salt and water, these regions are known as brine channels and efficiently transport brine into the ocean [3]. In industrial settings, such as metal casting, these 'chimneys' or 'freckles' are undesirable defects [1]. Henceforth, we shall refer to these features as channels.
The solidification of mushy layers and formation of channels has been considered in laboratory studies with a range of binary alloys, the most common being aqueous ammonium chloride. Experiments typically involve either cooling a tank of liquid from either above or below, known as 'fixed chill' [7,8], or translating a Hele-Shaw cell between two fixed heat exchangers; so called 'directional solidification' [9,10,11]. Fixed chill experiments reveal that the location of channels evolves transiently, with the spacing between channels increasing as the mushy layer grows deeper [12,13].
Theoretical descriptions of mushy layers have been developed based on mixture theory [14] and volume averaging techniques [15]. Steady-state solutions can be found in the absence of fluid motion [16], which provide the base state for linear [e.g 17] and weakly nonlinear [e.g 18,19,20] stability analyses. These predict that mushy layers are unstable to small perturbations above a critical Rayleigh number, which characterises the significance of buoyancy forcing relative to thermal and viscous dissipation. Such perturbations result in patterns of varying flow and porosity within the mushy layer, which may give rise to the liquid channels observed in experiments.
However, channel formation is strongly nonlinear; to accurately model the feedbacks involved, the governing equations must be integrated numerically. A split domain approach involves solving separately for liquid, solid, and 'mushy' regions along with the free boundaries which separate them. Various authors have used this method to find steady-state solutions [21,22,23] for a periodic array of channels, given a prescribed channel spacing. An alternative is the single domain 'enthalpy method' [24,25] in which a single set of equations are solved throughout the entire domain. This approach removes the need to explicitly solve for the free boundaries, which are often convoluted. Instead, mush-liquid boundaries are implicitly captured via variations of the porosity field. It is therefore a good choice for computing transient solutions where the channels and mush-liquid boundaries evolve over time. The enthalpy method formulation has been used extensively to simulate binary alloy solidification in a variety of physical settings, using both finite volume methods [10,26] and finite element methods [27,28].
Fully coupled nonlinear simulations are computationally expensive as they feature multiple scales. Narrow channels and diffusive boundary layers develop over time and require high resolution, but are restricted to a small portion of the domain. One approach for efficiently solving problems with multiple spatial scales is the use of adaptive mesh refinement (AMR). Various approaches to AMR exist, which can be broadly divided into either finite element methods on unstructured meshes [e.g. 29], or finite volume methods on structured meshes (following [30,31]). Adaptive mesh refinement has been used to study binary alloy solidification without fluid flow in the context of microstructural dendritic growth [32], and directional solidification of mushy layers [33]. Recently, various authors have studied channel formation and mushy layer convection using a similar model to the one considered here [34,35,36].
We develop an AMR code to simulate the solidification of binary alloys. In contrast to previous studies, which used finite elements, we implement a finite volume method based on block-structured meshes. Supported by the Chombo framework [37], which scales efficiently on thousands of processors [38], we implement adaptivity in time (so called 'subcycling') as well as space: coarser levels can be advanced by fewer, larger, timesteps as each refinement level is advanced at roughly the same Courant-Friedrich-Lewy (CFL) number.
Maintaining conservation laws with subcycling is not trivial. The solution on regions covered by a fine grid is calculated after the coarse grid, meaning that the more accurate fine fluxes cannot directly influence the coarse solution. Previous authors have accounted for this via 'reflux' corrections for hyperbolic problems [31] and elliptic problems with linear operators [39]. For nonlinear elliptic equations, such as those encountered here, Pretorius et. al. [40] applied Berger-Oliger timestepping to the hyperbolic parts of the equations, with the elliptic parts found by extrapolation, before performing an implicit solve for the elliptic parts over the entire hierarchy of grids. We use a different approach; by linearising the elliptic operators, we follow Almgren et. al. [39] and solve implicitly for the small diffusive flux corrections which emerge at coarse-fine boundaries.
In section 2 we introduce the mathematical model which is solved using the numerical method presented in sec-tion 3. We verify the accuracy of the method using benchmark problems for diffusive solidification (section 4.1) and convection in a fixed porous medium (section 4.2). Example solutions of the full model of binary alloy solidification are calculated in section 4.3 using both a uniform mesh and an adaptive mesh, which demonstrate the advantages of using an adaptive mesh for this problem. Finally, in section 5, we summarize our results and suggest future studies which may benefit from the methods presented here.

Problem formulation
We consider a volume-averaged model of a mushy layer [15] in two spatial dimensions (x, z). Whilst it is straightforward to extend the model and our method to three dimensions, we choose to focus on flow and solidification in a 2-D rectangular box of width and height h. We will consider various physical settings within this geometry, which might describe flow in a thin Hele-Shaw cell [10], concentrating primarily on a fluid which is cooled from the top boundary. Solidification increases the solute concentration in the fluid, and therefore the density, which can drive convection.
Adopting the concept of an 'ideal' mushy layer [41], we assume local thermodynamic equilibrium and can therefore determine the local temperature T , liquid fraction (or porosity) χ, liquid composition C l , and solid composition C s from the enthalpyĤ, bulk composition C = χC l + (1 − χ)C s , and an idealized phase diagram ( fig. 1). We denote the liquid (solid) thermal conductivity k l (k s ), specific heat capacity c p,l (c p,s ), species diffusivity D l (D s ), and density ρ l (ρ s ). For simplicity, we assume that solute diffusion is negligible in the solid with D s = 0 . We make the Boussinesq approximation [42], such that the liquid density is approximated as a constant ρ l , except where it affects the buoyancy. We make the additional approximation that the solid phase has density ρ s ≈ ρ l approximately equal to the liquid, with the result that we neglect flow driven by the expansion or contraction of a fluid upon solidification (see [43] for a discussion of this effect). The net effect of these approximations is that the small density differences in the liquid due to compositional and temperature variations control the magnitude of the buoyancy force, whilst the other material properties are approximated by the constant reference density ρ l in all other terms [41].
Given an initial composition C i , we choose characteristic scales for composition ∆C = C e − C i , temperature, ∆T = T L (C i ) − T e , and enthalpy ∆H = c p,l ∆T , where T e and C e are the eutectic temperature and composition, and T L (C i ) is the liquidus temperature at the initial composition (see the phase diagram fig. 1). We scale lengths by h, use the diffusive timescale h 2 /κ l , and scale velocities by κ l /h, where κ l = k l /c p,l is the thermal diffusivity. Finally we scale pressure with ρ l κ 2 l /h 2 and permeability with a reference value K 0 . The dimensionless temperature and composition are defined relative to the eutectic point, whilst the dimensionless enthalpy is given by where the specific heat ratio c p = c p,s /c p,l , whilst the Stefan number S = L/(c p,l ∆T ) gives the ratio between the latent heat release per unit mass of liquid upon freezing, and typical variations in the sensible heat.

Conservation equations
Conservation of energy and solute are given by [10, eqs. 2-4] and can be written in dimensionless form as ∂Θ ∂t where U is the Darcy velocity, Le = κ l /D l is the Lewis number, k = k s /k l is the heat conductivity ratio, and the dimensionless temperature θ(H, Θ) and liquid concentration Θ l (H, Θ) are functions of enthalpy and bulk concentration determined below from the phase diagram. The frame translation velocity V f rame is zero for fixed chill settings, and a constant vector aligned with the direction of sample translation for directional solidification. Recalling that we have made the Boussinesq approximation, and are assuming that ρ s = ρ l , mass conservation is We assume that the fluid density is well described by a linear equation of state, and that momentum conservation can be characterised by a Darcy-Brinkman equation [28] which transitions between the Navier-Stokes equation in pure liquid (χ = 1) and Darcy flow in a porous medium when χ 1. Non-dimensionalising [28, eqn 2.13] yields where the dimensionless permeability Π = Π(χ) is a function of the porosity, and we have defined the Prandtl number, Darcy number, and thermal and solutal Rayleigh numbers, Here, η is the (dynamic) fluid viscosity, α and β are the thermal and solutal expansivities, and g is the acceleration due to gravity. Following [10] we consider a modified Carman-Kozeny permeability function appropriate to solidification in a thin Hele-Shaw cell, In the pure liquid region χ → 1 and the permeability takes the value Π c = d 2 /(12K 0 ), which is the dimensionless permeability for a thin Hele-Shaw cell with narrow gap-width d and sets a limit on the total permeability of the system. In section 3.7 we consider the modifications to the algorithm if the Darcy-Brinkman equation (6) is replaced with the simpler Darcy flow law.

Phase diagram
Depending on the local enthalpy H and bulk composition Θ, a particular region of the domain may be in one of four phases; solid, liquid, eutectic solid, or mushy layer. Before the local temperature, liquid fraction, and liquid composition can be calculated, we must first determine which phase the system is in [24] using a phase diagram (see fig. 1). Following previous authors [e.g. 10] we linearise the liquidus H L (Θ) and solidus H S (Θ) which, along with the enthalpy at the top of the eutectic region H E (Θ), are sufficient to describe the phase diagram: where χ E (Θ) is the porosity on the mush-eutectic phase boundary and the partition coefficient p s = C s /C l is taken as 10 −5 in all our simulations. We have also introduced the so-called concentration ratio, C = (1 − p s )C e /∆C, which describes the size of the greatest solute concentration in the system relative to typical concentration variations. In writing eqs. (9) to (11) we have used the linear liquidus and eq. (1) to note that Θ l = −θ in the mushy layer, Θ = Θ l and χ = 1 along the liquidus, and θ = 0 in the eutectic solid phase. For full details of how we calculate θ, χ, Θ l and Θ s from H and Θ, see Appendix A. , the solid phase of the mushy-layer is composed of the solvent. For example, for an NaCl-H 2 O system with concentration below 230 g/kg the mushy-layer consists of water ice crystals and salty brine. Other alloys, such as NH 4 Cl-water, are typically considered to the right of the eutectic composition (shown with dashed lines). For temperature and composition above the liquidus T L (C) the alloy is purely liquid, and below the solidus T S (C) there is only solid. In between a mushy region forms, comprising a mixture of pure solid crystals on the solidus curve and solute enriched liquid on the liquidus curve, with liquid fraction χ and phase-weighted bulk concentration C. Considering H-C space, with liquidus H L (C) and solidus H S (C), highlights that there is an additional region of eutectic solid in addition to the mush. On the mush-eutectic phase boundary (H E (C) in (b)) T = T e and the porosity depends on the bulk concentration via χ = χ e (C). The eutectic solid is a phase the system often passes through during the transition from the mushy layer phase to the solid phase; the enthalpy difference H E (C) − H S (C) represents the energy which must be extracted for this transition to occur, with T = T e and χ = 0.

Method
We solve the governing equations (3) to (6) on a hierarchy of block-structured, non-uniform meshes using a finite volume discretization, facilitated by the Chombo framework [37]. We extend the AMR algorithm presented in [44] for solving the Navier-Stokes equations to include the Darcy-Brinkman correction, for describing porous media flow, and the effects of buoyancy. Additionally, we extend the Chombo toolbox to solve the nonlinear evolution equations for the enthalpy and bulk composition implicitly using the Full Approximation Scheme (FAS) [see e.g. 45]. Advective terms in the momentum equation and advection-diffusion equations are computed using a conservative unsplit Godunov method [46,47]. The scheme is second-order in space and first-order in time. We introduce our choice of AMR notation in section 3.1. We then introduce our algorithm (section 3.2), which consists of recursive updates on individual levels of refinement (section 3.3) followed by synchronization operations once a timestep has been completed (section 3.4) in order to maintain conservation and stability as needed.

AMR structure and notation
Our notation follows several previous works on AMR algorithms with Berger-Oliger timestepping [31,44,48], which we summarise here for a self contained presentation. A rectangular domain is discretized on a hierarchy of cell-centred grids Γ with grid spacing ∆x , where the level of refinement has index = 0, 1, 2, ..., max . On each level we define the level domain Ω , which consists of one or more rectangular subsets of Γ (see fig. 2a); it is on the level domains that we solve our equations. The coarsest level domain covers the entire problem domain Ω 0 = Γ 0 . Each level domain is properly nested [31]: no interfaces exist between Ω and Ω ±2 , only between Ω and Ω ±1 and posssibly Ω and the domain boundary δΩ.
The grid spacing is chosen such that the refinement ratio n ref = ∆x /∆x +1 satisfies n ref = 2 m , m = 1, 2, 3... . We write the refinement ratios for an entire AMR hierarchy as The timestep on each level ∆t scales with the grid spacing, n ref = ∆t /∆t +1 , and satisfies a Courant-Friedrichs-Lewy (CFL) condition which ensures that ∆t < σ ∆x /Max|U | where the velocity U is taken from the previous timestep and 0 < σ < 1. When solving problems where the fluid acceleration a ∼ Pr Max |U |/Da, Ra T , Ra C is large, but flow has not yet fully developed, we also ensure that ∆t < σ ∆x /a for stability. The level on which a variable is computed is denoted by a superscript , whilst the timestep t n l is indicated by a second superscript, e.g. ψ ,n l +1/2 = ψ(t n l + ∆t n l /2). Unless otherwise stated, time-centred fields are found by linear interpolation e.g. θ ,n+1/2 = (θ ,n + θ ,n+1 )/2. Variables are cell-centred unless otherwise indicated; a subscript f denotes face centred values. Where required, conversion from cell-centred to face-centred values, or vice versa, is achieved by linear interpolation and denoted by Av C→F , e.g.
bottom , require a 'reflux' correction to maintain conservation by ensuring that the fluxes used to update both coarse-and fine-level solutions are consistent.
The level domain Ω can be divided into valid and invalid regions, where Ω valid is the part of Ω not covered by Ω +1 , and we assume that the solution on the finer level Ω +1 is more accurate than the solution on Ω invalid , the part of Ω covered by the finer level mesh Ω +1 . We define single level operators as those which are applied to both valid and invalid regions of using a superscript, e.g. ∇ . In contrast, composite operators are applied to the solution on a union of valid regions and do not have a superscript.

AMR update
Following [44], we implement adaptivity in time via Berger-Oliger timestepping [30,31] as illustrated in fig. 3. To advance the AMR solution from time t to t + ∆t, we employ a recursive timestepping scheme. Starting with the coarsest level = 0, we advance the solution on level from time t to t + ∆t . The next finer level is then advanced n ref times with a timestep ∆t +1 = ∆t /n ref , and so on recursively, where the algorithm used on each level is identical and described in section 3.3. Once all levels finer than have reached t + ∆t , we perform a series of synchronization steps. The solution on a fine level + 1 will have been updated using more accurate fluxes than were used to update , which can lead to a discrepancy between the fluxes across a coarse-fine boundary as shown in fig. 2b, with a resulting loss of conservation. We compute any difference between the coarse fluxes F and the sum of the fine fluxes F +1 , δF = F − F +1 , and correct our solution on the coarse level grid Ω , to account for this as described in section 3.4. This step is called refluxing.
It is desirable for a numerical scheme to be freestream preserving; an initially constant scalar field λ should remain constant when advected by a divergence-free velocity field. However this is not the case in our scheme following refluxing, because the advection velocities are not discretely divergence-free across coarse-fine interfaces. Therefore we use the strategy suggested by [48] and evolve a scalar field λ according to where λ(x, t = 0) = 1. Deviations in λ from its initial value are then used to compute a lagged correction to the velocity field, U p , as part of the synchronisation step (see section 3.4 for more details). The pressure is computed using a projection method, as in [48] section 1.1 but modified to account for the porosity factor in the pressure gradient term. In general, below we solve problems of the form for a pressure p given the initial estimate of the velocity field U which does not, in general, satisfy the divergence constraint (eq. 5).

Single level update
On a single level of refinement , we begin a timestep at time t and wish to evolve the primary cell-centred variables, namely H (t), Θ (t), and the two velocity components U (t) = U (t), V (t) , by a time interval ∆t . In addition to the primary variables we have a scalar λ (t) and a vector U ,n p (t) which we use to correct for errors in freestream preservation at coarse-fine boundaries, as discussed later (section 3.4). The pressure is cell-centred, computed at the half time step t ± 1 2 ∆t , and decomposed into two parts, p = π + e, where π is the pressure computed during a single level update and e is a small correction, computed when synchronizing multiple levels of refinement to ensure that the velocity is divergence-free across coarse-fine boundaries. A separate cell-centred pressure field, φ l , is computed to correct the face centred advection velocities, as described in step 1 below.
Boundary conditions at the interface between a fine level + 1 and coarse level are computed via interpolation in space and time, as described in Appendix A of [44]. Elliptic operators require quadratic spatial interpolation, whilst for advection operators we use bilinear interpolation. In both cases, we use linear interpolation in time following [39,48]. Domain boundary conditions are first order in space and implemented via ghost cells. Because the domain edge is one dimension less than the domain itself, it is possible to drop one order of accuracy along the boundary, using a first order boundary-condition discretization whilst retaining O(∆x 2 ) global accuracy [49,50].
The update has 5 main steps.
Step 1 follows the work of [48], whilst steps 3 and 4 are based on [44].
1. Compute advection velocities. We require a face centred velocity U ,n+1/2 AD , in order to compute the advective terms (U·∇) in eqs. (3), (4) and (6). Rewriting the advection operator as U·∇ (U/χ) = (U/χ)·∇U− U/χ 2 U·∇χ, the momentum equation (eq. 6) becomes Following Appendix A.1 of [44], we treat the right hand side of eq. (14) as a source term evaluated at time t n and predict a face centred velocity U ,n+1/2 f using a hyperbolic tracing scheme. Note that advection on the LHS of eq. (14) used in the hyperbolic tracing scheme is transporting with an effective velocity U/χ. This 'advection' velocity is found by averaging from cell-centres to edges; (U/χ) ,n f = Av C→F (U/χ) ,n . In principle, if U/χ is large enough it may be necessary to use a CFL condition based on (U/χ)∆t/∆x < 1, however this situation was not encountered here. To ensure the advection velocities are divergence free, we solve for the pressure field φ , which is used along with the freestream preservation correction from the previous timestep to construct the final advection velocity: 2. Update scalar fields H, Θ, λ. For each scalar field s ∈ {H, Θ, λ}, upwinded face-centred scalar values s ,n+1/2 are predicted using the hyperbolic tracing scheme described in [48], from which we find the source term S ,n+1/2 s (u) = ∇ · u s ,n+1/2 due to the divergence of advected fluxes. The scalar λ is updated by solving eq. (12) explicitly, We treat the diffusion terms in the enthalpy and bulk composition evolution equations implicitly, where N H and N Θ are nonlinear operators for the diffusive terms.
The coupled set of eqs. (18) and (19) are solved by the nonlinear FAS multigrid algorithm [45], using a Gauss-Seidel routine for both the smoothing steps and for solving on the coarsest level. When evaluating the nonlinear operators (eq. 20), and following the solve, we compute θ, Θ l , Θ s and χ using the phase diagram described by eqs. (9) to (11). 3. Compute new velocity U ,n+1 . We first compute an unprojected velocity U ,n+1 * , treating the viscous and Darcy terms implicitly; where U ,n+1/2 AD,CC = Av F→C U ,n+1/2 AD and (U/χ) ,n+1/2 is found using the same hyperbolic tracing scheme used in step 2 for scalar fields. Following the procedure detailed in step 6 (section 2.2) of [44], we then apply a cell-centred projection operator to U ,n+1 * in order to obtain the approximately divergence-free U ,n+1 [48]; ∇ · χ ,n+1/2 ∇ π ,n+1/2 = 1 ∆t ∇ · U ,n+1 * + ∆t χ ,n+1/2 ∇ π ,n−1/2 , U ,n+1 = U ,n+1 * + ∆t χ ,n+1/2 ∇ π ,n−1/2 − ∇ π ,n+1/2 , where eq. (22) is solved first for π ,n+1/2 to yield the pressure correction in eq. (23). Due to the approximate nature of this projection, we find that it is necessary to manually enforce the physical condition that cells with zero liquid fraction cannot have any fluid flow; U ,n+1 = 0 if χ ,n+1 = 0. 4. Recursively update finer levels. If a finer level +1 exists, it is updated with a fractional timestep ∆t +1 = ∆t /n ref until it reaches the time of the current level, t l + ∆t . Where required, boundary conditions for the solve on level + 1 are computed by interpolation in space and time using the solution just found on level . 5. Synchronize with finer levels. If a finer level + 1 exists, the current level is synchronized with all finer levels as described in section 3.4

Synchronization
Because coarse-and fine-level solutions have been evolved independently, our solution will be continuous but will not preserve properties like conservation, freestream preservation, and the divergence constraint (eq. 5) on the composite multilevel mesh. To ensure that our scheme is conservative, and that the velocity field is divergence free at coarse-fine interfaces, we must perform a series of synchronization steps once all levels finer than some base level base , have reached the same time, t sync .
1. Refluxing. To ensure conservation, we correct for the flux mismatch δF at the interface between a coarse level and the fine level, + 1. Following [48,44], we introduce this mismatch via a source term in our original evolution equations and decompose the flux-corrected solution into the portion already computed and a small correction, e.g. U ,n+1 Corr = U ,n+1 + δU ,n+1 . The discretized equation for U, eq. (21), becomes where F O is the original source term and we have used eq. (23) to eliminate U * . Using eq. (21) we simplify eq. (24) to which is solved implicitly for δU l,n+1 where D R , the "reflux divergence", is the part of the coarse-level divergence operator operating on the flux mismatch on the coarse-fine interface [48]. For enthalpy and bulk composition, we must also solve implicitly for the correction (δH, δΘ) to ensure that our scheme is stable [39,44]. In order to diffuse the correction in a physically appropriate manner, we note that (δH, δΘ) is of the same order of magnitude as the flux mismatches (δF H , δF Θ ), which are themselves O(∆x 2 ) [51]. Hence, we make the approximation χ(H + δH, Θ + δΘ) ≈ χ(H, Θ) and use the linearised forms of the nonlinear diffusion operators defined in eq. (20), in order to compute the composite corrections Finally, the reflux correction for λ can simply be applied explicitly: 2. Apply multilevel projection. To ensure the velocity is divergence free at coarse-fine boundaries, we apply a composite projection as in [48]. A correction to the pressure e s is found by solving across the entire AMR hierarchy, where boundary conditions are as in [44]. The velocity is then corrected: 3. Compute freestream preservation correction. The freestream preservation correction is computed in a composite sense for use in the next timestep. As in [48], we solve for the correction e λ , For stability, the constant satisfies ζ < 1 such that the correction velocity, only partially corrects the solution in a given timestep. Additionally, we find from a linear stability analysis (Appendix B) that the condition ζ < 2(1 − σ) must also be met, in contrast to [44] who instead found that ζσ < 0.5 was a sufficient condition for stability. We typically choose ζ = 0.95, σ = 0.4. 4. Replace coarse grid solution with average of fine grid solution. For consistency, the solution on covered regions is replaced by the average of the solution on the overlying fine region, in a process applied recursively from max down to base .

Initialisation
At the start of a computation, or after the meshes have changed due to regridding, we must initialise the solution. We follow the steps described in [48] (section 3.7) and [44] (section 2.4), repeating the key details here for convenience.
1. At the start of a computation, the primary variables H, Θ and U are given some initial value. An approximation to the lagged pressure π is then found by taking a single non-subcycled timestep across all levels. 2. If regridding during a computation, newly defined fine meshes are filled by linear interpolation from the solution on the underlying coarse mesh. 3. Following any change of the grids, including at the start of a computation, the velocity field is projected in a composite sense to ensure that it is divergence free, and the freestream preservation correction U p is recalculated based on the newly-initialized λ field.

Refinement
Regridding occurs every n r timesteps. During this process, the solution at each grid cell on each level, whether valid or currently covered, is tested against some specified criteria to determine if refinement is required, in which case the cell is tagged for refinement. A new set of grids is then found where these tagged cells are covered by a finer level, whilst still satisfying the rules introduced in section 3.1 regarding proper nesting. This procedure enables the refinement or coarsening of the grids over time.
The appropriate refinement criteria will vary between different problems. Typical choices we make are to tag cells which lie on the mush-liquid phase boundary, where χ < 1 and χ = 1 in adjacent cells, or to tag cells where the undivided gradient of the enthalpy field is above some threshold value. The strategy employed for each problem we consider will be explained where appropriate.

Algorithm for Darcy flow
For flow in a Hele-Shaw cell with a narrow gap width, such that Π c Da 1, the momentum equation eq. (6) is well approximated by Darcy's law where we have rescaled the pressurep = (Da/Pr)p. In this case, the momentum U becomes a diagnostic rather than prognostic variable, and is entirely determined by the instantaneous enthalpy and bulk concentration. Solving eq. (35) instead of eq. (6) can be achieved via the changes to the algorithm described below. The resulting method requires less computational effort, making it advantageous in settings where eq. (35) is a good description of the physics such as flow in a narrow Hele-Shaw cell.
In the single level update (section 3.3) steps 2, 4, and 5 are unchanged, whilst step 3 is no longer required. During step 1, where we compute advection velocities, the predicted face centred velocity is now computed at time t n by averaging the known permeability, temperature, and liquid concentration from cell centres to faces; whilst the pressure correction is found by solving The final advection velocity is given by U AD = U ,n f − Π ,n f ∇ φ + U p . During synchronization (section 3.4) we no longer need to compute a reflux correction for the momentum and the multilevel projection (step 2) is redundant. During initialisation we no longer need to compute an approximation to the lagged pressure in step 1. Following a change of grids (step 3) we do not need to project the velocity field.

Results
We next demonstrate the accuracy and efficiency of our method. First, we will consider two benchmark problems which separately test the two key components of the algorithm: solving for the nonlinear coupling of enthalpy and bulk composition (section 4.1), and the solution of the Darcy-Brinkman momentum equation (section 4.2). Then we consider a series of transient solutions to the fully coupled problem on both uniform and adaptive meshes, which reveal the computational savings provided by the implementation of AMR.

Benchmark 1: diffusive solidification
An analytic solution exists for purely diffusive solidification in a semi-infinite domain [17], which we use to validate our solution to the thermodynamic aspect of the problem. Specifically, we solve our governing equations with no salt diffusion (Le = ∞), equal solid and fluid properties (k = c p = 1), no fluid motion (Ra T = Ra C = 0), and constant vertical frame advection V f rame = V f rame e z with V f rame = 1 on the domain illustrated in fig. 4a. The bottom boundary is held at a fixed temperature above the liquidus θ = 1.1; we then numerically invert the analytic solution to determine the far-field temperature θ(z = −∞) and the corresponding solution within the domain. The top boundary is held at the eutectic temperature and porosity, H top = H(θ e = 0, χ e ), and the x−direction is periodic. Initially, the domain is entirely liquid χ = 1 and at temperature θ = θ analytic . Simulations are stopped once a steady state is reached, defined by max θ ,n+1 − θ ,n < 10 −4 ∆t, for min ≤ ≤ max .
We make use of the adaptive mesh to refine around regions where we expect the error to be largest: the mushliquid interface and wherever the undivided temperature gradient is large (∆x |∇ θ | > 1/2). A typical steady state porosity field is shown in fig. 4b for a simulation with three levels of refinement ( max = 2) and refinement ratio of Defining θ err = θ − θ analytic , in (d) we plot 1D error profiles for three uniform mesh solutions with N z =16 (red), 32 (gold), and 64 (purple) alongside an AMR solution (blue) with a coarse grid N z = 16 and two refined levels. Convergence plots (e, f) confirm that our scheme is second order accurate in space, and that AMR enhances the accuracy of the solution. The refinement ratio between successive levels n ref is either a single value for the case with one refined level (red) or two values for the case of two refined levels (yellow).
2 between each level, which we denote by n ref = (2, 2). A one-dimensional temperature profile is plotted in fig. 4c, along with the error = θ − θ analytic for a selection of AMR and uniform mesh simulations in fig. 4d which reveal that the refined regions greatly reduce the local error compared to the coarser grid. We perform a convergence test using a uniform mesh, an adaptive mesh with one level of refinement (n ref = 2) and an adaptive mesh with two levels of refinement (n ref = (2, 2)), which confirms that the scheme has second order spatial accuracy ( fig. 4(e,f)).

Benchmark 2: convection in a porous medium
To test the coupled evolution of the momentum and heat transfer, we consider convective flow in a porous medium with isothermal heated and cooled sidewalls, and without phase change (S = 0). Following previous studies [28,53,52] we simulate a fluid with equal solid and liquid properties (k = c p = 1), a Prandtl number Pr = 1, and consider a variety of values for the thermal Rayleigh number Ra T and modified Darcy numberDa = Da χ 3 /(1 − χ) 2 whilst Ra C = 0. For the first test case we use a square domain (0 < x < 1, 0 < z < 1), with imposed porosity χ = 0.4 and boundary conditions as illustrated by case b in fig. 5a.
Following [39], we use this problem as a convergence test, and use a fixed variable mesh rather than an adaptive mesh. This choice ensures that all simulations have the same grid hierarchy at all timesteps, permitting an accurate test of the convergence properties. We choose a mesh with enhanced resolution on the left and right walls of the domain to better resolve the boundary layers (see fig. 5b). Once simulations have reached steady state, we compute the Nusselt number, averaged over the left and right boundaries of the domain. As seen in table 1, our results agree well with [28] and [52] over a wide range of parameters. For the cases withDa = 10 −6 , our solutions indicate the presence of a fine scale boundary layer; consequently the quadratic convergence regime is only achieved on the finest grids. To test the capability of the code to handle regions which are entirely solid (χ = 0) and to adaptively refine around regions of interest, we now impose a porosity which varies spatially from 1 at the centre of the domain (x c , y c ) = (0.5, 0.5) to ∼ 10 −45 at the boundaries according to The minimum value of eq. (40) (χ ∼ 10 −45 ) yields a permeability Π = 0 to within machine precision. The resulting convection is confined to the interior of the domain, as illustrated in fig. 5c, and we choose to refine around regions where the dimensionless velocity satisfies |U| > 15. The solution is advanced from an initial state given by θ = x, U = 0 until a dimensionless time t = 2 × 10 −4 , using a fixed timestep ∆t = 1.6 × 10 −4 ∆x. Following [44] we first demonstrate that our scheme is second-order accurate via Richardson error estimation, where the error is estimated as the difference between grids with successive resolutions. Then, we compute the error for both uniform and AMR simulations by considering the difference between them and a uniform 512×512 mesh.

Da
Ra T ∆x coarse [28] [52]   These results are shown in fig. 6a, where the grid spacing on the x−axis indicates the coarsest resolution used in a particular simulation. We see that the error for the 64×64 single level case is similar in magnitude to the 32×32 case with n ref = 2, and the 16×16 case with n ref = 4 or n ref = (2, 2), demonstrating the accuracy of our mesh refinement.
We now illustrate the computational performance of our algorithm by measuring the runtime and total number of grid cells advanced during the 256×256 single-level simulation, and equivalent AMR calculations. As shown in fig. 6b, the use of AMR provides significant savings both in terms of memory (approximately indicated by the number of cells advanced) and runtime.

Fully coupled system
We now apply our method to the fully coupled problem (eqs. (3) to (6)) featuring flow and phase change. First we consider the evolution of a highly porous inclusion within a relatively impermeable material (section 4.3.1), which confirms the accuracy and performance benefits of our scheme. Second, we simulate the solidification of an initially liquid solution from above ("fixed chill") (section 4.3.2), which demonstrates the capability of our method to solve more complex problems.

High porosity inclusion
Here, we simulate the evolution of a high porosity region within a relatively impermeable mushy layer on a periodic domain 0 < x < 1, 0 < y < 1. Such inclusions can act as the precursor to fully developed channels [54], and provide a simple test case for the full coupled system of equations (eqs. (3) to (6)).
We consider a material with differing solid and liquid properties (k = 2, c p = 0.5) which is described by the dimensionless parameters C = 10, Pr = 10, Ra T = 10, Ra C = 10 6 , Da = 10 −3 , S = 5, Le = 100, V f rame = 0 on a doubly periodic domain. An initial state is given by  which generates a small inclusion of higher porosity and higher fluid concentration, with a weak initial flow perturbation. This is evolved in time until t = 1.5 × 10 −4 with a fixed timestep ∆t = 8 × 10 −4 ∆x. In order to replicate a hydrostatic pressure gradient on a periodic domain, we introduce the additional body force where χ b , θ b , and Θ l,b are computed from the background enthalpy and bulk concentration (H b , Θ b ) = (S +0.5, −C /2). We plot the initial and final states in fig. 7. As the dense concentrated fluid sinks under gravity, it increases the local porosity by dissolution of the solid matrix. As in section 4.2, we consider simulations on both a uniform mesh and a variety of adaptive meshes. The refinement criteria is chosen to focus on regions where the undivided porosity gradients are large, i.e. ∆x |∇χ| > 0.1. As seen in fig. 8, we achieve second order spatial accuracy on a uniform mesh, whilst the AMR solution reduces the error to roughly the same magnitude as an equivalent uniform fine-level solution. The use of AMR with a refinement ratio n ref = 4 reduces the memory requirement by a factor of ≈ 10, and the run time by up a factor of ≈ 3 ( fig. 8b).

Fixed chill in a Hele-Shaw cell
To demonstrate the ability of our method to handle more complicated problems, we now simulate the solidification of an initially liquid region from above with V f rame = 0. We consider a unit domain 0 ≤ x ≤ 1, 0 ≤ z ≤ 1, with symmetry boundary conditions on the side walls (i.e. no lateral gradients) and fixed temperatures applied at the no slip top and bottom surfaces which are impermeable to salt according to: The initial state is defined by H(t = 0) = S + 1.3, Θ(t = 0) = −1.0, U(t = 0) = 0, and we solve for the resulting flow and solidification given the parameters Ra C = 10 6 , Ra T = 0, C = 2, Da = 10 −4 , Pr = 10, S = 5, Π c = 10 4 , k = 1, c p = 1.  . "Single-level Richardson" errors are calculated as the difference between computations where the grid resolution differs by a factor of 2. The "Single-level 512" error is the difference between each computation and an equivalent calculation performed on a uniform 512 × 512 mesh. In (b), we compare three simulations with the same finest resolution (∆x = 1/256) and a refinement ratio of 0 (no refinement), 2, or 4. Normalized values are relative to the single-level case, where the number of cells advanced was 3.15 × 10 6 and the CPU time was 180 seconds. For ease of comparison of timings, simulations were performed on a single 2.6 G Hz processor, with GNU compilers and Ubuntu 14.04 operating system.
In general, the choice of mesh refinement criteria will vary depending on the nature of the problem at hand. Here, simulations were run on a base mesh 32 × 32, with 2 levels of refinement using n =0 ref = 4, n =1 ref = 2. The first level of refinement is added wherever the porosity χ < 1, in order to capture dynamics in the mushy layer where the instability begins. The second level of refinement is added wherever the domain has already been refined to level 1 and ∆x|∇χ|(Θ + C )(U ·ẑ) < −0.5, which is chosen to ensure we add further levels of refinement around the channels where porosity gradients are large (hence the ∆x|∇χ| term), bulk solute concentration is high (hence the Θ + C term) and there is significant downwelling (hence the U ·ẑ term). As we are primarily interested in the mushy-layer dynamics we choose not to refine around the plumes once they leave the mushy region, which are left to mix and diffuse as they are advected downward through the liquid. Further levels of refinement could be added using eq. (48) to increase the resolution around the narrow channels. The evolution of the bulk concentration and porosity, along with the AMR grids, is shown in fig. 9. The mushylayer grows diffusively at early times, before the convective instability triggers the formation of channels through which solute is transported into the underlying fluid. The transport of solute out of the mushy layer reduces the local porosity, whilst the salty plumes sink into the liquid. As solidification continues some channels become extinct while others merge with their neighbours, causing the typical channel spacing to increase. The highest level of refinement tracks the channels across the domain, ensuring that they are well resolved with around 10 grid cells across each channel. This provides an efficient means to overcome the limitations on resolution for resolving flow in channels experienced in previous simulations in wide domains with many channels.

Conclusions
We have developed a method for simulating the transport of heat, mass, and solute during binary alloy solidification, using a finite volume approach with adaptive mesh refinement. Our algorithm features local refinement in space and time on block structured Cartesian grids, implemented via the Chombo framework, along with a series of synchronization steps to ensure flux conservation. We have extended the Chombo framework to include a solver for nonlinear elliptic equations using the Full Approximation Storage multigrid scheme. Using a hierarchy of test cases, we have demonstrated the second-order convergence of our method and the performance benefits which AMR can provide. The ability of our scheme to handle transient problems, with features evolving in time, has also been demonstrated.
Future work will look to extend the current algorithm to three dimensions, where we expect the benefits of using AMR to be amplified. In the current two-dimensional setting, this method could be particularly useful for efficiently simulating flows in a Hele-Shaw cell, as demonstrated in section 4.3.2. Here, we expect that our AMR scheme will be most useful for problems with strongly multiscale features, such as solidification of mushy layers with many high porosity channels where the channel width is small relative to their spacing. We anticipate that potential applications may include the studies of sea ice growth, metal casting, and processes in the Earth's interior. The ideas developed here may also be useful in designing algorithms for other related problems, such as flow through different types of porous materials.
where we have defined the CFL number σ = (U∆t)/∆x. Dropping terms of order 2 and dividing through by A n e ikx j leads to Taking the absolute value, we find We require |A| 2 < 1 for the scheme to be stable. By definition, ζ > 0 and σ > 0, whilst −2 < cos (k∆x) − 1 < 0. Therefore |A| 2 will be largest when ζ + σ > 1, such that the final term of eq. (60) is positive. In this case, the most unstable wavelength gives cos (k∆x) = −1 and the criteria for stability becomes Solving the resulting quadratic equation for ζ, and taking the root which gives ζ > 0, we find the stability condition: