\documentclass[reqno]{amsart}
\usepackage{hyperref}

\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2016 (2016), No. 258, pp. 1--8.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu}
\thanks{\copyright 2016 Texas State University.}
\vspace{8mm}}

\begin{document}
\title[\hfilneg EJDE-2016/258\hfil Sharp interface limit]
{Sharp interface limit of a homogenized phase field model for
phase transitions in porous media}

\author[M. H\"opker \hfil EJDE-2016/258\hfilneg]
{Martin H\"opker}

\address{Martin H\"opker \newline
Center for Industrial Mathematics,
University of Bremen,
Postfach 33 04 40, 28334 Bremen, Germany}
\email{hoepker@math.uni-bremen.de}

\thanks{Submitted August 19, 2016. Published September 22, 2016.}
\subjclass[2010]{80A22, 35Q79, 35B27}
\keywords{Stefan Problems; asymptotic expansions; phase field models;
\hfill\break\indent  partial differential equations; homogenization}

\begin{abstract}
 A homogenized phase field model for phase transitions in porous media
 is considered. By making use of the method of formal asymptotic expansion
 with respect to the interface thickness, a sharp interface limit problem
 is derived. This limit problem turns out to be similar to the classical
 Stefan problem with surface tension and kinetic undercooling.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{remark}[theorem]{Remark}
\allowdisplaybreaks

\section{Introduction}

The study of the phase change between water and ice in porous media plays 
a vital role in the understanding of important phenomena like the frost 
attack on concrete or the thawing of permafrost soil. 
In \cite{Hoepker2016}, a mathematical microscale model of this process,
based on the standard Caginalp phase field model for phase transitions, 
is presented and homogenized via two-scale convergence. 
Let $\Omega\subset\mathbb{R}^3$ be a bounded Lipschitz domain that represents 
a porous body, and let $S:=(0,T)$ be a time interval. The homogenized problem 
is of the form: Find $u$, $\theta$, and $\chi$ such that
\begin{subequations}\label{subeq:homogenized_classical}
 \begin{gather}
 |{Z^{\mathrm{s}}}|\rho_{\mathrm{s}} u'-\operatorname{div}(P^{\mathrm{s}} \nabla u) 
+ |\Gamma|\kappa_{RI}(u-\theta)=0, \label{eq:homogenized_classical_heat_solid_matrix}\\
  |{Z^{\mathrm{p}}}| \rho_{\mathrm{p}} \theta' 
+ |{Z^{\mathrm{p}}}| \frac{\lambda}{2} \chi'-\operatorname{div}
(P^{\mathrm{p}} \nabla \theta) +|\Gamma|\kappa_{RI}(\theta-u)=0, \label{eq:homogenized_classical_heat_porespace}\\
  |{Z^{\mathrm{p}}}| \alpha\xi^2 \chi' - \xi^2 
\operatorname{div}( P \nabla \chi) + |{Z^{\mathrm{p}}}| 
\frac{1}{2a} ( \chi^3 - \chi )
= |{Z^{\mathrm{p}}}| 2 \theta \label{eq:homogenized_classical_heat_phasefield}
 \end{gather}
\end{subequations}
in $S\times\Omega$, supplemented by exchange boundary conditions and initial 
conditions.

In the system \eqref{subeq:homogenized_classical}, $u$ and $\theta$ are the 
effective temperatures in the solid matrix and in the pore space of the 
porous medium, respectively, and $\chi$ is the effective phase field variable. 
The constant parameter $\xi$ represents the small interface thickness of the 
phase field.
The coefficients $\rho_{\mathrm{s}},\rho_{\mathrm{p}}, \kappa_{RI}, 
\lambda, \alpha, a$ are positive and constant parameters, and 
$P^{\mathrm{s}}$, $P^{\mathrm{p}}$, $P$ are homogenized, positive definite tensors. 
The constants $|{Z^{\mathrm{s}}}|$ and $|{Z^{\mathrm{p}}}|$ represent the volume 
of the solid matrix and of the pore space in the reference cell in the 
homogenization, and $|\Gamma|$ denotes the surface measure of their interface.
For the details of the modeling and the homogenization that lead to 
\eqref{subeq:homogenized_classical}, we refer to \cite{Hoepker2016}.
We also refer to \cite{Caginalp1986}, \cite{Caginalp1989}, and
\cite{Visintin1996} for an introduction to phase field models for phase
transitions, as well as to \cite{Showalter1982} for the similar concept
of \textit{mushy regions}.

In this article, we apply the method of formal asymptotic expansion, with 
respect to the interface thickness $\xi$, to derive a sharp interface model 
from the homogenized phase field model \eqref{subeq:homogenized_classical}. 
To this end, we follow the lines of \cite{Caginalp1989}, where the main
difference lies in the homogenized tensor $P$ in equation 
\eqref{eq:homogenized_classical_heat_phasefield}, which replaces the identity
 matrix from the standard phase field model. The sharp interface model turns 
out to be similar to the classical \textit{Stefan problem} 
(see, e.g., \cite{Stefan1891}, \cite{Meirmanov1992}, \cite{Visintin1996},
and \cite{Eck2008}) with surface tension and kinetic undercooling but involves
some additional terms.

This article is organized as follows:
 In \ref{section:scaling}, the system \eqref{subeq:homogenized_classical}
is appropriately scaled. In \ref{section:interface_coordinates}, we
introduce a local coordinate system near the phase interface.
Sections \ref{section:outer_expansions}, \ref{section:inner_expansions} are concerned 
with the asymptotic expansions of the variables $u$, $\theta$, and $\chi$ 
outside of and near the interface. Limit equations are derived from the 
system \eqref{subeq:homogenized_classical} by neglecting terms of high order. 
Finally, in \ref{section:full_sharp_interface_model}, the full sharp
interface limit model is stated.


\section{Scaling}\label{section:scaling}

Similarly as in \cite[Section IV]{Caginalp1989}, we will consider the limit 
$\xi\to 0$ and $a\to 0$ with fixed $\frac{\xi}{\sqrt{a}}$. 
To this end, we introduce the scaling $\epsilon=\xi^2$ and 
$\frac{\epsilon}{a}=C_0$ in equation \eqref{eq:homogenized_classical_heat_phasefield},
 which yields
\begin{equation}\label{eq:phase_field_only_scaled}
|{Z^{\mathrm{p}}}| \alpha \epsilon^2 \chi' 
- \epsilon^2 \operatorname{div}( P \nabla \chi) + |{Z^{\mathrm{p}}}| 
\frac{C_0}{2}( \chi^3 - \chi ) = \epsilon |{Z^{\mathrm{p}}}| 2 \theta.
\end{equation}
We define $g(\chi):= |{Z^{\mathrm{p}}}| \frac{C_0}{2} (\chi-\chi^3)$ and, 
for simplicity, also rename the constants in the equations 
\eqref{eq:homogenized_classical_heat_solid_matrix}, 
\eqref{eq:homogenized_classical_heat_porespace}, and 
\eqref{eq:phase_field_only_scaled} such that we obtain the scaled system
\begin{subequations}
\begin{gather}
 \rho_{\mathrm{s}} u'-\operatorname{div}(P^{\mathrm{s}} \nabla u) 
+ \kappa_{RI}(u-\theta) =0 \label{eq:sharp_matrix_heat_equation_renamed} \\
\rho_{\mathrm{p}} \theta' + \frac{\lambda}{2} \chi' 
- \operatorname{div}( P^{\mathrm{p}} \nabla \theta ) + \kappa_{RI}(\theta-u) = 0. \label{eq:sharp_heat_equation_renamed} \\
 \alpha \epsilon^2 \chi' - \epsilon^2 \operatorname{div}(P \nabla \chi) 
= g(\chi) + \epsilon \beta \theta \quad\text{in }\Omega\times S. 
\label{eq:sharp_phasefield_equation_renamed}
\end{gather}
\end{subequations}
Note that $\theta,\chi$, and $u$ depend on the scaling parameter $\epsilon$.

\begin{remark} \label{rmk2.1} \rm
Different scalings of phase field models may yield different asymptotic limits. 
See, e.g., \cite{Caginalp1989} and \cite{Caginalp1998}, where various
limit models are derived from the Caginalp phase field model.
\end{remark}

\section{Interface and Local Coordinates}\label{section:interface_coordinates}

Using the phase field variable $\chi$, we define the interface between the 
solid and the liquid phase at time $t$ to be given by the set
$$
\Gamma(t,\epsilon):=\{x\in\Omega: \chi(t,x)=0\}.
$$
In the following, we assume that $\Gamma(t,\epsilon)$ is a closed, connected, 
and sufficiently smooth surface. A point $x\in\Omega$ near $\Gamma(t,\epsilon)$ 
can then be represented in the local orthogonal coordinate system $(r,s_1,s_2)$, 
which is often used when deriving sharp limits of phase field models 
(see, e.g., \cite{Caginalp1989,Caginalp2005, Fife1995}):
By $r(t,x,\epsilon)$, we denote the signed distance of $x$ to the interface 
$\Gamma(t,\epsilon)$ (which is chosen positive if $x$ lies in the liquid phase), 
and $s_1(t,x,\epsilon)$ and $s_2(t,x,\epsilon)$ are measures of arc length along 
$\Gamma(t,\epsilon)$ from a reference point. We write $s=(s_1~s_2)$. 
The coordinate system $(r,s)$ satisfies
$|\nabla r|=1$ and $\Delta r =\kappa$ near $\Gamma(t,\epsilon)$, 
where $\kappa$ is the mean curvature of $\Gamma(t,\epsilon)$
(The sign of the mean curvature is chosen such that $\kappa>0$ 
for a strictly convex solid phase.).
Furthermore, $\nabla r \cdot \nabla s_i = 0$ for $i=1,2$.
We assume the following asymptotic expansions of the coordinates:
\begin{gather*}
 r(t,x,\epsilon)=r^0(t,x)+\epsilon r^1(t,x)+\epsilon^2r^2(t,x)+\dots, \\
 s_i(t,x,\epsilon)=s_i^0(t,x)+\epsilon s_i^1(t,x)+\epsilon^2 s_i^2(t,x)+\dots, 
\quad\text{for }i=1,2.
\end{gather*}
We assume that, for $\varepsilon\to 0$, the interfaces $\Gamma(t,\epsilon)$ 
converge (uniformly in distance) to a closed, connected, and smooth 
interface $\Gamma^0(t)$ that has a neighbourhood parametrized by $r^0$ 
and $s^0=(s_1^0~s_2^0)$. The liquid and the solid phase, which are seperated
 by $\Gamma^0(t)$, are denoted by $\Omega_l(t)$ and $\Omega_s(t)$, respectively. 
By $\vec{n}(t,x)$, we denote the unit normal vector on $\Gamma^0(t)$ 
at point $x$ into the liquid phase $\Omega_l(t)$.

Let $A\in\mathbb{R}^{3\times 3}$ be a constant and symmetric matrix, and 
let $u\colon S\times\mathbb{R}^3\times\mathbb{R}\to\mathbb{R}$. 
We consider the coordinate change
$$
u(t,x,\epsilon)=\tilde{u}(t,r,s,\epsilon).
$$
For later reference, we note that this transformation yields
\begin{equation}\label{eq:divergence_transformation}
\begin{aligned}
 \operatorname{div}(A\nabla u)
&= \frac{\partial^2 \tilde{u}}{\partial r^2} \nabla r^\mathrm{T}A\nabla r
 +\frac{\partial \tilde{u}}{\partial r}\operatorname{div}(A\nabla r)
 +2 \sum_{l=1}^{2} \frac{\partial^2 \tilde{u}}{\partial r \partial s_l}
\nabla s_l^\mathrm{T}A\nabla r \\
&\quad + \sum_{l=1}^{2} \frac{\partial \tilde{u}}{\partial s_l} 
\operatorname{div}( A \nabla s_l)
 +\sum_{l,m=1}^{2} \frac{\partial^2 \tilde{u}}{\partial s_l \partial s_m} 
\nabla s_m^\mathrm{T}A\nabla s_l
\end{aligned}
\end{equation}
and
\begin{equation}\label{eq:time_derivative_transformation}
 \frac{\partial u}{\partial t}
=\frac{\partial \tilde{u}}{\partial t}+\frac{\partial r}{\partial t}
\frac{\partial \tilde{u}}{\partial r}+\sum_{l=1}^{2}
\frac{\partial s_l}{\partial t}\frac{\partial \tilde{u}}{\partial s_l}.
\end{equation}


\section{Outer Expansions}\label{section:outer_expansions}

Outside of $\Gamma(t,\epsilon)$, we assume the following \emph{outer expansions} 
in the cartesian coordinate system
\begin{subequations}
\begin{gather}
 u(t,x,\epsilon)=u^0(t,x)+\epsilon u^1(t,x)+\epsilon^2 u^2(t,x)+\dots, 
\label{eq:outer_expansion_u}\\
 \theta(t,x,\epsilon)=\theta^0(t,x)+\epsilon\theta^1(t,x)
+\epsilon^2\theta^2(t,x)+\dots, \label{eq:outer_expansion_theta}\\
 \chi(t,x,\epsilon)=\chi^0(t,x)+\epsilon\chi^1(t,x)+\epsilon^2\chi^2(t,x)+\dots.
 \label{eq:outer_expansion_chi}.
\end{gather}
\end{subequations}
The terms on the right-hand side of the equations \eqref{eq:outer_expansion_theta} 
and \eqref{eq:outer_expansion_chi} are assumed to be twice differentiable in 
$\Omega\backslash\Gamma^0(t)$ but may be discontinuous across the limit phase 
interface $\Gamma^0(t)$. Since the variable $u$ is not directly connected to 
the phase change, we assume the $u_k$, for all $k\in\mathbb{N}$, to be twice 
differentiable in the whole domain $\Omega$.
From the order $\mathcal{O}(1)$ of \eqref{eq:sharp_matrix_heat_equation_renamed}, 
we deduce the equation
\begin{equation}\label{eq:sharp_matrix_heat_equation_renamed_o1}
 \rho_{\mathrm{s}} {u^0}'-\operatorname{div}(P^{\mathrm{s}} \nabla u^0) 
+ \kappa_{RI}(u^0-\theta^0)=0,
\end{equation}
which is valid in the whole domain $\Omega$. By using the expansion of $\chi$, 
it is easy to obtain that
\begin{equation}\label{eq:expansion_of_g}
 \begin{aligned}
 g(\chi)&=|{Z^{\mathrm{p}}}| \frac{C_0}{2}(\chi^0 - (\chi^0)^3)
 +\epsilon |{Z^{\mathrm{p}}}| \frac{C_0}{2} ( \chi^1 - 3 (\chi^0)^2\chi^1) 
+ \mathcal{O}(\epsilon^2) \\
 &=g(\chi^0) +\epsilon g'(\chi^0)\chi^1 + \mathcal{O}(\epsilon^2).
 \end{aligned}
\end{equation}
Equations \eqref{eq:sharp_heat_equation_renamed} and 
\eqref{eq:sharp_phasefield_equation_renamed} thus yield the order 
$\mathcal{O}(1)$ equations
\begin{equation}\label{eq:sharp_heat_equation_renamed_o1}
\rho_{\mathrm{p}} {\theta^0}' + \frac{\lambda}{2} {\chi^0}' 
- \operatorname{div}( P^{\mathrm{p}} \nabla \theta^0 ) 
+ \kappa_{RI}(\theta^0-u^0) = 0
\end{equation}
and
\begin{equation}\label{eq:gchizero}
 g(\chi^0)=0,
\end{equation}
respectively. Note that equation \eqref{eq:sharp_heat_equation_renamed_o1} 
is only valid on each side of $\Gamma^0(t)$, due to the possible discontinuity 
of $\theta^0$ and $\chi^0$ across $\Gamma^0(t)$.
Equation \eqref{eq:gchizero} implies that $\chi^0$ takes one of the constant 
values $0,1$, or $-1$ on each side of $\Gamma^0(t)$.
Hence, ${\chi^0}'=0$, and we can thus deduce
\begin{equation}\label{eq:sharp_heat_equation_renamed_o1_final}
\rho_{\mathrm{p}} {\theta^0}' - \operatorname{div}( P^{\mathrm{p}} \nabla \theta^0 ) 
+ \kappa_{RI}(\theta^0-u^0) = 0
\end{equation}
from \eqref{eq:sharp_heat_equation_renamed_o1}, outside of $\Gamma^0(t)$.

Furthermore, we consider the coordinate change
\begin{gather*}
 u(t,x,\epsilon)=\tilde{u}(t,r,s,\epsilon), \\
 \theta(t,x,\epsilon)=\tilde{\theta}(t,r,s,\epsilon), \\
 \chi(t,x,\epsilon)=\tilde{\chi}(t,r,s,\epsilon)
\end{gather*}
and assume the following expansions to hold:
\begin{gather*}
 \tilde{u}(t,r,s,\epsilon)=\tilde{u}^0(t,r,s)
+\epsilon \tilde{u}^1(t,r,s)+\epsilon^2 \tilde{u}^2(t,r,s)+\dots, \\
 \tilde{\theta}(t,r,s,\epsilon)=\tilde{\theta}^0(t,r,s)
+\epsilon\tilde{\theta}^1(t,r,s)+\epsilon^2\tilde{\theta}^2(t,r,s)+\dots, \\
 \tilde{\chi}(t,r,s,\epsilon)=\tilde{\chi}^0(t,r,s)
+\epsilon\tilde{\chi}^1(t,r,s)+\epsilon^2\tilde{\chi}^2(t,r,s)+\dots.
\end{gather*}


\section{Inner Expansions}\label{section:inner_expansions}

In a neighbourhood of $\Gamma(t,\epsilon)$, we introduce a coordinate change, 
using the stretched variable $z:=\frac{r}{\epsilon}$:
\begin{gather*}
 \tilde{u}(t,r,s,\epsilon) =U(t,z,s,\epsilon), \\
 \tilde{\theta}(t,r,s,\epsilon) =\vartheta(t,z,s,\epsilon), \\
 \tilde{\chi}(t,r,s,\epsilon) =X(t,z,s,\epsilon),
\end{gather*}
and assume the \emph{inner expansions}
\begin{gather*}
 U(t,z,s,\epsilon)=U^0(t,z,s)+\epsilon U^1(t,z,s)+\epsilon^2 U^2(t,z,s)+\dots, \\
 \vartheta(t,z,s,\epsilon)=\vartheta^0(t,z,s)+\epsilon\vartheta^1(t,z,s)
+\epsilon^2\vartheta^2(t,z,s)+\dots, \\
 X(t,z,s,\epsilon)=X^0(t,z,s)+\epsilon X^1(t,z,s)+\epsilon^2 X^2(t,z,s)+\dots.
\end{gather*}
In the context of the above coordinate change, we note that 
$\frac{\partial}{\partial r}=\frac{1}{\epsilon}\frac{\partial}{\partial z}$.
By using the identities \eqref{eq:divergence_transformation} and 
\eqref{eq:time_derivative_transformation}, we can thus rewrite equation 
\eqref{eq:sharp_heat_equation_renamed} as
\begin{align*}
&\rho_{\mathrm{p}} \frac{\partial \vartheta}{\partial t}
 +\frac{1}{\epsilon}\rho_{\mathrm{p}}\frac{\partial r}{\partial t}
 \frac{\partial \vartheta}{\partial z}+\rho_{\mathrm{p}}\sum_{l=1}^{2}
 \frac{\partial s_l}{\partial t}\frac{\partial \vartheta}{\partial s_l}
+ \frac{\lambda}{2} \frac{\partial X}{\partial t}
 +\frac{1}{\epsilon}\frac{\lambda}{2}\frac{\partial r}{\partial t}
 \frac{\partial X}{\partial z}+\frac{\lambda}{2}\sum_{l=1}^{2}
 \frac{\partial s_l}{\partial t}\frac{\partial X}{\partial s_l} \\
&-\frac{1}{\epsilon^2}\frac{\partial^2 \vartheta}{\partial z^2} 
 \nabla r^\mathrm{T}P^{\mathrm{p}}\nabla r
 -\frac{1}{\epsilon}\frac{\partial \vartheta}{\partial z}
 \operatorname{div}(P^{\mathrm{p}}\nabla r)
 -\frac{1}{\epsilon}2 \sum_{l=1}^{2} 
 \frac{\partial^2 \vartheta}{\partial z \partial s_l}
 \nabla s_l^\mathrm{T}P^{\mathrm{p}}\nabla r \\
& - \sum_{l=1}^{2} \frac{\partial \vartheta}{\partial s_l} \operatorname{div}
 ( P^{\mathrm{p}} \nabla s_l)
 -\sum_{l,m=1}^{2} \frac{\partial^2 \vartheta}{\partial s_l \partial s_m} 
 \nabla s_m^\mathrm{T}P^{\mathrm{p}}\nabla s_l
+ \kappa_{RI}(\vartheta-U) = 0.
\end{align*}
Multiplying by $\epsilon^2$ yields
\begin{equation}\label{eq:porespace_heat_order_equation}
\begin{aligned}
&\epsilon^2 \rho_{\mathrm{p}} \frac{\partial \vartheta}{\partial t}
 +\epsilon\rho_{\mathrm{p}}\frac{\partial r}{\partial t}
 \frac{\partial \vartheta}{\partial z}
+\epsilon^2 \rho_{\mathrm{p}}\sum_{l=1}^{2}
 \frac{\partial s_l}{\partial t}\frac{\partial \vartheta}{\partial s_l}
+ \epsilon^2 \frac{\lambda}{2} \frac{\partial X}{\partial t}
 +\epsilon \frac{\lambda}{2}\frac{\partial r}{\partial t}
 \frac{\partial X}{\partial z}\\
&+\epsilon^2\frac{\lambda}{2}\sum_{l=1}^{2}
 \frac{\partial s_l}{\partial t}\frac{\partial X}{\partial s_l} 
 -\frac{\partial^2 \vartheta}{\partial z^2} \nabla r^\mathrm{T}
 P^{\mathrm{p}}\nabla r
 -\epsilon\frac{\partial \vartheta}{\partial z}\operatorname{div}
 (P^{\mathrm{p}}\nabla r) \\
& - \epsilon 2\sum_{l=1}^{2} \frac{\partial^2 \vartheta}{\partial z 
 \partial s_l}\nabla s_l^\mathrm{T}P^{\mathrm{p}}\nabla r 
 - \epsilon^2\sum_{l=1}^{2} \frac{\partial \vartheta}{\partial s_l} 
 \operatorname{div}( P^{\mathrm{p}} \nabla s_l)\\
& -\epsilon^2\sum_{l,m=1}^{2} \frac{\partial^2 \vartheta}{\partial s_l 
 \partial s_m} \nabla s_m^\mathrm{T}P^{\mathrm{p}}\nabla s_l
+ \epsilon^2\kappa_{RI}(\vartheta-U) = 0.
\end{aligned}
\end{equation}
By using the inner expansions, we deduce from the phase field equation 
\eqref{eq:sharp_phasefield_equation_renamed} that
\begin{equation}\label{eq:phasefield_order_equation}
\begin{aligned}
&\epsilon^2 \alpha \frac{\partial X}{\partial t} +\epsilon \alpha \frac{\partial r}{\partial t}\frac{\partial X}{\partial z}+\epsilon^2 \alpha \sum_{l=1}^{2}\frac{\partial s_l}{\partial t}\frac{\partial X}{\partial s_l}
\\
&-\frac{\partial^2 X}{\partial z^2} \nabla r^\mathrm{T}P\nabla r
 -\epsilon\frac{\partial X}{\partial z}\operatorname{div}(P\nabla r)
 -\epsilon 2 \sum_{l=1}^{2} \frac{\partial^2 X}{\partial z 
 \partial s_l}\nabla s_l^\mathrm{T}P\nabla r \\
&-\epsilon^2 \sum_{l=1}^{2} \frac{\partial X}{\partial s_l} 
 \operatorname{div}( P \nabla s_l)
 -\epsilon^2\sum_{l,m=1}^{2} \frac{\partial^2 X}{\partial s_l \partial s_m}
 \nabla s_m^\mathrm{T}P\nabla s_l \\
&= g(X)+ \epsilon \beta \vartheta.
\end{aligned}
\end{equation}
Similarly as with equation \eqref{eq:expansion_of_g}, we obtain
\begin{equation}
 g(X)=g(X^0) +\epsilon g'(X^0)X^1 + \mathcal{O}(\epsilon^2).
\end{equation}
The order $\mathcal{O}(1)$ of equation \eqref{eq:porespace_heat_order_equation} 
yields
$$
- \frac{\partial^2 \vartheta^0}{\partial z^2} \nabla {r^0}^\mathrm{T}
P^{\mathrm{p}}\nabla r^0=0.
$$
From $|\nabla r|=1$, it follows that $\nabla r^0\neq 0$. 
Since $P^{\mathrm{p}}$ is positive definite, we can hence deduce
$$
- \frac{\partial^2 \vartheta^0}{\partial z^2}=0.
$$
This yields $\vartheta^0(t,z,s)=a(t,s)z+b(t,s)$ with $a(t,s),b(t,s)\in\mathbb{R}$.
We apply the matching condition (see, e.g., \cite{Caginalp1989})
\begin{equation}\label{eq:first_matching_condition_theta}
 \tilde{\theta}^0(t,\Gamma^0_{\pm}) = \lim_{z\to\pm\infty} \vartheta^0(t,z,s),
\end{equation}
which connects the inner and the outer expansions. Here, we denote by
 $\tilde{\theta}^0(t,\Gamma^0_{\pm})$ the value of $\tilde{\theta}^0$ when 
approaching $\Gamma^0(t)$ from the liquid and the solid phase, respectively. 
From the natural assumption that $\tilde{\theta}^0$ is bounded at the interface 
$\Gamma^0(t)$, we deduce $a=0$ from the equation 
\eqref{eq:first_matching_condition_theta}. Thus,
\begin{equation}
 \vartheta^0(t,z,s)=b(t,s),
\end{equation}
independently of $z$.
Analogously to \eqref{eq:first_matching_condition_theta}, we apply the
 matching condition
\begin{equation}\label{eq:first_matching_condition_chi}
 \lim_{z\to\pm\infty} X^0(t,z,s) = \tilde{\chi}^0(t,\Gamma^0_{\pm} )=\pm 1.
\end{equation}
We note that $X(t,0,s,\epsilon)=0$ by the definition of $\Gamma(t,\epsilon)$. 
It, thus, directly follows that
\begin{equation*}
 X^0(t,0,s)=0.
\end{equation*}
We now consider the $\mathcal{O}(\epsilon)$ balance of equation 
\eqref{eq:porespace_heat_order_equation}:
\begin{align*}
&\rho_{\mathrm{p}}\frac{\partial r^0}{\partial t}
 \frac{\partial \vartheta^0}{\partial z}
+\frac{\lambda}{2}\frac{\partial r^0}{\partial t}\frac{\partial X^0}{\partial z}
- 2\frac{\partial^2 \vartheta^0}{\partial z^2} \nabla {r^1}^\mathrm{T}
 P^{\mathrm{p}}\nabla r^0
- \frac{\partial^2 \vartheta^1}{\partial z^2} \nabla {r^0}^\mathrm{T}
 P^{\mathrm{p}}\nabla r^0 \\
&-\frac{\partial \vartheta^0}{\partial z}\operatorname{div}(P^{\mathrm{p}}\nabla r^0)
-2 \sum_{l=1}^{2} \frac{\partial^2 \vartheta^0}{\partial z \partial s_l}
\nabla {s_l^0}^\mathrm{T}P^{\mathrm{p}}\nabla r^0
=0.
\end{align*}
Since $\vartheta^0$ is constant with respect to the variable $z$, we obtain
\[
\frac{\lambda}{2}\frac{\partial r^0}{\partial t}\frac{\partial X^0}{\partial z}
= \frac{\partial^2 \vartheta^1}{\partial z^2} \nabla {r^0}^\mathrm{T}
P^{\mathrm{p}}\nabla r^0.
\]
Thus, integration with respect to $z$ yields
\[
\frac{\lambda}{2}\frac{\partial r^0}{\partial t} X^0 + c
= \frac{\partial \vartheta^1}{\partial z} \nabla {r^0}^\mathrm{T}
P^{\mathrm{p}}\nabla r^0,
\]
where $c$ may depend on $t$ and $s$ but not on $z$. In the following, we will 
apply the matching condition (see, e.g., \cite{Caginalp1989})
$$
\frac{\partial \tilde{\theta}^0}{\partial r}(t,\Gamma^0_{\pm} )
=\lim_{z\to\pm\infty}\frac{\partial \vartheta^1}{\partial z}(t,z,s).
$$
We calculate
\begin{align*}
\frac{\partial \tilde{\theta}^0}{\partial r}(t,\Gamma^0_{\pm} ) 
\nabla {r^0}^\mathrm{T}P^{\mathrm{p}}\nabla r^0
&=\lim_{z\to\pm\infty}\frac{\partial \vartheta^1}{\partial z}
 \nabla {r^0}^\mathrm{T}P^{\mathrm{p}}\nabla r^0
=\frac{\lambda}{2}\frac{\partial r^0}{\partial t} \lim_{z\to\pm\infty} X^0 + c \\
&= \pm \frac{\lambda}{2}\frac{\partial r^0}{\partial t} + c.
\end{align*}
Thus,
\begin{equation}\label{eq:jump_condition_gamma0}
\nabla {r^0}^\mathrm{T}P^{\mathrm{p}}\nabla r^0 
\Big[\frac{\partial \tilde{\theta}^0}{\partial r}\Big]_{\Gamma^0_\pm}
=-\lambda V \quad\text{on } \Gamma^0(t),
\end{equation}
where $V=-\frac{\partial r^0}{\partial t}$ is the normal velocity of 
$\Gamma^0(t)$ (in the direction $\vec{n}$) and 
$[\frac{\partial \tilde{\theta}^0}{\partial r}]_{\Gamma^0_\pm}
=[\frac{\partial \theta^0}{\partial \vec{n}}]_{\Gamma^0_\pm}$ denotes 
the jump of the normal derivative $\frac{\partial \theta^0}{\partial \vec{n}}$ 
across $\Gamma^0(t)$.

We now consider the problem of order $\mathcal{O}(\epsilon)$ of equation 
\eqref{eq:phasefield_order_equation}:
\begin{equation}\label{eq:phasefield_order_equation_epsilon}
\begin{aligned}
&\alpha \frac{\partial r^0}{\partial t}\frac{\partial X^0}{\partial z}
-2\frac{\partial^2 X^0}{\partial z^2} \nabla {r^1}^\mathrm{T}P\nabla r^0
-\frac{\partial^2 X^1}{\partial z^2} \nabla {r^0}^\mathrm{T}P\nabla r^0 \\
&-\frac{\partial X^0}{\partial z}\operatorname{div}(P\nabla r^0)
- 2 \sum_{l=1}^{2} \frac{\partial^2 X^0}{\partial z \partial s_l}
 \nabla {{s_l}^0}^\mathrm{T}P\nabla r^0 \\
&= g'(X^0)X^1 + \beta \vartheta^0,
\end{aligned}
\end{equation}
to which we (see, e.g., \cite{Eck2004}) add the boundary condition
\begin{equation}\label{eq:boundary_condition_X1}
\lim_{z\to\pm\infty}X^1(t,z,s)=0,
\end{equation}
meaning that the phase field variable equals zero up to order $\varepsilon$.
The order $\mathcal{O}(1)$ problem of equation \eqref{eq:phasefield_order_equation}
 reads
\[
\frac{\partial^2 X^0}{\partial z^2} \nabla {r^0}^\mathrm{T}P\nabla r^0+ g(X^0)=0,
\]
from which we deduce, by differentiating,
\[
\frac{\partial}{\partial z}\frac{\partial^2 X^0}{\partial z^2} 
\nabla {r^0}^\mathrm{T}P\nabla r^0+ g'(X^0)\frac{\partial X^0}{\partial z}=0.
\]
Thus, using integration by parts and the condition \eqref{eq:boundary_condition_X1} 
(compare to the approach in \cite[Section 7.9]{Eck2008}), we obtain
\begin{align*}
&\int_\mathbb{R} \Big(\frac{\partial^2 X^1}{\partial z^2} \nabla 
{r^0}^\mathrm{T}P\nabla r^0+ g'(X^0)X^1\Big)\frac{\partial X^0}{\partial z}\,\mathrm{d}z \\
&=\int_\mathbb{R} \frac{\partial^2 X^1}{\partial z^2} \nabla 
{r^0}^\mathrm{T}P\nabla r^0 \frac{\partial X^0}{\partial z} 
+ g'(X^0)X^1 \frac{\partial X^0}{\partial z}\,\mathrm{d}z \\
&=\int_\mathbb{R} X^1 \nabla {r^0}^\mathrm{T}P\nabla r^0 
 \frac{\partial}{\partial z}\frac{\partial^2 X^0}{\partial z^2} 
 + g'(X^0)X^1 \frac{\partial X^0}{\partial z}\,\mathrm{d}z \\
&=\int_\mathbb{R} X^1 \Big( \nabla {r^0}^\mathrm{T}P\nabla r^0 
\frac{\partial}{\partial z}\frac{\partial^2 X^0}{\partial z^2} 
+ g'(X^0)\frac{\partial X^0}{\partial z}\Big)\,\mathrm{d}z=0.
\end{align*}
By using this relation, we directly deduce from the multiplication of 
equation \eqref{eq:phasefield_order_equation_epsilon} by 
$\frac{\partial X^0}{\partial z}$ and integration that
\begin{equation}\label{eq:final_equation_interface_temperature}
\begin{aligned}
&\Big(\alpha \frac{\partial r^0}{\partial t} - \operatorname{div}(P\nabla r^0) \Big) 
\int_\mathbb{R} |\frac{\partial X^0}{\partial z}|^2 \,\mathrm{d}z \\
&- \int_\mathbb{R}2\frac{\partial^2 X^0}{\partial z^2}
  \nabla {r^1}^\mathrm{T}P\nabla r^0 \frac{\partial X^0}{\partial z}\,\mathrm{d}z
- \int_\mathbb{R} 2 \sum_{l=1}^{2} \frac{\partial^2 X^0}{\partial z \partial s_l}
 \nabla {{s_l}^0}^\mathrm{T}P\nabla r^0 \frac{\partial X^0}{\partial z}\,\mathrm{d}z \\
&= \beta \vartheta^0 \int_\mathbb{R} \frac{\partial X^0}{\partial z}\,\mathrm{d}z.
\end{aligned}
\end{equation}

\section{The Sharp Interface Model}\label{section:full_sharp_interface_model}

We note that (see \cite{Caginalp1989}) the surface tension is given by
\begin{equation*}
 \sigma=\int_\mathbb{R} |\frac{\partial X^0}{\partial z}|^2\,\mathrm{d}z.
\end{equation*}
Furthermore, the mean curvature of $\Gamma^0(t)$ is given by $\kappa^0=\Delta r^0$. 
Hence, $\overline{\kappa}:=\operatorname{div}(P\nabla r^0)$, which equals 
$\kappa^0$ in the case $P=I$, is related to the mean curvature. 
From the relation \eqref{eq:first_matching_condition_chi}, we deduce
 $\int_\mathbb{R} \frac{\partial X^0}{\partial z}\,\mathrm{d}z=2$, 
and we can thus rewrite equation \eqref{eq:final_equation_interface_temperature} as
\begin{equation}
\begin{aligned}
&(- \alpha V - \overline{\kappa} ) \sigma
-\int_\mathbb{R}2\frac{\partial^2 X^0}{\partial z^2} \nabla {r^1}^\mathrm{T}
 P\nabla r^0 \frac{\partial X^0}{\partial z}\,\mathrm{d}z \\
&-\int_\mathbb{R} 2 \sum_{l=1}^{2} \frac{\partial^2 X^0}{\partial z \partial s_l}
 \nabla {{s_l}^0}^\mathrm{T}P\nabla r^0 \frac{\partial X^0}{\partial z}\,\mathrm{d}z\\
&= 2\beta \vartheta^0.
\end{aligned}
\end{equation}
We collect the terms that would vanish in the derivation if $P=I$ as
\[
\text{additional}:=-\int_\mathbb{R}2\frac{\partial^2 X^0}{\partial z^2} 
\nabla {r^1}^\mathrm{T}P\nabla r^0 \frac{\partial X^0}{\partial z}\,\mathrm{d}z 
  -\int_\mathbb{R} 2 \sum_{l=1}^{2} 
\frac{\partial^2 X^0}{\partial z \partial s_l}\nabla {{s_l}^0}^\mathrm{T}
P\nabla r^0 \frac{\partial X^0}{\partial z}\,\mathrm{d}z.
\]
By recalling that $\vartheta^0$ does not depend on $z$ and by applying the 
matching condition \eqref{eq:first_matching_condition_theta}, we deduce
\begin{equation*}
\tilde{\theta}^0(t,\Gamma^0_{\pm})=\frac{(-\alpha V -\overline{\kappa} ) 
\sigma}{2\beta}+\frac{\text{additional}}{2\beta}
\end{equation*}
and, thus,
\begin{equation*}
 \theta^0(t,x) = \frac{(-\alpha V -\overline{\kappa} ) \sigma}{2\beta}
+\frac{\text{additional}}{2\beta} \quad\text{on}~\Gamma^0(t).
\end{equation*}
This relation and the identities \eqref{eq:sharp_matrix_heat_equation_renamed_o1}, 
\eqref{eq:sharp_heat_equation_renamed_o1_final}, and
 \eqref{eq:jump_condition_gamma0} now yield the sharp interface model
\begin{gather*}
\rho_{\mathrm{s}} {u^0}'-\operatorname{div}(P^{\mathrm{s}} \nabla u^0) 
+ \kappa_{RI}(u^0-\theta^0)=0 \quad\text{in }S\times\Omega, \\
\rho_{\mathrm{p}} {\theta^0}' - \operatorname{div}( P^{\mathrm{p}} \nabla \theta^0 ) 
+ \kappa_{RI}(\theta^0-u^0) = 0 \quad\text{for } t\in S,\; 
 x\in \Omega_l(t)\cup\Omega_s(t), \\
\nabla {r^0}^\mathrm{T}P^{\mathrm{p}}\nabla r^0 
\big[\frac{\partial \theta^0}{\partial \vec{n}}\big]_{\Gamma^0_\pm}
=-\lambda V \quad\text{for }t\in S,\; x\in\Gamma^0(t)
\end{gather*}
with
\begin{equation*}
 \theta^0=\frac{(-\alpha V -\overline{\kappa} ) \sigma}{2\beta}
+\frac{\text{additional}}{2\beta} \text{ for }t\in S,\; x\in\Gamma^0(t),
\end{equation*}
which is of similar structure as the Stefan problem with surface tension 
and kinetic undercooling, compare to \cite{Caginalp1989}.

\begin{thebibliography}{10}

\bibitem{Caginalp1989}
G. Caginalp;
 \emph{{S}tefan and {H}ele-{S}haw type models as asymptotic limits
 of the phase-field equations}, Phys.~Rev.~A (3) \textbf{39} (1989), no.~11,
 5887--5896.

\bibitem{Caginalp1998} G.~Caginalp, X.~Chen;
 \emph{Convergence of the phase field model to its
 sharp interface limits}, European J.~Appl.~Math. \textbf{9} (1998), 417--445.

\bibitem{Caginalp2005} G.~Caginalp, C.~Eck;
 \emph{Rapidly converging phase field models via second
 order asymptotics}, Discrete Contin.~Dyn.~Syst. \textbf{Supplement Volume
 2005} (2005), 142--152.

\bibitem{Caginalp1986} Gunduz Caginalp;
 \emph{An analysis of a phase field model of a free boundary},
 Arch.~Rational Mech.~Anal. \textbf{92} (1986), no.~3, 205--245 (English).

\bibitem{Eck2008} C.~Eck, H.~Garcke,  P.~Knabner;
 \emph{Mathematische Modellierung}, Springer,  2008.

\bibitem{Eck2004} Christof Eck;
 \emph{A two-scale phase field model for liquid-solid phase
 transitions of binary mixtures with dendritic microstructure}, Habilitation
 Thesis, Universit\"at Erlangen-N\"urnberg, July 2004.

\bibitem{Fife1995} P.~C. Fife, O.~Penrose;
\emph{Interfacial dynamics for thermodynamically
 consistent phase-field models with nonconserved order parameter},
 Electron.~J.~Differential Equations \textbf{1995} (1995), no.~16, 1--49.

\bibitem{Hoepker2016} M.~H\"opker;
 \emph{Extension operators for {S}obolev spaces on periodic
 domains, their applications, and homogenization of a phase field model for
 phase transitions in porous media}, Ph.D. thesis, Universit\"at Bremen, 2016.

\bibitem{Meirmanov1992} Anvarbek M. Meirmanov;
 \emph{The {S}tefan problem}, De Gruyter expositions in
 mathematics, De Gruyter, 1992.

\bibitem{Showalter1982} R.~E. Showalter;
 \emph{Mathematical formulation of the {S}tefan problem}, Int.
 J. Engng Sci. \textbf{20} (1982), no.~8, 909--912.

\bibitem{Stefan1891} J. Stefan;
 \emph{Ueber die {T}heorie der {E}isbildung, insbesondere \"uber die
 {E}isbildung im {P}olarmeere}, Ann.~Phys. \textbf{278} (1891), no.~2,
 269--286.

\bibitem{Visintin1996} Augusto Visintin;
 \emph{Models of phase transitions}, Progress in Nonlinear
 Differential Equations, vol.~28, Birkh\"auser, 1996.

\end{thebibliography}

\end{document}