A SPACE-TIME DISCONTINUOUS GALERKIN SPECTRAL ELEMENT METHOD FOR THE STEFAN PROBLEM

. A novel space-time discontinuous Galerkin (DG) spectral element method is presented to solve the one dimensional Stefan problem in an Euler- ian coordinate system. This method employs the level set procedure to describe the time-evolving interface. To deal with the prior unknown interface, a backward transformation and a forward transformation are introduced in the space-time mesh. By combining an Eulerian description, i.e., a ﬁxed frame of reference, with a Lagrangian description, i.e., a moving frame of reference, the issue of dealing with implicitly deﬁned arbitrary shaped space-time elements is avoided. The backward transformation maps the unknown time-varying interface in the ﬁxed frame of reference to a known stationary interface in the moving frame of reference. In the moving frame of reference, the transformed governing equations, written in the space-time framework, are discretized by a DG spectral element method in each space-time slab. The forward transformation is used to update the level set function and then to project the solution in each phase back from the moving frame of reference to the ﬁxed Eulerian grid. Two options for calculating the interface velocity are presented, and both options exhibit spectral accuracy. Benchmark tests indicate that the method converges with spectral accuracy in both space and time for the temperature distribution and the interface velocity. A Picard iteration algorithm is introduced in order to solve the nonlinear algebraic system of equations and it is found that just a few iterations lead to convergence.

1. Introduction. The Stefan problem is a moving boundary problem that models phase change; e.g., the freezing and thawing process for a solid-liquid system [49,38,43,7]. A key issue in the Stefan problem is that the interface between different phases is evolving in time. The heat equation for each phase is solved in a timedependent domain. The velocity of the phase boundaries in turn depends on the jump in the normal heat flux at the interface separating the two phases. In other words, the problem requires one to find the solutions in a prior unknown domain and to compute the shape of the unknown domain as part of the solution, which makes the problem nonlinear. Such a problem arises in numerous industrial and technological applications, for example, metal solidification, food freezing, dendritic 2. Mathematical model. The Stefan problem describes the temperature distribution of a pure material undergoing a phase transition. Let us consider a square domain, Ω ∈ R d , of a pure material but at different states, such as ice and water. Let x = (x 1 , x 2 , ..., x d ) be the spatial variable, where d is the spatial dimension. The domain Ω is composed of a time-dependent solid region denoted as Ω s t , and a time-dependent liquid region Ω l t . The interface Γ(t), separating different phases, is a prior unknown free boundary. The description of the Stefan problem on a Cartesian domain, Ω ∈ R 2 , is illustrated in Fig. 1. 2.1. Governing equations. We assume that convection, thermal expansion and buoyancy are not considered. Since there is no flow in the liquid region, the governing equation in each phase is written as follows, and Ω l t are time-dependent domains at time t, while domain Ω is fixed. Γ(t) is the moving interface at time t, and n = (n 1 , n 2 ) is the interface normal pointing in the direction of interface movement. φ is the level set function introduced in section 2.2.
where ρ is the density, c v denotes the specific heat capacity, k is the thermal conductivity, θ is the temperature distribution, and indices "s", "l" denote the solid and liquid phase, respectively. The jump condition, also referred to as the Stefan condition, at the interface is (with the assumption of ρ s = ρ l = ρ and c s v = c l v = c v ), where L is a positive constant representing the latent heat of fusion, n = (n 1 , n 2 , ..., n d ) is the interface normal pointing in the direction of interface movement, and V n = V · n is the normal velocity at the interface. On ∂Ω, Dirichlet boundary conditions are specified. At the interface, we set For a classical Stefan problem, θ Γ = θ m where θ m denotes the melting temperature. However, for a realistic modeling, such as dendritic solidification, it is important to take surface tension effects and interface curvature into account, which leads to the Gibbs-Thomson relation [11]. We refer the reader to [24,23,58,57,17] and the references therein for information on freezing models and numerical methods that take into account more detailed physics than we do in this article.
A SPACE-TIME DGSEM FOR THE STEFAN PROBLEM 3599 The evolution of the level set function is determined by the following equation, where W is a continuous extension of the normal velocity V n of the interface, such that, where x * is the closest point on Γ(t) to the coordinates of a point (x, t). The outward unit normal pointing in the direction of interface movement, n, can be determined by the level set function φ as follows, The level function φ is initialized as a signed distance function where d Γ is the distance of a point x to the interface Γ. Updating the level set function by Eq. (5), may change the level set function away from being a distance function. Thus, it is necessary to perform a reinitialization procedure to keep the level set function as a distance function, which satisfies |∇φ| = 1 (cf. [56]).
3. Space-time discretization. In this section, we follow the notation of, e.g., [53,46], to introduce the definitions of the space-time domain, space-time slabs and space-time elements. In a space-time discretization, we introduce a space-time domain E as E = Ω × [t 0 , T ] ∈ R d+1 . The coordinates of a pointx ∈ E are defined asx = (x, x d+1 ) with the spatial variables x = (x 1 , x 2 , ..., x d ) and the time variable t = x d+1 . The space-time domain boundary ∂E consists of Ω 0 := {x ∈ ∂E | x d+1 = t 0 }, Ω T := {x ∈ ∂E | x d+1 = T }, and Q := {x ∈ ∂E | t 0 < x d+1 < T }. First, we partition the time interval [t 0 , T ] by the time levels 0 = t 0 < t 1 < ... < t E (t) = T . The space-time domain E is then divided into E (t) space-time slabs. The n-th space-time slab is denoted as E n = E ∩ I n , where I n = [t n , t n+1 ] is the n-th time interval with length ∆t n = t n+1 − t n . Now, we describe the construction of the space-time elements in the space-time slab E n . Divide the spatial domain Ω tn and Ω tn+1 into E (x) nonoverlapping spatial elements K n and K n+1 with a uniformed size h. A space-time element K n is then obtained by connecting K n and K n+1 . The element boundary ∂K contains K n , K n+1 and Q n K . On ∂K , the outward unit space-time normal is denoted byn K = (n, n d+1 ).
The Stefan problem on a Cartesian domain Ω involves two different phases, i.e., Ω s t and Ω l t . Thus, the space-time elements are classified into two categories: regular cells and cut-cells.
• Regular cells. A space-time element is referred to as a regular cell if it contains only one phase. We denote a regular cell as K n,s j if the space-time element K n j is in the solid phase. Similarly, K n,l j denotes a regular cell which contains only liquid phase.
• Cut-cells. A cell is called a cut-cell if a regular cell intersects with the moving front. In particular, a cut-cell contains different phases, which are separated by the phase boundary. Here, we denote a cut-cell as K n,sl j if the space-time element K n j intersects with the front Γ. The subelements associated with the cut-cell K n,sl j are denoted as K n,s,Γ j and K n,l,Γ j .
To generate a cut-cell in a space-time slab E n , one can partition the time interval [t n , t n+1 ] by a set of Gauss-Lobatto points (τ 0 , τ 1 , ..., τ p ), and then determine the zero level set crossing values Γ(τ i ) (i = 0, ..., p). By connecting Γ(τ i ) using high order interpolation in time, we introduce a curve Q n Γ cutting a regular cell into two subelements, i.e., K n,s,Γ j and K n,l,Γ j . (See Fig. 2). The space-time slab discretization is a tessellation of a domain E n , which implies a division into non-overlapping, space-time-filling subsets. In the space-time slab E n , the tessellation of the solid (liquid) phase is defined as T n,s h (T n,l h ), which consists of all regular cells, K n,s (K n,l ), and subelements of cut-cells, K n,s,Γ (K n,l,Γ ). Fig. 2 illustrates a tessellation of the space-time slab E n in a space-time domain, which is composed of tiles [16]. As we see in Fig. 2, regular cell, K n,s j , is one tile, while cut cells, K n,s,Γ and K n,l,Γ , make up two tiles. Then the tessellation of the space-time slab E n is defined as T n h , such that T n h = {T n,s h , T n,l h }. In addition, the tessellation of the space-time domain E is denoted as T h = ∪ n T n h . In a space-time domain E ∈ R 2 , a sketch of a space-time slab E n with the moving interface Γ(t) is demonstrated in Fig. 2. 4. Fixed frame of reference versus moving frame of reference. In a spacetime slab E n , a cut-cell K n,sl j is generated due to the moving front Q n Γ intersecting or embedding within a space-time element. Some examples of a cut-cell are shown in Fig. 3. A subelement generated in a cut-cell could be a triangle, a quadrilateral or a pentagon. Note that for the Stefan problem, a key issue is that the phase boundary is transported with a velocity that is proportional to the jump of the normal heat flux at the evolving and prior unknown boundary. Instead of dealing with implicitly defined arbitrary shaped space-time elements, we transform a prior unknown time-varying front into a known stationary front by applying the idea of Arbitrary Lagrangian-Eulerian Methods [37,41]. First, two coordinate systems are introduced.
• Fixed frame of reference. The fixed frame of reference is the frame that we use to define the physical problem. The coordinates of the fixed frame of reference are denoted as (x, t), where x = (x 1 , x 2 , ..., x d ) is the spatial variables. (See Fig. 4 (a, c)). • Moving frame of reference. The moving frame of reference is the coordinate system that we use to define the transformed physical problem. In the moving frame of reference, the coordinates are defined as (x, t). (See Fig. 4 (b, d)).
Next, we construct a mappingx = M b (x, t) between the moving frame of reference and the fixed frame of reference. In particular, the mapping will take a point on the time-varying front in the fixed frame of reference and transform it to a point on the stationary front relative to the moving frame of reference. For example, Γ(t n ) and Γ(t n+1 ) in Fig. 4 (a or c) denote the zero level set at time t n and t n+1 , respectively. The transformation M b (x, t) transforms Γ(t n ) and Γ(t n+1 ) toΓ(t n ) andΓ(t n+1 ) in Fig. 4 (b or d), such thatΓ(t n ) = Γ(t n ) andΓ(t n+1 ) = Γ(t n ). Thus, , where x = x 1 is the spatial variable and x 2 denotes the time variable. The square with bold lines is an example of a space-time slab E n , such that E n = E ∩ I n , where I n = [t n , t n+1 ]. We assume that the solid phase is on the left hand side of the interface Q n Γ , while the liquid is on the other side. The space-time element K n,s j is referred to as a regular cell in the solid phase. However, K n,sl j+1 and K n,sl j+2 are examples of a cut-cell due to the time evolving interface, Q n Γ , defined on a nonconforming computational grid. In this example, the curve Q n Γ obtained by 4th order interpolation cuts K n,sl j+1 into two subelements, i.e., a pentagon K n,s,Γ j+1 and a triangle K n,l,Γ j+1 in a space-time slab E n , our transformationx = M b (x, t) is determined by where ∇ denotes the gradient operator with respect to x, and W is defined in Eq. (6). Since the zero level set of the transformed front is determined by the front at time t n , we refer to this transformation as the backward mapping transformation. Note that matrix ∇M b is the Jacobian of the transformation, which is defined by The moving front is evolving from left to right in each case. The shape of the generated subelements could be triangle, quadrilateral or pentagon.
After updating the temperature distribution in the moving frame of reference and calculating the updated interface curve Q n Γ , we need to project the solutions back onto the fixed frame of reference. This projection x = M f (x, t) in a space-time slab E n is called the forward mapping transformation and is defined by where∇ denotes the gradient operator with respect tox. In one spatial dimension, we apply the method of characteristics to solve the transformations (8) and (10), and obtain, • the backward transformation, • the forward transformation, x(x, t) =x + t tn W (s)ds.
By applying a mapping technique, the phase front becomes a straight line in the time direction in the moving frame of reference (see Fig. 4), which reduces the cost of quadrature on the implicitly defined surfaces and volumes. In one spatial dimension, the tessellation of the space-time domain becomes rectangular meshes in the moving frame of reference. In each case, the curve Q n Γ is mapped to a straight bold lineQ n Γ , which is either aligned with an inter-element boundary or cutting a regular cell into two subelements.
4.1. The transformed equations in the moving frame of reference. In this section, we describe a general approach to transform the heat equations in the fixed frame of reference, Eqs. (1) and (2), into the moving frame of reference. In the moving frame of reference, the coordinates of a point is denoted as (x, t). As the mapping procedure is identical in both phases, we consider only the solid phase, Eq. (1), and drop the index "s", Given a mappingx = M b (x, t), we find the effect of the mapping on Eq. (13) by use of the chain rule. First, the time derivative of θ transforms as, whereθ denotes the temperature distribution in the moving frame of reference. For the gradient operator ∇θ, we have where superscript " " denotes the matrix transpose, and ∇x is the Jacobian of the transformation, which is defined in Eq. (9). Similarly, we havẽ It then follows that From (15), (16) and (17), we obtain The transformation for the divergence operator, ∇ · F , is defined as where the symbol ":" is defined as A : By applying the transformations Equations (14) and (18) to (20), the transformed heat equation is then written as In one spatial dimension, the transformation (10) is linear due to the velocity extension, W , is a function of time only. We then have∇W = 0 and ∇x = (∇x) −1 = 1 based on the transformations (11) and (12). The transformed heat equation is then simplified as Hence, the space-time formulation of Eq. (22) is written as, where ∇ is the space-time gradient operator with respect tox = (x, t), B ∈ R 1+1 , and A ∈ R (1+1)×(1+1) ; B and A are defined as, 4.1.1. Proposed extension of space-time mapping algorithm to multiple dimensions.
In multiple dimensions, one would solve (10) by the method of characteristics where a Picard iteration technique is applied to find the arrival points of the characteristics. At each iteration, one solves with the initial guess x (0) (x, t) =x, which is the initial condition of Eq. (10). It is possible that the transformations become singular in multiple dimensions. It would then be necessary to reduce the spatial and temporal order as well as reduce ∆t n if the condition number of∇x exceeds a threshold. One can also use the condition number of∇x as an indicator for adaptive mesh refinement.
Note that using vector identities, the transformed heat equation (21) in multiple dimensions is written as where Lθ is a non-conservative term defined as follows The term Lθ is linearized so that the nonlinear parts are frozen for each separate Picard iterate.

5.1.
Outline of the space-time DG spectral element method. We present an outline of our space-time DG spectral element method for solving the Stefan problem in Algorithm 1. Note that the transformed heat equation (22) in the moving frame of reference is a nonlinear equation since W in (6) and B in (24) are functions of the interfacial temperature gradients evaluated from both sides of the interface. We use the following expression to indicate the nonlinear property, To solve such a nonlinear equation, we use a Picard iteration scheme for which at each Picard iteration a linear advection-diffusion equation is solved in the spacetime framework locally. More details of the nonlinear system solver is discussed in section 5.1.1. In multiple dimensions especially, the level set equation (5) needs to be solved with spectral accuracy, which can be done by the method proposed by Sussman and Hussaini [55].
5.1.1. Nonlinear system solver. We use a Picard iteration scheme to solve the nonlinear system. The matrix of the linearized system in step 4 of Algorithm 1 is nonsymmetric. The solver we applied in our numerical tests is the direct solver from the LAPACK library. LAPACK is short for Linear Algebra PACKage, which provides routines for solving systems of linear equations, eigenvalue problems and singular value problems [4]. We predict that many iterative solvers for non-symmetric matrices can be applied to solve the linearized system, i.e., GMRES and BiCGSTAB. However, an efficient preconditioner needs to be developed for those iterative solvers. Instead of a Picard iteration method, a Newton type method can be applied to solve the nonlinear system, One obtains the following by applying the Newton method, where ∂f ∂θ is defined as follows, This method will be the same as the Picard iteration when ∂A(θ) ∂θ is zero. To avoid the evaluation of this derivative, we chose the Picard iteration rather than a Newton type method. Also, the nonlinearity in the transformed equations is not strong so that a few iterations in the Picard iteration scheme leads to convergence. Algorithm 1 the space-time discontinuous Galerkin spectral element method for solving the Stefan problem in a space-time slab E n .
1: Given θ n,s , θ n,l and φ n at time t n in the space-time slab E n . 2: Compute the normal velocity V n n at the interface from Eq. (3). Details are given in section 5.4. 3: Define a space-time slabẼ n in the moving frame of reference by introducing the backward transformationx = M b (x, t) (Eq. (8)) and forward transformation (10)). Here, the temperature values in the solid and liquid phases are denoted asθ n,s andθ n,l ,respectively. 4: InẼ n , apply a Picard iteration scheme to solve (23) in each phase, and let the initial guess of the normal velocity to be V > tol) do i Use the DG spectral element method to solve the linear equation (27) in each phase. The discretization of (27) will be discussed in section 5.3. ii Applying the forward transformation 6: end while 7: Update the level set function, φ j (j = 1, ..., p (t) ), using the forward transformation. 8: return θ n+1,s , θ n+1,l , V n+1 n and φ n+1 . 5.2. Function spaces and notations. First, we introduce a transformation,x = X(ξ, η), which defines a mapping that connects a one spatial dimension space-time elementK (x = (x, t) ∈K ) in the moving frame of reference and the reference square R = [−1, 1] × [−1, 1] ((ξ, η) ∈ R). In the space-time domainẼ ∈ R 2 , each space-time element in the moving frame of reference is a square. So the transformation (x, t) = X(ξ, η) [36] is linear in each coordinate direction and is written as where {x j } j=1,...,4 are four corners of a space-time elementK , which are numbered counter clockwise. Let's introduce two function spaces, V h and Σ h , associated with the tessellatioñ T h defined in the moving frame of reference, where P p (R) is the set of all polynomials of degree at most p = (p (x) , p (t) ) on R, with p (x) in the spatial direction and p (t) in the time direction.
Due to the discontinuous approximation spaces, the temperature and the temperature gradient approximations are double valued across element boundaries. Thus, we introduce the average {{·}} and jump · operators. Considering two adjacent elementsK n j andK n j+1 such thatQ n j,j+1 =K n j ∩K n j+1 . LetnK n j andnK n j+1 denote the corresponding outward unit normal ofK n j andK n j+1 onQ n j,j+1 . Then the average {{·}} and jump · operators are defined as, where ν +K  23), we follow the same approach presented in [15]. By introducing an auxiliary variableσ, we rewrite Eq. (23) into a first order system Multiply (37) with a test function τ ∈ Σ h , and then integrate over a space-time elementK s , we obtain After integration by parts twice, we have the following formulation, wheren is the unit normal at the interface pointing in the direction of interface movement in a space-time slabẼ n ,θ D denotes the numerical flux, which should be carefully defined to ensure stability of the method. Similarly, the weak formulation of (38) for a space-time elementK s is obtained by first multiplying (38) with a test function ν ∈ V h , and then performing integration by parts once.
whereθ A andσ D are the numerical fluxes. We separate the numerical fluxes into an advection fluxθ A and the diffusion fluxes (θ D ,σ D ). For the advection flux, we apply the upwind flux, as described in [53].
For the diffusion fluxes, we choose the LDG [15] fluxes with the stabilization parameter α j = 0,θ where β · n = sign(v · n)/2 [14], with v = (1, 1) . Note that special care is needed when imposing the boundary conditions througĥ θ D . The reason is that the boundary conditions must be the exact values imposed in the moving frame of reference, which are obtained by applying the backward transformation on the boundary conditions imposed in the fixed frame of reference. Since the backward transformation (Eq. (8)) depends on the interface velocity, we have 5.4. Calculating the interface velocity. In the Stefan problem, the phase boundary is evolving in time with the velocity defined in Eq. (3). The accuracy of the interface velocity directly affects the accuracy of the interface position. The interface velocity is a function of the temperature gradient evaluated at both sides of the interface. The order of accuracy for approximating the interface velocity must be commensurate with the order for approximating the solutions of the heat equations, otherwise the overall order of accuracy will reduce to the minimum of the velocity or temperature order. In order to obtain overall spectral accuracy, the interface velocity must be computed by a spectrally accurate method. We propose two options for computing the interface velocity. We find that each option leads to a spectrally accurate method for computing the interface velocity.
Since the velocity is determined by the temperature gradient evaluated at the interface in both phases, one option is to approximate the interface normal velocity using the following expression wheren is the unit normal at the interface pointing in the direction of interface movement in a space-time slabẼ n . Note thatσ is the auxiliary variable introduced in (37 where n is the unit normal at the interface pointing in the direction of interface movement in the fixed frame of reference. We refer to Eq. (46) as the weak form method.
Another option is to apply the normal probe method in the fixed frame of reference. In a space-time slab E n ∈ R 2 , the interface curve Q n Γ is updated by the forward transformation, while the temperature values θ(x, t) in the fixed frame of reference are obtained by the backward transformation.
Note that the transformations (11) and (12) are discretized by a (p (t) + 1)-point Gaussian quadrature rule, where p (t) is the polynomial order in the time direction. In order to obtain the temperature values θ(x, t) in the fixed frame of reference, we first find the departure points (i.e.,x) of the Legendre Gauss-Lobatto nodes (i.e., x) in the fixed frame of reference by the backward transformation (11), then the temperature values θ(x, t) are obtained via a p (x) -th order accurate interpolation, such that where I m is the interpolation operator due to the mapping, and p (x) is the polynomial order in the spatial direction. According to the zero level set at time t j (j = 1, ..., p (t) ), we construct a fictitious element in each phase with the size h that is the length of a regular space-time element in the spatial direction, i.e., Fig. 5 is an example of the normal probe method applied at time t p (t) . In each fictitious element, the temperature values, θ f ic (x, t), at its Gauss-Lobatto nodes are obtained by where I n is the interpolation operator due to the normal probe method. The normal interface velocity is then computed by, V n = 1 ρL ∇ I n (I m (θ n+1,s )) − ∇ I n (I m (θ n+1,l )) · n.
The aim of using fictitious elements is to avoid computing the temperature gradient in a small cell, which can eliminate the instability caused by computing the temperature gradient directly (i.e., spectral collocation discretization) in a small cell. The description of the normal probe method is illustrated in Fig. 5.   6. Numerical experiments. In this section, we test the space-time DG spectral element method for solving the Stefan problem on two benchmark problems in one spatial dimension. The numerical results are compared with the exact solutions in order to examine the convergence rate both with respect to space and time. The time step size is restricted by the following CFL criterion where C is a constant and set to be less than 1.0, and h is the length of a regular space-time element. This CFL condition (51) makes sure that the interface does not move across two regular space-time elements in one space-time slab. As we mentioned in section 5.1, we apply a direct solver from the LAPACK library in order to invert the resulting linear system in each space-time slab. In Figure 5. Description of the normal probe method applied in a space-time slab E n ∈ R 2 in the fixed frame of reference. The polynomial order, p = (p (x) , p (t) ), in each space-time element is chosen to be (3,2). We assume that the solid phase is on the left hand side of the interface Q n Γ while the liquid is on the other side. The fictitious element in the solid phase at time t n+1 is denoted as e n+1,s f ic , and e n+1,l f ic is the fictitious element for the liquid phase at time t n+1 . The length of both fictitious elements is h, which is the same as the length of a regular space-time element, i.e., the length of K n,s j−1 . The number of Gauss-Lobatto nodes in each fictitious element is p + 1, where p = p (x) . The temperature gradient is evaluated at the node with symbol " " in each fictitious element. the last space-time slab E E (t) −1 , the errors are measured in the L ∞ norms at time t E (t) , i.e., t E (t) = T , where T is the finial computational time.
where θ E (t) i denotes the exact temperature solution evaluated at node x i , and (θ h ) E (t) i is the approximation temperature evaluated at the same node.
where the normal velocity V n is a constant, and we take V n = 1.0. The interface is governed by Γ(t) = V n t.
In this test, the ice phase is on the left and motionless, existing at the melting temperature, i.e., θ s = θ m = 0, while the liquid phase is on the right and initialized as a supercooled liquid. The interface between the ice and liquid phase is initialized at the origin, Γ(0), and moves to the right with the normal velocity V n defined in Eq. First, we demonstrate the spectral accuracy of the approximation by plotting the maximum errors in the temperature and the interface velocity as a function of polynomial order, p = (p (x) , p (t) ), respectively. The simulation is computed over the time t = 0 to t = 0.5 with the time step ∆t = 3.84 × 10 −2 . The polynomial orders in space and in time are chosen to be the same, i.e., p (x) = p (t) . The computational domain is divided into E (x) = 5 and E (x) = 20 elements. Two options for approximating the interface velocity, i.e., the weak form method (46) and the normal probe method (50) are compared in Fig. 6. Through the comparison of the temperature (left) and the interface velocity (right), we see that both options show exponential convergence of the computed temperature to the exact solution, as well as the interface velocity.
The overall accuracy of our method is spectrally accurate. However, the convergence rates of temperature and the interface velocity are different as shown in Fig. 6. This is because the velocity is determined by the temperature gradient evaluated at the interface in both phases. The convergence rate of interface velocity is then the same as the one of temperature gradient, which is one order less than the convergence rate of temperature.
The errors stagnate after polynomial order 6 for the test cases with E (x) = 20 (Fig.6) because the error hits the asymptotic region where the truncation error dominates the computations. The error stagnates around 10 −13 ; we hypothesize that this is due to the accumulation of round-off errors which interact with the error in the position and weights for the Legendre-Gauss points. The Legendre-Gauss points are computed by the method described by Kopriva [36]. The errors in both Note that the DG spectral element method is applied on the transformed heat equations in both phases in the moving frame of reference (Algorithm 1). If we transform the exact solution (53) into the moving frame of reference by the backward transformation (11), we havẽ So the convergence rate does not depend on the polynomial order p (t) in this test problem, which is shown in Figs. 7 and 8. We choose three values of the polynomial order in time, i.e., p (t) = p (x) , p (t) = 2 and p (t) = 10. All cases lead to the same error behavior, which is consistent with the exact solution (54) in the moving frame of reference. Fig. 7 applies the weak form method (46) to calculate the interface velocity, while Fig. 8 applies the normal probe method (50). In Table 1 and Table 2, we list the maximum errors of the temperature, interface position and interface velocity. We also report the number of required Picard iterations in the last space-time slab as a function of the polynomial order in Table 1 and Table 2. The computational domain is divided into E (x) = 5 elements in the spatial direction, and the time step is chosen to be ∆t = 3.84×10 −2 . We choose the tolerance of the Picard iteration to be 10 −13 (Algorithm 1). The errors in Table 1 are obtained using the weak form velocity algorithm (46), while Table 2 uses the normal probe velocity algorithm (50). From Table 1 and Table 2, we see that only a few (i.e., less than 15) iterations lead to convergence no matter which option we applied to compute the interface velocity.  Figure 6. Errors in the temperature (left) and the interface velocity (right) as a function of polynomial order p (x) in the liquid phase. The polynomial order in time direction, p (t) , is chosen to be the same as p (x) . The comparison is made between the weak form method (46) and the normal probe method (50). The simulation is computed over the time t = 0 to t = 0.5 with the time step ∆t = 3.84 × 10 −2 .  Figure 7. Errors in the temperature (left) and the interface velocity (right) as a function of polynomial order p (x) in the liquid phase. The simulation is computed over the time t = 0 to t = 0.5. The computational domain is divided into E (x) = 5 in spatial direction, with E (x),l = 4 in the liquid in the beginning and E (x),l = 2 in the liquid at the end. The interface velocity is computed by the weak form method (46). The time step is chosen to be ∆t = 3.84 × 10 −2 .  Figure 8. Errors in the temperature (left) and the interface velocity (right) as a function of polynomial order p (x) in the liquid phase. The simulation is computed over the time t = 0 to t = 0.5. The computational domain is divided into E (x) = 5 in spatial direction, with E (x),l = 4 in the liquid in the beginning and E (x),l = 2 in the liquid at the end. The interface velocity is computed by the normal probe method (50). The time step is chosen to be ∆t = 3.84×10 −2 . Table 1. Errors of the temperature ( Err θ ∞ ), interface position ( Err Γ ∞ ) and interface velocity ( Err Vn ∞ ) in the last spacetime slab. The interface velocity is computed by applying the weak form method (46). The number of the Picard iterations is listed in the last column with tol = 10 −13 .
3.89E-003 5.29E-003 1.34E-002 12 (3,1) 1.55E-004 1.85E-004 5.08E-004 12 (4,1) 3.73E-006 4.16E-006 1.25E-005 10 (5,1) 6.80E-008 7.32E-008 2.29E-007 6 (6,1) 9.97E-010 1.05E-009 3.47E-009 5 Now, we fix the polynomial order p but vary the grid size h. The results are shown in Table 3. From Table 3, we observe that the L ∞ norm of the temperature distribution is of order p when polynomials of degree (p, p) are used. This is referred to as "sub-optimal" convergence since the convergence rate is not the classical rate for a pth order approximation which is (p + 1, p + 1). The reason for sub-optimal convergence is that the accuracy of the temperature distribution is impacted by the accuracy of the interface velocity, which in turn is determined by the temperature gradient evaluated at the interface separating the two phases. In other words, the interface velocity, which is computed by (3), is a polynomial of degree p − 1 if the temperature is approximated by a polynomial of degree p. Thus, the overall accuracy is of order p when polynomials of degree (p, p) are used. Table 3. Order of convergence of the temperature ( Err θ ∞ ) for both options in the last space-time slab.
Order of convergence p = (p (x) , p (t) ) E (x) Weak form (46) [3,62], where α = k s /(ρ s c s v ), θ wall is the Dirichlet boundary condition at x = 0, erf(·) is the error function and λ is a solution to the transcendental equation.
In this test, the ice phase is on the left, while the liquid phase is on the right and motionless, existing at the melting temperature, i.e., θ l = θ m . The interface between the ice and liquid phase is initialized at Γ(0) = 0.1 and moves to the right with the normal velocity V n defined in Eq. (3). Dirichlet boundary conditions for temperature distributions are enforced on the domain's boundary ∂Ω using the exact solution.
We set ρ l = ρ s and c l v = c s v . Now, we derive the following non-dimensional form, with Dirichlet boundary conditions Here ζ = x/l, τ = α(t/l 2 ), Θ = (θ − θ m )/∆θ, ζ Γ = Γ/l, K ls = k l /k s and St = (c v ∆θ)/L, where l is the length of the domain and set to be 1, ∆θ = θ wall − θ m and Γ is the phase boundary position. Note that St is called the Stefan number. First, we demonstrate the spectral accuracy of our method by plotting the maximum errors in the temperature and the interface velocity as a function of polynomial order, p = (p (x) , p (t) ), respectively. The computational domain is divided into E (x) = 6, and E (x) = 24 elements. The Stefan number St is set to be 0.02, and the interface moves from x = 0.1 to x = 0.3. Two options for approximating the interface velocity, i.e., the weak form method (46) and the normal probe method (50) are compared in Fig. 9. In Fig. 9, the left part reveals the exponential convergence of the temperature, while the exponential convergence of the interface velocity is observed from the right part. We see that both options exhibit spectral accuracy.
Normalprobe, E (x) =24 Figure 9. Errors in the temperature (left) and the interface velocity (right) as a function of polynomial order p (x) in the ice phase. The comparison is made between the weak form method (46) and the normal probe method (50). The interface starts at x = 0.1 and moves to x = 0.3. The polynomial order in time direction, p (t) , is chosen to be 10 so that the temporal errors are negligible. The time step is chosen to be ∆t = 4.9 × 10 −2 . The Stefan number St is set to be 0.02. Now, we examine maximum errors in the temperature and the interface velocity as a function of temporal order p (t) . Fig. 10 applies the weak form method (46) to compute the interface velocity, while Fig. 11 uses the normal probe method (50). The number of the space-time slab is denoted by E (t) , which is chosen to E (t) = 41, E (t) = 81 and E (t) = 162. The interface starts at x = 0.1 and moves to x = 0.3. Fig. 10 and Fig. 11 are shown the exponential convergence in time by fixing the polynomial order in space, p (x) , to be 10 so that the spatial errors are negligible. In Table 4 and Table 5, we show the maximum errors of the temperature, interface position and interface velocity. We also report the number of required Picard iterations in the last space-time slab as a function of the polynomial order in Table 4 and Table 5. The computational domain is divided into E (x) = 6 elements in the spatial direction, and the time step is chosen to be ∆t = 4.9 × 10 −2 . We choose the tolerance for the Picard iteration to be 10 −15 (Algorithm 1). The errors in Table 4 are obtained using the weak form method (46) for computing the interface velocity, while Table 5 uses the normal probe method (50). From Table 1 and Table  2, it shows that only a few (i.e., less than 10) iterations lead to convergence for different combinations of the polynomial order no matter which option we applied to compute the interface velocity.
Note that the solution to the Stefan problem is quasi-steady. For example, the interface moves relatively rapidly if the Stefan number St is large (see [48]). We vary the Stefan number, i.e., St = 0.001, St = 0.02, St = 0.05 and St = 0.08. Fig. 12 applies the weak form method (46) to compute the interface velocity, while Fig. 13 uses the normal probe method (50). The interface starts at x = 0.1 and moves to x = 0.4. Fig. 12 and Fig. 13 are shown the exponential convergence for different values of St.
7. Summary and conclusions. A space-time DG spectral element method for solving the Stefan problem has been presented in one spatial dimension. The Normalprobe, E (t) =81 Normalprobe, E (t) =162 Figure 11. Errors in the temperature (left) and the interface velocity (right) as a function of polynomial order p (t) in the ice phase. The interface velocity is computed by the normal probe method (50). The interface starts at x = 0.1 and moves to x = 0.3. The comparison is made among different numbers of the space-time slab, i.e., E (t) = 41, E (t) = 81 and E (t) = 162. The corresponding time step is ∆t = 4.9 × 10 −2 , ∆t = 2.5 × 10 −2 and ∆t = 1.2 × 10 −2 . The polynomial order in spatial direction, p (x) , is chosen to be 10 so that the spatial errors are negligible. The computational domain is divided into E (x) = 5 in the spatial direction. The Stefan number St is set to be 0.02. Table 4. Errors of the temperature ( Err θ ∞ ), interface position ( Err Γ ∞ ) and interface velocity ( Err Vn ∞ ) in the last spacetime slab. The interface velocity is computed by applying the weak form method (46). The computational domain is divided into E (x) = 6 in the spatial direction. The interface starts at x = 0.1 and moves to x = 0.3. The time step is chosen to be ∆t = 4.9 × 10 −2 . The number of the Picard iteration is listed in the last column with tol = 10 −15 .  Table 5. Errors of the temperature ( Err θ ∞ ), interface position ( Err Γ ∞ ) and interface velocity ( Err Vn ∞ ) in the last spacetime slab. The interface velocity is computed by applying the normal probe method (50). The computational domain is divided into E (x) = 6 in the spatial direction. The interface starts at x = 0.1 and moves to x = 0.3. The time step is chosen to be ∆t = 4.9 × 10 −2 . The number of the Picard iteration is listed in the last column with tol = 10 −15 .   curved subelements in a cut-cell are handled by mapping the fixed frame of reference to a moving frame of reference. In the moving frame of reference, we have a rectangular mesh in each space-time slab. In each space-time slab, a DG spectral element method is applied to solve the transformed heat equations in the space-time formulation in each phase independently. We present two options for computing the interface velocity: the weak form option and the normal probe option. Numerical experiments on standard one-dimensional Stefan problems demonstrate the spectral convergence in both space and time for the temperature and the interface velocity. Both options for calculating the interface velocity exhibit spectral convergence in both space and time. Moreover, the Picard iteration convergences only in a just iterations (less than 15) for both options. We have also analyzed the present method for different values of the Stefan number, St, and found that the space-time spectral accuracy property of the present method is not sensitive to the values of St. The present method is not restricted to one spatial dimension, and its extension to multiple space dimensions is the subject of future activity. In multiple space dimensions, a high-order accurate numerical quadrature scheme [50] needs to be developed for computing the integrals over curved and implicitly defined surfaces and volumes in the spatial element, generated by the phase interface, in the weak formulations. Future work will also explore an alternative direct solver methodology [42] for the sparse linear systems resulting from each Picard iteration.