Local existence of strong solutions to micro–macro models for reactive transport in evolving porous media

Abstract

Two-scale models pose a promising approach in simulating reactive flow and transport in evolving porous media. Classically, homogenised flow and transport equations are solved on the macroscopic scale, while effective parameters are obtained from auxiliary cell problems on possibly evolving reference geometries (micro-scale). Despite their perspective success in rendering lab/field-scale simulations computationally feasible, analytic results regarding the arising two-scale bilaterally coupled system often restrict to simplified models. In this paper, we first derive smooth dependence results concerning the partial coupling from the underlying geometry to macroscopic quantities. Therefore, alterations of the representative fluid domain are described by smooth paths of diffeomorphisms. Exploiting the gained regularity of the effective space- and time-dependent macroscopic coefficients, we present local-in-time existence results for strong solutions to the partially coupled micro–macro system using fixed-point arguments. What is more, we extend our results to the bilaterally coupled diffusive transport model including a level-set description of the evolving geometry.

1. Introduction

Multi-scale techniques are known as a powerful tool to model reactive transport in porous media [27]. They enjoyed large interest over the last decades due to their wide range of applicability in groundwater treatment, bioremediation, soil formation due to precipitation and dissolution of minerals, or weathering as well as technical/engineering applications and geothermal energy systems [6, 37]. For many of the aforementioned applications, the structure of the porous medium is altered. Therefore, such problems inherently belong to the class of free-boundary problems, i.e. moving boundary problems with a-priori unknown interface evolution. Moving/free-boundary problems enjoy continuous mathematical concern and development [6, 8]. An overview of methods to analyse moving boundary problems such as energy methods or tools from operator and semi-group theory is found in [33]. Multi-scale models are typically derived by means of volume averaging [45] or homogenisation [19] and can be extended to the situation of evolving microstructures [27] in the context of reactive transport.

In two-scale models for reactive transport in evolving porous media such as derived in [5, 44] flow and reactive transport equations are typically solved on the macroscopic domain. These Partial Differential Equations (PDEs) encompass several effective parameters as coefficients such as porosity or diffusivity that are connected to the underlying microscopic geometry. Due to the evolution of the porous medium, the arising parameters depend on both space and time. As this evolution is for instance often driven by chemical reactions, i.e. dependent on the solution of the transport equation, the type of models considered here inherently features a two-way coupling between the scales complicating analytical treatment. We illustrate the coupling of the macroscopic equations to the underlying geometry resolved in representative unit cells in Figure 1. Henceforth, we distinguish two different types of coupling. The two-way coupling between both scales will be referred to as full coupling. Commonly, also a simplified coupling structure is investigated in the literature, disregarding the back-coupling from the macro to the micro-scale, cf. Figure 1. We refer to the arising one-sided coupling as a partial coupling scenario.

Figure 1. Schematic presentation of micro–macro coupling in multi-scale reactive transport models. Geometry-dependent effective parameters influence the macroscopic flow and solute transport, see black-coloured arrows. In the fully coupled scenario, solute concentrations on the macroscopic domain $\Omega$ additionally prescribe the evolution of the underlying microscopic geometry, see red-coloured arrows.

An additional challenge is posed by a suitable framework to capture the evolving geometries. A variety of convenient methods is outlined in [1] regarding the description of evolving microstructures. Commonly, level-set methods or phase-field approaches are used, in which case the macroscopic concentrations either prescribe a normal interface velocity, cf. [44], or induce a source term to the phase field [5], respectively. However, especially for highly symmetrical shapes, geometry evolution is often modelled and simulated in a simplified way via Ordinary Differential Equations (ODEs) for typical characteristic parameters such as the porosity in case of [35, 36] or the thickness of the precipitation layer as in [43]. Furthermore, models can also handle the geometry evolution implicitly by computing porosity from conservation of mass and deducing all geometry-related effective quantities therefrom by heuristic laws, cf. [20, 41] facilitating the numerical treatment. In [14, 30] micro–macro models including fully resolved microscopic geometries have been investigated numerically. Regarding the most general setting, the associated evolution PDE (level-set or phase-field equation) is solved on reference geometries virtually attached to each degree of freedom of the macroscopic discretisation. We note that the analysis of upscaled diffuse interface (phase-field) models is relatively straightforward as the phase-field variable describing the geometry directly enters the governing equations as coefficients. However, the direct analysis typically involves estimates which do not hold in the corresponding sharp-interface limit, i.e. blow up for diffusive interface width tending to zero, cf. [30] requiring an additional examination of the limit model.

In the literature, different approaches are present to obtain existence results for the effective reactive transport model including geometry alterations. In [12, 32] perforated microscopic domains are mapped onto a periodic reference domain using diffeomorphisms. Existence results for the transformed microscopic equations on the reference domain are then leveraged to the effective model by means of two-scale convergence. Typically, the upscaling process of the transformed model is of increased complexity due to the appearance of additional factors in the highest-order terms arising from the transformation itself. Restricting to diffusion–reaction systems, [31] presents existence results to the fully coupled system with geometry evolution modelled via an ODE for the determinant of the deformation gradient. As such, only effects emerging from changes in the pore-space volume are reflected. Likewise, diffusion–advection–reaction equations are treated in [11] assuming an a-priori given geometry evolution (partial coupling).

On the other hand, existence results can also be derived by investigating the effective model itself. Existence of weak solutions to a homogenised diffusion-driven model with partial and full coupling between the scales was shown in [24] using transformations of the model equations onto fixed reference domains. However, the analysis performed requires the smoothness of the effective diffusion tensor as an additional assumption (partial coupling) or neglects its evolution with time completely (full coupling). Likewise, in [34], fully coupled systems with diffusion-driven transport are investigated. Furthermore, [33] considered partially coupled systems with advection including degenerating hydrodynamic parameters. Yet, a key assumption is the a-priori knowledge of the relation between effective parameters and the porosity which plays the role of an order parameter.

In this paper, we follow the latter approach performing analysis directly on the effective micro–macro model. As a key result, we prove smooth dependence of effective parameters on the underlying geometry. More precisely, we consider smoothly bounded underlying microscopic geometries and restrict to setups precluding degeneracy of effective parameters. By describing geometry alterations by smoothly parameterised paths of diffeomorphisms, we make use of an additional parameter characterising its state. As such, the presented framework covers a broad range of geometry evolution which does not rely on parameterisability by a single physical quantity such as the porosity. On the one hand, we use these results to investigate the existence of strong local-in-time solutions to the partially coupled model potentially including advective solute transport. On the other hand, we treat the scenario of full coupling and diffusion-driven transport. In that case, the macroscopic concentrations are coupled to the full level-set equation for geometry evolution. Thereby, we extend similar results known in the literature for simple and restrictive geometries like cubes or spheres, cf. [34]. Moreover, we discuss sufficient conditions for the long-term existence of solutions in the presented framework.

Our paper is outlined as follows: in Section 2, we present an established model for reactive flow and transport in evolving porous media. Restricting to solute transport by diffusion only in Section 3, we derive smooth dependence results for the diffusion tensor on the geometry and prove local-in-time existence to the partially and fully coupled model. Establishing analogous results for the permeability tensors in Section 4, smooth dependence of the Darcy velocity field on the underlying geometry is shown. Finally, existence of solutions to the partially coupled model including advective solute transport is proven in Section 4.4.

2. Model

This research is based on a micro–macro model for reactive flow and transport in evolving porous media introduced by [42] where its derivation from a detailed pore-scale model by formal homogenisation arguments was performed in two spatial dimensions. A generalisation to three dimensions was derived in [34] by modification of the original deduction. More precisely, the model under consideration consists of a set of coupled PDEs describing solute transport, fluid flow and geometry evolution. First, we introduce the transport equation for a solute chemical species $c$ on the macroscopic domain of interest $\Omega \subset \mathbb{R}^d$ with $d\in \{2,3\}$ :

(2.1) \begin{align} \partial _t (\phi c) + \nabla \cdot ({v} c) - \nabla \cdot (\mathbb{D}\nabla c) = \sigma f(c) \quad \text{in } (0,T)\times \Omega, \end{align}

where $\phi$ is the porosity, $\sigma$ is the specific surface area, $v$ is the flow velocity, $\mathbb{D}$ is the effective diffusion tensor and $f$ is a source/sink term due to heterogeneous reactions. As such, the parabolic equation (2.1) models solute transport by diffusion and advection processes as well as chemical reactions at the fluid–solid interface. In this paper, we consider equation (2.1) with Dirichlet boundary conditions as well as initial conditions of the form:

\begin{align*} c(t,x) &= C^0(t,x) && t\in (0,T), \; x\in \partial \Omega, \\ c(0,x) &= c_0(x) && x\in \Omega, \end{align*}

cf. Theorem 3, 6.

The effective parameters in (2.1) are derived from unit cells $Y= \left ({-}\frac{1}{2}, \frac{1}{2} \right )^d$ representing a local reference elementary volume of the underlying geometry. The respective exterior boundary is denoted by $\partial Y$ . For the following, we consider solid inclusions compactly contained within $Y$ . Let us denote the remaining fluid domain therein by $\mathcal{P}$ and the interior boundary by $\partial ^{\text{int}}\mathcal{P}$ , cf. Figure 1. The diffusion tensor is then given as

(2.2) \begin{align} \mathbb{D}_{i,j}(t,x) \,:\!=\, \int _{\mathcal{P}(t,x)}\left (\partial _{y_i}\zeta _j + \delta _{ij}\right )\,dy, \quad i,j\in \lbrace 1,\dots,d\rbrace, \end{align}

with Kronecker delta $\delta _{ij}$ , where $\zeta _j$ are the solutions to the following elliptic problems:

(2.3) \begin{align} - \nabla _y\cdot (\nabla _y \zeta _j) &= 0 && \text{in }{\mathcal{P}(t,x)}, \nonumber \\ \nabla _y \zeta _j\cdot \nu &= - e_j\cdot \nu && \text{on }{\partial ^{\text{int}}\mathcal{P}(t,x)},\\ \zeta _j \text{ periodic in } y, &\quad \int \limits _{\mathcal{P}(t,x)} \zeta _j \; dy = 0, \nonumber \end{align}

with outer unit normal $\nu =\nu (t,x,y)$ . For the derivation in the context of evolving geometries, see [42]. In this case, the evolution is reflected in time-dependent domains. Note that the shape of the elliptic problem (2.3) is identical to the one derived under the assumption of fixed underlying geometries in [18], leading to fixed domains $\mathcal{P}$ . We refer to (2.3) as diffusion cell problems.

Second, the advective flow field $v$ and the associated pressure field $p$ are given by Darcy’s equation:

(2.4) \begin{align} v &= -\frac{\mathbb{K}}{\mu } \nabla p &&\text{in } \Omega,\; t\in (0,T), \\ \nabla \cdot v &=0 && \text{in } \Omega,\; t\in (0,T), \nonumber \end{align}

with viscosity $\mu$ and permeability tensor $\mathbb{K}$ . As $\mu$ is a constant being characteristic to the solvent, we set it to one for convenience. Darcy’s equation is supplemented by the following boundary condition:

\begin{align*} p(t,x) &= p_0(t,x) && t\in (0,T), \; x \in \partial \Omega _{\text{Dir}}, \\ v(t,x) &= v_0(t,x) && t\in (0,T), \; x \in \partial \Omega _{\text{flux}}, \end{align*}

with a partition $\partial \Omega _{\text{Dir}} \dot{\cup } \partial \Omega _{\text{flux}} = \partial \Omega$ . For uniqueness of the pressure solution, a positive measure of $\partial \Omega _{\text{Dir}}$ is required.

Note that the condition of a divergence-free velocity field in (2.4) is a common simplification as discussed in [14]. Due to the much larger time-scale of geometry evolution compared to fluid flow, it is justified to disregard the flow induced by fluid displacement arising from a variable pore-space volume.

The permeability tensor in (2.4) is defined as

(2.5) \begin{align} \mathbb{K}_{i,j}(t,x) \,:\!=\, \int _{\mathcal{P}(t,x)} \omega _j^i \; dy, \quad i,j\in \lbrace 1,\dots,d\rbrace, \end{align}

where $(\omega _j, \pi _j)$ are the solutions to the Stokes-type problems, cf. [42]:

(2.6) \begin{align} -\Delta _y \omega _j + \nabla _y \pi _j &= e_j && \text{in }{\mathcal{P}(t,x)},\nonumber \\ \nabla _y \cdot \omega _j &=0 && \text{in }{\mathcal{P}(t,x)},\\ \omega _j &=0 &&\text{on }{\partial ^{\text{int}}\mathcal{P}(t,x)}, \nonumber \\ \omega _j, \; \pi _j \text{ periodic in } y, &\quad \int \limits _{\mathcal{P}(t,x)}\pi _j \; dy = 0. \nonumber \end{align}

Likewise, the shape of the Stokes-type problem is identical to the one derived under the assumption of fixed underlying geometries in [18]. Again, the evolving geometry is reflected in the time-dependent domains. We refer to (2.6) as permeability cell problems. Note that both cell problems (2.3), (2.6) result in symmetric positive semi-definite tensors $\mathbb{D}$ , $\mathbb{K}$ .

Finally, the level-set method is used to describe the evolving underlying geometries [37]. Therefore, we assume the existence of a level-set function $\Phi _0\,:\,\Omega \times Y \to \mathbb{R}$ characterising the solid part within the unit-cell attached to the macroscopic point $x\in \Omega$ at initial time by $\{\Phi _0(x,\cdot )\gt 0\}$ . Consequently, $\{\Phi _0(x,\cdot )\lt 0\}$ refers to the fluid domain $\mathcal{P}$ and $\{\Phi _0(x,\cdot )=0\}$ denotes the fluid–solid interface $\partial ^{\text{int}}\mathcal{P}$ . We require that the gradient of $\Phi _0$ does not vanish along the zero-level-set to ensure the representation of a submanifold of codimension one. For a normal interface velocity field $v_n\,:\,(0, T)\times \Omega \times Y \to \mathbb{R}$ , the evolution of $\Phi _0$ is described by the level-set equation for $\Phi \,:\,(0, T)\times \Omega \times Y \to \mathbb{R}$ , cf. [37]:

(2.7) \begin{align} \frac{\partial \Phi }{\partial t} + v_n |\nabla _y \Phi | &= 0 &&\quad \text{in } (0, T)\times \Omega \times Y, \\ \Phi (0,\cdot,\cdot ) &= \Phi _0 &&\quad \text{in } \Omega \times Y.\nonumber \end{align}

Due to the finite propagation speed of information in (2.7), cf. [10], we refrain from posing explicit boundary conditions on $(0, T)\times \Omega \times \partial Y$ here, see also Section 3.4. The different sub-domains of $Y$ (fluid domain, solid domain, interface) at a certain time $t$ are encoded by the sign of $\Phi (t,\cdot,\cdot )$ according to the convention for the initial condition $\Phi _0$ above. As such, the level-set function uniquely determines porosity $\phi$ and specific surface $\sigma$ as required in (2.1) via

(2.8) \begin{align} \phi (t,x) = \mathcal{H}^{d}(\{\Phi (t,x,\cdot ) \lt 0\}), \quad \sigma (t,x) = \mathcal{H}^{d-1}(\{\Phi (t,x,\cdot ) =0\}), \end{align}

using the d-dimensional Hausdorff measure $\mathcal{H}^{d}$ . The interface velocity is coupled to the chemical reaction in a mass conserving way, e.g.

(2.9) \begin{align} v_n(t,x,y) =- v_{\text{mod}}(y) f(c(t,x)), \end{align}

potentially using a scalar speed modification function $v_{\text{mod}}$ which allows for a varying normal interface velocity within a unit-cell, cf. [13]. Depending on the sign of $f$ , the model allows for shrinking as well as growing interfaces. Necessary restrictions on the shape of $f$ are discussed in Sections 3.3, 3.5.

3. Smooth parameter dependence and existence for diffusive transport

In this section, we consider a special case of the model introduced in Section 2. Neglecting advective transport for the solute species, we consider the partially and fully coupled case, cf. Figure 1. As such, the coupling from the micro to the macro-scale is conveyed by $\phi,\sigma$ and $\mathbb{D}$ only, facilitating the analysis. After discussing the setup in more detail, this section first considers the smoothness of the partial coupling. Therefore, we investigate the dependence of the diffusion tensor $\mathbb{D}$ on deformations of the microscopic geometry via diffeomorphisms in Section 3.2 as the principle step to establish existence results for the partially coupled problem in Section 3.3. Moreover, we show the induction of suitable diffeomorphisms by the level-set equation in Section 3.4. By doing so, we establish regular dependence of the required effective quantities $\phi,\sigma,\mathbb{D}$ on the level-set equation’s time variable, enabling the application of fixed-point arguments. This ultimately leads to local-in-time existence results for the fully coupled case, cf. Figure 1, in Section 3.5. The roadmap for this section as described above is illustrated in Figure 2.

Figure 2. Roadmap for Section 3: Key steps to perform fixed-point argument on existence of solutions of the model specified in Section 3.1. Main statements are given in blue boxes, linking lemmas, theorems and corollaries are highlighted in red boxes. Black arrows indicate the steps necessary in the partial coupling scenario, red arrows refer to additional steps necessary for the fully coupled case, cf. Figure 1.

3.1. Setting

For the following, we consider solute transport by diffusion only. As the term $\partial _t (\phi c)$ in (2.1) is difficult to handle analytically, it is shifted to the right-hand side, cf. [34]. Therefore, we write

(3.1) \begin{align} \phi \partial _t c - \nabla \cdot (\mathbb{D}\nabla c) = \sigma f(c) - \partial _t \phi c \quad \text{in } (0,T)\times \Omega, \end{align}

rendering the equations for flow and permeability determination (2.4), (2.6) superfluous. Hence, it is sufficient to close the model regarded in this section by (2.3), (2.7) and (2.9). We emphasise that the following considerations assume good-natured conditions such as the smoothly bounded solid geometry being compactly contained within the unit-cell $Y$ . By considering local-in-time estimates, this setting is maintained by an appropriate choice of initial conditions. Using diffeomorphisms to describe solid alteration, we particularly exclude clogging scenarios and degenerating equations.

3.2. Continuous dependence of diffusion tensors

In order to prove existence to the model described in Section 3.1, we make extensive use of existence theory for linear parabolic equations, cf. Theorem A.3 in the appendix. As such, we require moderate regularity for $\mathbb{D}$ as the coefficient of the leading order term in (3.1).

Therefore, this section is concerned with the dependence of the diffusion tensor $\mathbb{D}$ on the evolving geometry. To do so, we first focus on a single unit cell and therefore on a single cell problem posed on the fluid domain $\mathcal{P}$ . While regarding $(t,x)$ as independent parameters for now, we present in Section 3.3 how regularity with respect to the macroscopic variables is established. The method presented consists of three steps: At first, we establish higher regularity for weak solutions to the diffusion cell problem (2.3). Although using standard methods, we state the results in detail due to the uncommon periodic boundary conditions. Based on that, a mapping between the geometry and the elliptic PDE’s solution of desired regularity is constructed using the implicit function theorem following the technique of [16]. Finally, we extend our smoothness results from the solutions $\zeta _j$ , $j\in \{1,\dots,d\}$ , of (2.3) to the diffusion tensor $\mathbb{D}$ which is given as an affine-linear functional of $\zeta _j$ , cf. (2.2).

For the formulation of problem (2.3) and our regularity result Lemma 3.1, we assume $\partial ^{\text{int}}\mathcal{P}$ to be $C^{2,1}$ -regular. As it becomes apparent in Theorem 1, it is necessary to consider the full class of elliptic PDEs of type (2.3) with general source term and Neumann boundary conditions. In a first step towards a suitable weak formulation of problem (2.3), we introduce periodic Sobolev spaces according to [8]. Let $\mathcal{P}_{\text{ext}} \subset \mathbb{R}^d$ denote the perforated domain obtained by periodic extension of $\mathcal{P}$ in $\mathbb{R}^d$ . Then we define $H_\#^k (\mathcal{P})$ for $k\in \mathbb{N}$ as the closure of Y-periodic functions in $C^\infty (\mathcal{P}_{\text{ext}})$ with respect to the $H^k$ -norm. For the weak formulation of the diffusion cell problem with Neumann boundary conditions, we introduce the following function spaces:

\begin{align*} H^k_{\#,0}(\mathcal{P}) = \left \lbrace v\in H_\#^k ( \mathcal{P})\,:\, \int _{\mathcal{P}} v \; dy =0 \right \rbrace, \end{align*}

equipped with the norm $||u||_{ H^k_{\#,0}(\mathcal{P})} = ||u_{|\mathcal{P}}||_{ H^k(\mathcal{P})}$ .

Accordingly, the weak problem with general source term $f$ and Neumann boundary condition $g$ reads: Find $u\in H^1_{\#,0}(\mathcal{P})$ such that

(3.2) \begin{align} \int _{\mathcal{P}} \nabla u \cdot \nabla v \; dy - \int _{\partial ^{\text{int}}\mathcal{P}} g v \; d\sigma = \int _{\mathcal{P}} f v \; dy, \quad \forall v\in H^1_{\#}(\mathcal{P}), \end{align}

with $g\in L^2(\partial ^{\text{int}}\mathcal{P})$ , $f\in L^2(\mathcal{P})$ . Apparently, the compatibility condition $\int \limits _{\mathcal{P}} f \; dx =- \int \limits _{\partial ^{\text{int}}\mathcal{P}} g \; d\sigma$ is necessary for solvability.

Next, we establish unique solvability for the above problem and present conditions ensuring solutions to be of higher regularity. Particularly, a $C^{2,1}$ -regular interior boundary $\partial ^{\text{int}}\mathcal{P}$ proves to be sufficient for equations (2.3) to hold in the $L^2$ -sense and therefore point-wise almost everywhere in $\mathcal{P}$ .

Lemma 3.1. Elliptic regularity. Problem ( 3.2 ) as the weak generalisation of ( 2.3 ) has a unique solution $u\in H^1_{\#,0}(\mathcal{P})$ . Let $k$ be an integer number $k\geq 0$ . If the interior boundary $\partial ^{\text{int}}\mathcal{P}$ is furthermore $C^{k+2,1}$ -regular and $g\in H^{k+\frac{1}{2}}(\partial ^{\text{int}}\mathcal{P})$ , $f\in H_\#^k(\mathcal{P})$ fulfilling the compatibility condition $\int \limits _{\mathcal{P}} f \; dx = -\int \limits _{\partial ^{\text{int}}\mathcal{P}} g \; d\sigma$ , then $u\in H^{k+2}_{\#,0}(\mathcal{P})$ .

Proof. The coercivity of the bilinear form on the Hilbert space $H^1_{\#,0}(\mathcal{P})$ is guaranteed by Poincaré’s inequality. The Lax–Milgram theorem therefore implies the unique existence of a solution $u \in H^1_{\#,0}(\mathcal{P})$ to the weak formulation, cf. [8].

Let us now suppose that $\partial ^{\text{int}}\mathcal{P}$ is $C^{k+2,1}$ -regular. By the surjectivity of the higher-order trace operator continuously extending the following function, cf. [27],

\begin{align*}{\text{Tr}}_k &\,:\, H^{k+2}(\mathcal{P}) \to \prod _{l=0}^{k+1} H^{k-l+\frac{3}{2}} (\partial ^{\text{int}}\mathcal{P}), \quad \forall u \in H^{k+2}(\mathcal{P})\cap C^{k+1}(\bar{\mathcal{P}}), \\{\text{Tr}}_k u &= \left (u|_{\partial ^{\text{int}}\mathcal{P}}, \partial _\nu u|_{\partial ^{\text{int}}\mathcal{P}},\cdots, \partial _\nu ^{k+1} u|_{\partial ^{\text{int}}\mathcal{P}} \right ), \nonumber \end{align*}

we find a function $\Psi \in H^{k+2}_{\#,0}(\mathcal{P})$ with $\dfrac{\partial \Psi }{\partial \nu } = \dfrac{\partial u}{\partial \nu }$ on $\partial ^{\text{int}}\mathcal{P}$ in the trace sense vanishing in a neighbourhood of $\partial Y$ . Exploiting linearity of the problem and carrying out the subsequent argument for $u-\Psi$ and the corresponding source term $\tilde{f}=f-\Delta \Psi \in H^k(\mathcal{P})$ , we can assume homogeneous Neumann boundary conditions. By Theorem 3 of [25], the required interior higher regularity is established. Following the lines of Theorem 4 of [25], higher regularity is also obtained in a neighbourhood of every interior boundary point $y\in \partial ^{\text{int}}\mathcal{P}$ . Using an open covering argument, we obtain $u \in H^{k+2}_{\#,0}(\mathcal{P})$ .

Remark 3.2. Neumann boundary conditions. As $\partial ^{\text{int}}\mathcal{P}$ is of class $C^{2,1}$ , $\nu \cdot e_1 \in H^{\frac{1}{2}}(\partial ^{\text{int}}\mathcal{P})$ holds. Furthermore, we have $\int _{\partial ^{\text{int}}\mathcal{P}} \nu \cdot e_1 \; d\sigma = 0$ by Gauss’s theorem. As such, the last lemma covers the unique solvability of problem ( 2.3 ) in $H^{2}_{\#,0}(\mathcal{P})$ and the equation holds in a point-wise almost-everywhere sense.

Remark 3.3. Higher regularity for homogenised model. Note that the higher regularity result of Lemma 3.1 holds for the diffusion cell problems of the homogenised model. However, regarding the upscaling process from the associated pore-scale diffusion equation to the effective one considered here, the convergence of the sequence of transport problems $c_\varepsilon$ defined on the $\varepsilon$ -periodic domains $\Omega _\varepsilon$ against the homogenised solution $c$ is not valid with respect to these stronger norms. This is essentially due to the absence of uniform boundedness of $c_\varepsilon$ with respect to $\varepsilon$ .

Remark 3.2 enables us to define mappings from the solution space of equation (3.2) to the bulk and boundary data spaces with a point-wise interpretation. More precisely, we consider such mappings which involve the geometry alteration as a parameter. This idea poses the main ingredient in the subsequent investigation of the dependence of solutions to the weak problem (3.2) and functionals thereof on the geometry $\mathcal{P}$ . Following the approach presented by [16], a Lagrangian description of the initial problem on varying domains is taken. The main step is to rewrite the equation on a fixed domain of reference $\mathcal{P}$ . This technique has also been successfully applied to the homogenisation of PDEs on non-uniformly periodic or evolving domains in the context of porous media [11, 17, 30] or shape optimisation minimising energy functionals depending on a PDE’s solution [39]. In order to re-define functions mapping from the altered domains as functions on a fixed domain of reference, we make use of the concept of diffeomorphisms and pullbacks:

Definition 1. Diffeomorphism. Let $h\,:\,Y \to h(Y)\subset \mathbb{R}^d$ be a bijective mapping with $h\in C^{k,\alpha }(Y, \mathbb{R}^d)$ for $k\geq 1,$ $\alpha \in [0,1]$ . We call $h$ a diffeomorphism of class ${\text{Diff}}^{k,\alpha }(Y, \mathbb{R}^d)$ iff the inverse $h^{-1}$ satisfies $h^{-1}\in C^{k,\alpha }(h(Y),Y)$ .

The action of a diffeomorphism on the set $\mathcal{P}\subset Y$ is illustrated in Figure 3.

Figure 3. Smooth deformation of domain $\mathcal{P}$ with circular inclusion mediated by a diffeomorphism $h\in{\text{Diff}}_\Box ^{1} (\bar{Y})$ . As $h$ preserves the exterior boundary $\partial Y$ , the image-set $h(\mathcal{P})$ is an admissible unit cell pore-space geometry.

Definition 2. Pullback. Let $h$ be a diffeomorphism $h\in{\text{Diff}}^{1}(Y)\subset C^1(Y, \mathbb{R}^d)$ and $l\,:\,h(Y) \to \mathbb{R}^n$ , $n\in \mathbb{N}$ . We define the pullback $h^*$ by

\begin{align*} h^*(l)\,:\,Y \to \mathbb{R}^n, \quad h^* l(x) \,:\!=\, l(h(x)). \end{align*}

Note that the pullback is a linear operation and its inverse is given as $(h^*)^{-1}(l) = (h^{-1})^*(l)$ . Furthermore, for sufficiently smooth diffeomorphisms, the pullback is a bounded operator between the $H^m$ function spaces on the original and deformed set, cf. [16, 23]. In Example 3.5, we present an explicit construction of a family of diffeomorphisms mapping circular inclusion of different radii to one another and provide illustrations of the associated pullback of a distance function in Figure 4.

Figure 4. Visualisation of diffeomorphism (3.7) for $r_2=0.3,\; r_1=0.1$ : (a) illustrates the related circles ( $r_2$ -black, $r_1$ -red) posing the interior boundary of the domain. (b) The displacement field is shown, i.e. $h_{r_1}-\text{id}_Y$ . As enforced by the interpolation function $\xi$ in (3.7), the displacement smoothly vanishes close to the exterior boundary $\partial Y$ . (c), The graph of $h_{r_1}$ , cf. (3.7), is shown along the $y_1$ axis illustrating the three different sections (uniform contraction, transition and identity). (d) Displays the pullback $h_{r_1}^*(\Phi )$ with $\Phi (y_1,y_2)=r_1-||(y_1,y_2)||_2$ . Contour lines uniformly spaced in increments of $0.1$ are added in white. The zero-level-set of $h_{r_1}^*(\Phi )$ highlighted in red corresponds to a circle of radius $r_2$ .

Using the tools introduced in Definitions 1 and 2, we are now able to characterise solutions to (3.2) on the domain $h(\mathcal{P})$ as roots of a function $F$ . Taking a Lagrangian point of view, we work on a fixed domain of reference $\mathcal{P}$ . As such, functions first need to be conveyed to the deformed domain $h(\mathcal{P})$ where the differential operators according to (2.3) are applied. Performing a pullback with the inverse deformation, functions are translated back onto the domain of reference. Accordingly, let us consider the following mapping:

(3.3)

where $\nu _{h(\mathcal{P})}$ denotes the outer unit normal with respect to the domain $h(\mathcal{P})$ and ${\text{Tr}}_{h(\mathcal{P})}$ the standard trace operator on $h(\mathcal{P})$ . Note that the normal vectors can be extended within a tubular neighbourhood of $\partial ^{\text{int}}\mathcal{P}$ , cf. Theorem A.1. By the trace theorem and change of variable rule, $F$ is well-defined as a mapping between the stated spaces. Note that in the first component $F_1$ a vanishing mean value is enforced. Thereby, we eliminate an additional degree of freedom to ensure surjectivity of $F$ .

Since the image of the unit-cell under an arbitrary diffeomorphism may not be a unit cell, we must also introduce a restricted class of deformations. Therefore, let ${\text{Diff}}_\Box (\bar{Y})$ denote the set of diffeomorphism preserving the exterior boundary $\partial Y$ defined by:

(3.4) \begin{align} {\text{Diff}}_\Box ^{k,\alpha } (\bar{Y}) = \left \lbrace h \in{\text{Diff}}^{k,\alpha } (\bar{Y},\bar{Y})\,:\, h_{|U} = \text{id}_U \right \rbrace, \end{align}

for some fixed neighbourhood $\partial Y \subset U\subset \bar{Y}$ . In fact, the diffeomorphism illustrated in Figure 3 belongs to this specified class. The construction (3.4) is inspired by the Hanzawa transformation, cf. [32], and requires the diffeomorphism to decay smoothly towards the identity at the exterior boundary. As such, periodic functions with respect to $Y$ admit a periodic pullback for $h\in{\text{Diff}}^{k,\alpha }_\Box (\bar{Y})$ . The importance of the mapping $F$ defined in (3.3) is now reflected in the following characterisation property. For $u\in H^2_{\#,0}(\mathcal{P})$ , it holds

for all $h\in{\text{Diff}}^{2,1}_\Box (\bar{Y})$ . As such, we established a one-on-one relation between the roots of $F$ and the PDE’s solution on the deformed domain $h(\mathcal{P})$ .

We will now apply the implicit function theorem to $F$ to obtain a continuous mapping $h \mapsto u$ in the respective spaces as summarised in Theorem 1 below. As a result, a relation between the deformation of $\mathcal{P}$ and the associated solution to (2.3) is established. Therefore, we check the following properties:

Continuity of F: Let us denote the pullback of $u$ via $h^{-1}$ by $v=h^{*-1}u$ constituting a function on $h(\bar{Y})$ and the image point on $h(\bar{Y})$ by $y=h(x)$ . First, we consider the first term of the first component $F_1$ of the mapping $F$ . By applying the chain rule, we have

(3.5) \begin{align} & h^*\Delta h^{*-1} u (x) = \Delta v(y) = \Delta u(h^{-1}(y)) = \\ &\sum _{i,j=1}^d \frac{\partial ^2 u}{\partial x_i \partial x_j} \left (\sum _{k=1}^d \frac{\partial h^{-1}_i}{\partial y_k} \frac{\partial h^{-1}_j}{\partial y_k} \right ) +\sum _{i=1}^d \frac{\partial u}{\partial x_i} \left (\sum _{k=1}^d \frac{\partial ^2 h^{-1}_i}{\partial y_k \partial y_k} \right ). \nonumber \end{align}

Note that by the inverse function theorem and the representation of a matrix’ inverse via the cofactor matrix, we can rewrite all partial derivatives of $h_i^{-1}(y)$ as a function of derivatives of $h_i(x)$ only involving multiplications and division by $\text{det}(\nabla h(x))$ . By the uniform boundedness of the last expression away from zero, this map is in particular locally Lipschitz continuous. For any sequence $(u_i,h_i)_{i \in \mathbb{N}}$ in $ H^2_{\#,0}(\mathcal{P}) \times{\text{Diff}}^{2,1}(\bar{Y})$ converging to $(u,h)$ in the product topology, we have the convergence of the derivatives of $u_i$ in $L^2$ and the convergence of the derivatives of $h_i$ uniformly. Hence, expression (3.5) is continuous in $(u,h)\in H^2_{\#,0}(\mathcal{P}) \times{\text{Diff}}^{2,1}(\bar{Y})$ . As the integral is a linear and bounded operator $\int \,:\, L^2(\mathcal{P})\to \mathbb{R}$ , we obtain continuity for $F_1$ .

We can apply the same strategy to prove continuity of the second component $F_2$ of $F$ . Noting the representation

\begin{align*} \nu _{h(\mathcal{P})} (y) = (\nabla h)^{-T}(x) \; \nu _{\mathcal{P}}(x) ||(\nabla h)^{-T}(x) \; \nu _{\mathcal{P}}(x)||^{-1}_2 \end{align*}

derived in [16] using the inverse transposed Jacobian matrix $(\nabla h)^{-T}$ , we compute

\begin{align*} F_{2}(u,h) (x) &= h^{*}{\text{Tr}}_{h(\mathcal{P})}\left ((\nu _{h(\mathcal{P})} \cdot \nabla h^{*-1}u) (y) - \nu _{h(\mathcal{P})}(y) \cdot e_1\right ) \nonumber \\ &={\text{Tr}}_{\mathcal{P}} \left (\left ( (\nabla h)^{-T} \nu _{\mathcal{P}} \right ) \cdot \left [\sum _{i=1}^d\frac{\partial u}{\partial x_i} \frac{\partial h^{-1}_i}{\partial y_j}-\delta _{j,1}\right ]_j \cdot ||(\nabla h)^{-T}(x) \; \nu _{\mathcal{P}}(x)||_2^{-1}\right ) \end{align*}

for $x\in \partial ^{\text{int}}\mathcal{P}$ . As ${\text{Diff}}^{2,1}(\bar{Y}) \subset C^{2,1}(\bar{Y},\mathbb{R}^d)$ is open, we conclude the continuity of $F$ on a neighbourhood $V$ of $(u,\text{id}_{\bar{Y}})$ .

Continuity of $F^{\prime }_u$ : By Definition 2 the pullback operator is linear. Using the linearity of the differential and trace operators involved, we conclude

for all $w\in H^2_{\#,0}(\mathcal{P})$ . Following the arguments from above, we obtain continuous Fréchet differentiability with respect to the first argument on a neighbourhood $V$ of $(u,\text{id}_{\bar{Y}})$ .

Bijectivity of $F^{\prime }_u$ : The bijectivity of $F^{\prime }_u(u,\text{id}_{\bar{Y}})$ onto $L_0^2(\mathcal{P}) \times H^{\frac{1}{2}}(\partial ^{\text{int}}\mathcal{P})$ is equivalent to finding a unique solution to the elliptic problem (2.3) on $\mathcal{P}$ . More precisely, for a given point $(f,g)$ in the image space of $F$ , we search for a weak solution $u$ for the Neumann boundary condition $g$ and source term $f+c$ for a constant $c\in \mathbb{R}$ . By Lemma 3.1, there exists exactly one $c$ such that the problem admits a solution (compatibility condition) in $H^2_{\#,0}(\mathcal{P})$ . In that case, the solution is unique. Note that due to the set of invertible bounded linear operators between Banach spaces being open and the continuity of $F^{\prime }_u$ , the bijectivity property of $F^{\prime }_u$ in fact holds on a neighbourhood of $(u,\text{id}_{\bar{Y}})$ .

Summarising the above arguments, we conclude the following statement.

Theorem 1. Continuous dependence of weak solutions to ( 2.3 ) on diffeomorphisms. Assume a $C^{2,1}$ open set $Y\setminus \bar{\mathcal{P}} \subset Y$ being compactly contained in $Y$ and $u\in H^2_{\#,0}(\mathcal{P})$ such that $F(u,\text{id}_{\bar{Y}})= (0,0)$ , i.e. $u$ is a solution to problem ( 3.2 ). Then there exists a neighbourhood $V\subset C^{2,1}(\bar{Y})$ of $\text{id}_{\bar{Y}}$ and a continuous function $g\,:\,V\to H^2_{\#,0}(\mathcal{P})$ , $g(\text{id}_{\bar{Y}})=u$ , such that

\begin{align*} F(g(h),h)= (0,0), \quad \forall h\in V. \end{align*}

Particularly, $h^{*-1}g(h)$ solves ( 3.2 ) up to an additive constant on $h(\mathcal{P})$ for all $h\in{\text{Diff}}^{2,1}_\Box (\bar{Y})\cap V$ .

Proof. This is an immediate consequence of the implicit function theorem for Banach spaces as given in [41] and the arguments above.

By the previous theorem, we established a continuous relation between the diffeomorphism $h$ describing the alteration of the domain $\mathcal{P}$ and the pullback of the solution $g(h)$ to the associated problem (2.3). In a final step, we show that the continuous dependence carries over to the desired quantity $\mathbb{D}$ defined in (2.3):

Corollary 3.4. Continuous dependence of $\mathbb{D}$ on diffeomorphisms. Under the assumptions of Theorem 1, the mapping

\begin{align*} R\,:\, V \to \mathbb{R},\quad h\mapsto \int _{h(\mathcal{P})} \nabla h^{*-1}g(h) \; dy \end{align*}

is continuous. Therefore, the diffusion tensor depends continuously on ${\text{Diff}}^{2,1}_\Box (\bar{Y})$ -variations of the domain $\mathcal{P}$ .

Proof. By the change of variables theorem and chain rule, we have

(3.6) \begin{align} \int _{h(\mathcal{P})} \nabla h^{*-1}g(h) \; dy &= \int _{\mathcal{P}} \nabla g_h(x) \cdot \nabla h^{-1}(h(x)) \cdot \mid \text{det} (\nabla h(x))\mid \; dx \\ &=\int _{\mathcal{P}} \nabla g_h(x) \cdot \nabla h(x)^{-1} \cdot \mid \text{det} (\nabla h(x))\mid \; dx. \nonumber \end{align}

As $g$ is continuous with respect to the $H^2$ -norm in the image space, continuity of the functional is proven.

By the previous corollary, we established the continuous behaviour of $\mathbb{D}$ on $h\in{\text{Diff}}_\Box ^{2,1}(\bar{Y})$ in the topology of $C^{2,1}(\bar{Y})$ . However, this degree of regularity is insufficient for our later purposes, cf. Theorem A.2. In order to obtain stronger results, we specify the setting more tailored to our later application. Let us now consider a smooth mapping $h\,:\,({-}S,S) \to{\text{Diff}}_\Box ^{2,1}(\bar{Y})$ for some artificial time horizon $S\gt 0$ and $h(0)=\text{id}_{\bar{Y}}$ . This relates to a 1-parametric deformation of the initial geometry.

Example 3.5. Continuous path of diffeomorphisms. Consider inclusions $Y\setminus \bar{\mathcal{P}}$ of circular shape and of different radii. In this case, a smooth diffeomorphism on $Y$ can be easily constructed by radially compressing/expanding annuli within a compact subset of $Y$ . Given two radii $0\lt r_1\leq r_2\lt \frac{1}{2}$ , a deformation mapping a circle of radius $r_2$ to a circle of radius $r_1$ is defined by

(3.7) \begin{align} h_{r_1}(y)=\begin{cases} y, & |y|\gt \frac{1}{2}, \\[4pt] (1-\xi ( |y|))y + \xi (|y|)\left(\frac{r_1}{r_2} y\right), & r_2 \leq |y|\leq \frac{1}{2}, \\[4pt] \frac{r_1}{r_2} y, & |y|\leq r_2 \end{cases} \end{align}

choosing a suitable $\xi \in C^\infty \left ([r_2,\frac{1}{2}]\right )$ with $\xi ^\prime \in C_0^\infty \left ((r_2,\frac{1}{2})\right )$ , $\xi (r_2)=1$ , $\xi (\frac{1}{2}) = 0$ , see [9]. As such, the domain remains unchanged for $|y|\gt \frac{1}{2}$ and is uniformly contracted for $|y|\leq r_2$ with a smooth convex-combination layer in between. We illustrate this procedure for the choice $r_2=0.3,\; r_1=0.1$ in Figure 4. For $r_2$ fixed, we can consider the path

\begin{align*} h\in C^0\left (0,r_2;{\text{Diff}}^{2,1}_\Box (\bar{Y})\right ), \quad h\,:\, s\mapsto h_s. \end{align*}

Then $R\circ h$ with $R$ being defined in Corollary 3.4 is also a continuous mapping. Consequently, the diffusion tensor depends continuously on $s$ .

In the following, we prove differentiability of the diffusion tensors along such 1-parametric curves. More precisely, $C^m$ -mappings $h\,:\,({-}S,S) \to{\text{Diff}}_\Box ^{2,1}(\bar{Y})$ will be considered. Accordingly, we switch the above setting to the following:

tracking along a fixed path of diffeomorphisms in comparison to (3.3). In order to obtain higher differentiability of the resolution function $g$ in Theorem 1, higher regularity of $F$ needs to be established. Revisiting the calculations performed in (3.5) we immediately see a transfer of regularity to $F$ with respect to the second variable. As the mapping is linear with respect to the first variable, we again obtain $C^m$ -regularity for $F$ . More precisely, the following theorem holds extending the smooth dependence results for simple parametric families of shapes as derived in [9, 34].

Theorem 2. Smooth dependence of solutions to ( 2.3 ) on geometric parameter. Consider a $C^m$ -curve $s\mapsto h_s$ of ${\text{Diff}}_{\Box }^{2,1}(\bar{Y})$ -embeddings and $h_0=\text{id}_{\bar{Y}}$ for $m\geq 1$ . Assume a $C^{2,1}$ open set $Y\setminus \bar{\mathcal{P}} \subset Y$ being compactly contained in $Y$ and $u\in H^2_{\#,0}(\mathcal{P})$ such that $F(u,0)= (0,0)$ , i.e. $u$ is a solution to problem ( 3.2 ). Then there exists a neighbourhood $V$ of zero and a $C^m$ -function $g\,:\,V\to H^2_{\#,0}(\mathcal{P})$ , $g(0)=u$ , such that

\begin{align*} F(g(s),s)= (0,0), \quad \forall s\in V. \end{align*}

Particularly, $h_s^{*-1}g(s)$ solves ( 3.2 ) up to an additive constant on $h(\mathcal{P})$ for all $s\in V$ .

Proof. This is again an immediate consequence of the implicit function theorem for Banach spaces as given in [41].

Again, we can leverage this regularity to the diffusivity tensors (2.2).

Corollary 3.6. Smooth dependence of $\mathbb{D}$ on geometric parameter. Under the assumptions of Theorem 2, the mapping

\begin{align*} R\,:\, V \to \mathbb{R},\quad s\mapsto \int _{h_s(\mathcal{P})} \nabla h_s^{*-1}g(s) \; dy \end{align*}

is $m$ -times continuously differentiable. Therefore, the diffusion tensor ( 2.2 ) depends $C^m$ -regularly on variations of the domain $\mathcal{P}$ along a path of diffeomorphisms of specified regularity.

Proof. Since the map $s \mapsto h_s$ is $C^m$ -regular, we establish the same degree of Fréchet differentiability in the spatial derivatives $s \mapsto \nabla h_s$ as mappings $({-}S,S) \to C^1(\bar{Y},\bar{Y})$ with respect to the corresponding norms. Revisiting the calculations performed in (3.6) and using the product rule for Fréchet differentiable functions, we conclude the assertion.

Remark 3.7. Finite dimensional parametrisation. Due to the mathematical structure of the problem, Theorem 2 holds analogously for families of diffeomorphisms that are parameterisable by a finite number of parameters. As such, the real-valued order parameter $s$ can be replaced by its vector-valued analogue, allowing for more sophisticated couplings and more complex geometries. A natural field of application is posed by two-mineral-phase solids, cf. [14], where the two interacting phases obey their distinct evolution laws.

3.3. Existence for partial coupling

In this section, we present an existence result for strong local-in-time solutions to (3.1) under the assumption of no back-coupling from the macro- to the micro-scale, cf. Figure 1. That is, we assume the evolution of the underlying pore geometry to be known a-priori. As such, the treatment of the problem is accessible more easily in comparison to the fully coupled system which we will discuss in Section 3.5. In [34], a similar fully coupled problem was solved under the assumption of effective parameters being parameterisable by the porosity $\phi$ in a smooth and a-priori known way. As opposed to this restriction, we will introduce a new order parameter $s$ which corresponds to the parametrisation in the path of diffeomorphisms used to describe evolved initial geometries. Particularly, this approach allows for the opposite description of multiple geometries which admit the same porosity. Due to the identical structure of the model, we will use the methods of [34] for our following analysis.

To this point, we considered an initial geometry $\mathcal{P}$ of class $C^{2,1}$ and a parameterised path of diffeomorphisms $h\,:\,({-}S,S)\to{\text{Diff}}_\Box ^{2,1}(\bar{Y})$ of class $C^1$ with $S\in \mathbb{R}^+$ and $h_0 = \text{id}_{\bar{Y}}$ . In order to allow for spatial variations of effective parameters, we realise the prescription of the geometry evolution by specifying $s\,:\, \Omega _T \to ({-}S,S)$ . That is, we assume the state of each microscopic unit cell $Y(t,x)$ to be given as the state of a single master unit cell $Y^*$ at time-parameter $s(t,x)$ . The corresponding setup is visualised in Figure 5 illustrating the assignment of initial conditions on $\Omega$ . Accordingly, the effective parameters of (3.1) are given by $\phi (s),\; \sigma (s),\; \mathbb{D} (s)$ . Note that we perform the necessary redefinition of the functions $\phi,\; \sigma,\; \mathbb{D}$ to mappings from an open interval of $\mathbb{R}$ without change of notation. We furthermore restrict to the non-degenerative case, i.e. we assume

(3.8) \begin{align} \forall s\in ({-}S,S)\,:\, 0\lt \phi (s) \lt 1, \quad \sigma (s) \gt 0, \quad \mathbb{D}(s)\gt 0, \end{align}

where the last inequality holds in the sense of matrices (Loewner partial ordering). This assumption is naturally fulfilled for small times $t\gt 0$ by suitably prepared initial conditions. Moreover, we restrict our consideration to locally Lipschitz reaction rates $f(c)$ generalising the linear reaction rates prescribed in [34]. Finally, we state the following anisotropic Sobolev spaces:

\begin{align*} \mathcal{X}_1\,:\!=\,W^{1,2}_r(\Omega _T) = L^r(0,T;W^{2,r}(\Omega )) \cap W^{1,r}(0,T;L^r(\Omega )), \end{align*}

with $r\gt d+2$ which play a crucial role in the subsequent existence result. Following the major steps of [34], we have

Figure 5. Left: Master unit cell $Y^*$ with interior boundary $\partial ^{\text{int}}\mathcal{P}^*$ in black corresponding to $s=0$ . Each colour corresponds to the interior boundary of a deformed cell which is reachable from $\mathcal{P}^*$ along a path of diffeomorphisms parameterised by $s$ . Right: Macroscopic domain $\Omega$ coloured according to the initial underlying geometry displayed in the left image. For the two exemplary macroscopic points $x_1,x_2\in \Omega$ the associated microscopic geometry is displayed.

Theorem 3. Existence of strong solutions, partial coupling, diffusive transport. Let $\Omega \subset \mathbb{R}^d,\; d\in \{2,3\}$ , be a $C^2$ -domain, $r\gt d+2$ with initial conditions $c_0\in W^{2-\frac{2}{r},r}(\Omega )$ and Dirichlet boundary conditions $C_0 \in W^{1-\frac{1}{2r}, 2-\frac{1}{r}}_r(\partial \Omega _T)$ being compatible in the sense of $C_0(0,\cdot )=c_0$ on $\partial \Omega$ . Moreover, let the evolution of the pore-space geometry be given by an order parameter $s\in C^1(\overline{\Omega _{T_1}})$ , $s(t,x)\in ({-}S+\epsilon,S-\epsilon ),\; \epsilon \gt 0$ , $\forall (t,x)\in \Omega _T$ , and a path $h\in C^1({-}S,S;{\text{Diff}}_{\Box }^{2,1}(\bar{Y}))$ of diffeomorphisms such that $h_0 = \text{id}_{\bar{Y}}$ . Assume the initial inclusion $Y^*\setminus \bar{\mathcal{P}} \subset Y^*$ to be a compactly contained $C^{2,1}$ open set. Let ( 3.8 ) hold true and $f$ be locally Lipschitz. Then there exists a time $0\lt T\leq T_1$ such that ( 3.1 ) admits a unique solution $c\in \mathcal{X}_1$ .

Proof. First, we note that using the results of Theorem 2 the mapping $\mathbb{D}\,:\, ({-}S,S) \to \mathbb{R}^{d,d}$ is of class $C^1$ and accordingly $\mathbb{D}\circ s \in C^1(\overline{\Omega _{T_1}})$ . Similarly, we obtain $\phi \in C^1(({-}S,S))$ as a consequence of the following representation:

\begin{align*} \phi (s) =\int \limits _{h_s(\mathcal{P})} 1\; dy = \int \limits _{\mathcal{P}} \lvert \text{det} \left (\nabla h_s\right )\rvert \; dx \end{align*}

using the change of variables theorem. In order to establish regularity for $\sigma$ , we consider $h(\partial ^{\text{int}}\mathcal{P})$ as an evolving manifold. Let $\varphi \,:\,U\subset \mathbb{R}^{d-1}\to \mathbb{R}^d$ be a suitable local parameterisation of $\partial ^{\text{int}}\mathcal{P}$ . Then, $\hat{\varphi }\,:\, ({-}S,S) \times U \to \mathbb{R}^d$ defined as

\begin{align*} \hat{\varphi }(s,x) \,:\!=\, h_s(\varphi (x)) \end{align*}

is a parameterisation of $h_s(\partial ^{\text{int}}\mathcal{P})$ for each $s\in ({-}S,S)$ . Using the regularity of $\hat{\varphi }$ , we infer continuity of the mapping

\begin{align*} s \mapsto \int \limits _{h_s(\varphi (U))} 1 \; d\sigma = \int \limits _{U} \sqrt{\text{det}\left (\nabla \hat{\varphi } (\nabla \hat{\varphi } )^T \right ) } \; dx \end{align*}

which translates to the continuity of $\sigma$ by using a partition of unity subordinate to the domains of parameterisation, cf. [22]. As a result of the continuity of all effective parameters with respect to $s$ , there exists a $\delta \in (0,1)$ such that:

(3.9) \begin{align} \forall s\in ({-}S+\epsilon,S-\epsilon ) \,:\, \delta \lt \phi (s) \lt 1-\delta, \quad \sigma (s) \gt \delta, \quad \mathbb{D}(s)\gt \delta \mathbb{1}_d. \end{align}

Note that the continuity of the eigenvalues is inherited from the continuity of $\mathbb{D}$ , cf. [20]. Following the technique presented in [34], the proof is now based on Schauder’s fixed-point theorem applied to the set

\begin{align*} \mathcal{K}_1 = \left \lbrace c\in \mathcal{X}_1 \,:\, ||c||_{\mathcal{X}_1}\leq K \right \rbrace \end{align*}

for a constant $K\geq 1$ chosen appropriately later. Apparently, $\mathcal{K}_1$ is a convex, closed and bounded subset of $\mathcal{X}_1$ . In order to apply standard linear solution theory of parabolic equations, we rewrite equation (3.1) in the following fixed-point form:

(3.10) \begin{align} \partial _t c -\nabla \cdot \left ( \frac{\mathbb{D}(s)}{\phi (s)} \nabla c\right ) = -\frac{\partial _t [\phi (s)]}{\phi (s)} \tilde{c} + \frac{\sigma (s) }{\phi (s)} f(\tilde{c}) + \frac{\mathbb{D}(s)}{\phi (s)^2} \nabla [\phi (s)] \cdot \nabla \tilde{c}. \end{align}

Now consider the mapping $\mathcal{F}_1\,:\,\mathcal{K}_1\to L^{r}(\Omega _T)$ which maps a concentration $\tilde{c}\in \mathcal{K}_1$ to the right-hand side of (3.10). Using the compact embeddings

(3.11) \begin{align} \mathcal{X}_1 \hookrightarrow W^{\frac{3}{4},\frac{3}{2}}_{2r}(\Omega _T), \quad \mathcal{X}_1 \hookrightarrow C^{\frac{1}{2},1}(\overline{\Omega _T}), \end{align}

$\mathcal{F}_1$ shows to be compact, cf. [34]. Furthermore, the parabolic theory of Theorem A.3 delivers a continuous solution operator $\mathcal{F}_2\,:\, L^{r}(\Omega _T) \to W^{1,2}_r(\Omega _T)$ to (3.10). More precisely, we apply Theorem A.3 to the above-prescribed initial and boundary conditions with coefficients defined as

\begin{align*} a_{i,j}(t,x) = \frac{\mathbb{D}_{i,j}(s(t,x))}{ \phi (s(t,x))},\quad a_i(t,x) = -\sum \limits _{j=1}^d \partial _j \left ( \frac{\mathbb{D}_{i,j}(s(t,x))}{ \phi (s(t,x))}\right ), \quad a(t,x) = 0 \end{align*}

and source term $f$ according to the right-hand side of (3.10).

As such, for a given $\tilde{c}\in \mathcal{X}_1$ , we have

\begin{align*} ||c||_{\mathcal{X}_1} \leq C_p \left (\text{DC}+ \left \lVert -\frac{\partial _t [\phi (s)]}{\phi (s)} \tilde{c} + \frac{\sigma (s)}{\phi (s)} f(\tilde{c}) + \frac{\mathbb{D}(s)}{\phi (s)^2} \nabla [\phi (s)] \cdot \nabla \tilde{c}\right \rVert _{L^r(\Omega _T)} \right ), \end{align*}

abbreviating the contributions from initial and boundary data by

\begin{align*} \text{DC}\,:\!=\,||c_0||_{W^{2-\frac{2}{r},r} (\Omega )} + ||C_0||_{W^{1-\frac{1}{2r}, 2-\frac{1}{r}}_r(\partial \Omega _T)}. \end{align*}

Due to $f$ being locally Lipschitz, there exists a monotone function $\tilde{f}$ satisfying

(3.12) \begin{align} \tilde{f}\,:\, [0,\infty ) \to \mathbb{R}, \quad |f(x)| \leq \tilde{f}(|x|), \end{align}

Similar to the estimates established in [34], we obtain using Hölder’s inequality and (3.12)

\begin{align*} ||c||_{\mathcal{X}_1} &\leq C_p \left (\text{DC} + C_s T^{\frac{1}{r}} |\Omega |^{\frac{1}{r}}\left (||\tilde{c}(t)||_{L^\infty (\Omega _T)} + \tilde{f}(||\tilde{c}||_{L^\infty (\Omega )}) \right ) \right ) \\ &+ C_p C_s \left (T^{\frac{1}{2r}} ||\nabla [\phi (s)] ||_{L^\infty (0,T; L^{2r}(\Omega ))} ||\nabla \tilde{c} ||_{L^{2r}(0,T; L^{2r}(\Omega ))} \right ), \nonumber \end{align*}

with constant $C_s$ depending on the bounds of $\sigma, \phi,\partial _t \phi,\mathbb{D}$ , cf. (3.9). By the embeddings (3.11), every appearing norm of $\tilde{c}$ is bounded by a multiple of $||\tilde{c}||_{\mathcal{X}_1}$ and therefore by a multiple of $K$ . As such, we obtain the self-mapping property of $\mathcal{F}=\mathcal{F}_2\circ \mathcal{F}_1\,:\, \mathcal{K}_1 \to \mathcal{K}_1$ for sufficiently large $K$ and sufficiently small $T$ , finishing the existence proof. Uniqueness follows analogously to [34] by considering two solutions $c_1,c_2\in \mathcal{X}_1$ as well as their difference $\bar{c}=c_2-c_1$ . Subtracting both associated equations in fixed-point form (3.10), we obtain

(3.13) \begin{align} \partial _t \bar{c} -\nabla \cdot \left ( \frac{\mathbb{D}(s)}{\phi (s)} \nabla \bar{c}\right ) = -\frac{\partial _t [\phi (s)]}{\phi (s)} \bar{c} + \frac{\sigma (s) }{\phi (s)} \left (f(c_2)-f(c_1)\right ) + \frac{\mathbb{D}(s)}{\phi (s)^2} \nabla [\phi (s)] \cdot \nabla \bar{c}. \end{align}

Testing (3.13) with $\bar{c}$ and estimating the right-hand side with Hölder’s and Young’s inequality show

(3.14) \begin{align} \frac{1}{2} \partial _t ||\bar{c}||^2_{L^2(\Omega )}&\leq \left ( \left \lVert \frac{\partial _t [\phi (s)]}{\phi (s)} \right \rVert _{L^\infty (\Omega )} + L \left \lVert \frac{\sigma (s) }{\phi (s)} \right \rVert _{L^\infty (\Omega )} + C(\epsilon )\left \lVert \frac{\mathbb{D}(s)}{\phi (s)^2} \nabla [\phi (s)] \right \rVert _{L^\infty (\Omega )}\right ) ||\bar{c}||^2_{L^2(\Omega )}\nonumber \\ &- \int \limits _\Omega \frac{\mathbb{D}(s)}{\phi (s)}\nabla \bar{c} \cdot \nabla \bar{c} \; dx + \epsilon \left \lVert \frac{\mathbb{D}(s)}{\phi (s)^2} \nabla [\phi (s)] \right \rVert _{L^\infty (\Omega )} ||\nabla \bar{c}||^2_{L^2(\Omega )}, \end{align}

where $L$ denotes the Lipschitz constant of $f$ with respect to the compact interval

\begin{align*} \left [-\max \limits _{i\in \{1,2\}} \left \lbrace ||c_i||_{C^0(\overline{\Omega _T})}\right \rbrace,\; \max \limits _{i\in \{1,2\}} \left \lbrace ||c_i||_{C^0(\overline{\Omega _T})}\right \rbrace \right ]. \end{align*}

By the uniform coercivity of $\mathbb{D}$ , we can absorb the last addend of (3.14) into the diffusion term for sufficiently small $\epsilon \gt 0$ . Uniqueness now follows from Gronwall’s inequality.

3.4. Level-set equation induced diffeomorphisms

The smooth dependence results derived in Section 3.2 are based on geometry deformation by smooth paths of diffeomorphisms. In the following, we investigate under which conditions on the initial geometry and normal velocity field alterations performed by the level-set equation (2.7) induce such paths. To do so, we make use of the method of characteristics. As a result, we can replace the assumption of a prescribed path of diffeomorphisms in Theorem 3 by a prescribed normal velocity field $v_n(t,x,y)$ and let the geometry evolve according to the level-set equation which is a much more natural setup from the viewpoint of applications. At first, we must fix the class of real-valued functions whose level sets are guaranteed to be smooth submanifolds of codimension one:

Definition 3. Regular level-set function. Let $\Gamma \subset Y$ denote the boundary of an open set of class $C^{k,\alpha }$ , $k\geq 2$ , compactly contained in $Y$ . Then we call a function $\Phi \in C^{k,\alpha }(Y)$ regular level-set function associated with $\Gamma$ iff $\Gamma = \{y\in Y\,:\, \Phi (y) = 0\}$ and $\nabla \Phi = -\nu$ on $\Gamma$ .

Remark 3.8. Equivalence of representation. For every manifold $\Gamma$ as given in Definition 3, there exists an associated regular level-set function, cf. [16] Chapter 1 or Theorem A.1 .

Remark 3.9. Level-set solution concepts. We note that viscosity solutions pose a powerful and popular weak solution concept for a wide range of non-linear first-order PDEs like the level-set equation, for a basic introduction to this topic, see e.g. [4]. However, the solutions obtained by the related theory are only continuous by definition providing insufficient regularity of the level sets, cf. Section 3.2. We therefore directly consider classical solutions in the following. Note that due to the hyperbolic nature of the non-linear level-set equation, the evolution of singularities such as shocks is to be expected within finite time. As such, this framework is particularly unsuited for obtaining solutions existing globally in time.

In order to apply the theory for non-linear first-order PDEs, we rewrite the level-set equation (2.7) in the form $F(\vec{x},\Phi,D\Phi )=0$ with $F\,:\,(\vec{x},z,\vec{p})\mapsto \mathbb{R}$ . Note that in this notation, $\vec{x}=(x,t)$ , $\vec{p}=(p,p_{d+1})$ corresponds to the spatial and temporal variable, and $D=(\nabla _x,\partial _t)$ denotes the related differential operator. Apparently, we have

\begin{align*} F(\vec{x},z,\vec{p}) = p_{d+1}+v_n(\vec{x}) |p|. \end{align*}

with the normal interface velocity $v_n$ of (2.7). Prescribing $v_n$ smooth such that it vanishes in a neighbourhood of $\{|p|=0 \}$ , we have $F\in C^2(\mathbb{R}^{2d+3},\mathbb{R})$ . Switching to the characteristic system of ODEs, the equations read

(3.15) \begin{align} \dot{\vec{p}}(s) &=- \partial _{\vec{x}} F (\vec{x}(s),z(s),\vec{p}(s)) - \partial _{z} F (\vec{x}(s),z(s),\vec{p}(s)) \vec{p}(s) = -Dv_n(\vec{x}) |p|, \nonumber \\ \dot{z}(s) &= \partial _{\vec{p}} F(\vec{x}(s),z(s),\vec{p}(s))\cdot \vec{p}(s) = v_n(\vec{x}) \frac{p}{|p|}\cdot p + p_{d+1} = 0,\\ \dot{\vec{x}}(s) &= \partial _{\vec{p}} F (\vec{x}(s),z(s),\vec{p}(s)) = \left (v_n(\vec{x}) \frac{p}{|p|}, 1\right ), \nonumber \end{align}

with initial conditions

(3.16) \begin{align} \vec{p}(0) = \left (\nabla \Phi (0,y), -v_n(y)|\nabla \Phi (0,y)| \right ), \quad z(0) = \Phi (0,y), \quad \vec{x}(0)=(y,0). \end{align}

According to (3.15), the projected characteristics of a solution to (2.7) move in normal direction to the interface with speed $v_n(\vec{x})$ . Furthermore, the function value of $\Phi$ remains constant along trajectories. As such, the zero-level-set describing the position of the fluid–solid interfaces is transported along $x$ . Note that every admissible point of the characteristic ODE system is non-characteristic, i.e. the characteristics admit a strictly monotone distance to the set of prescribed initial data $\bar{Y}\times \{0\}$ . As such, equation (2.7) admits a unique local-in-time $C^2$ -solution by standard theory [10]. Furthermore, the parameterised trajectories associated with $x$ induce a smooth path of diffeomorphisms and the artificial parameter $s$ coincides with the actual physical time $t$ . More precisely, we have the following statement:

Lemma 3.10. Smooth path of diffeomorphisms generated by level-set solution. Let a normal velocity field $v_n\in C^{4,1}(\bar{Y},\mathbb{R})$ be given. Furthermore, let a regular level-set function $\Phi _0\in C^4(\bar{Y},\mathbb{R})$ be given, such that $\bar{\mathcal{P}}=\left \lbrace \Phi _0\leq 0 \right \rbrace$ with $Y\setminus \mathcal{P}$ being compact in $Y$ . Then the local-in-time solution to the level-set equation ( 2.7 ) with respect to initial conditions $\Phi _0$ induces a $C^1$ -path $h$ in ${\text{Diff}}_\Box ^{2,1}(\bar{Y})$ such that $h_t(\bar{\mathcal{P}}) = \left \lbrace \Phi (t,\cdot )\leq 0 \right \rbrace$ for all times $t$ sufficiently small.

Proof. In the following, we investigate the regularity of the trajectories associated with $\vec{x}$ solving the systems of ODEs (3.15). More precisely, we consider the smaller system in the spatial variables $x$ and $p$ which is closed due to the time independence of $v_n$ . As a suitable Banach space for $(x,p)$ , we introduce

\begin{align*} \mathcal{X} = (C^3(\bar{Y}, \mathbb{R}))^{d} \times (C^3(\bar{Y}, \mathbb{R}))^{d}. \end{align*}

Then the reduced initial conditions (3.16) are element of $\mathcal{X}$ . It is straightforward to check that the structure function $f\,:\,\mathcal{X} \supset V \to \mathcal{X}$ of the reduced ODE system

(3.17) \begin{align} f(x,p) = \left (v_n(x) \frac{p}{|p|},\;-\nabla v_n(x) |p| \right ) \end{align}

is well-defined and locally Lipschitz continuous with respect to the norm $|.|_{\mathcal{X}}$ for a small neighbourhood $V$ of $(x(0),p(0))$ , as long as $v_n = 0$ in a neighbourhood of $\{y\in Y \,:\, |\nabla \Phi _0(y)| = 0\}$ . Due to the tubular neighbourhood theorem (Theorem A.1), this can be achieved by modifying $v_n$ appropriately without changing the solution locally at $\partial ^{\text{int}} \mathcal{P}$ . More precisely, we force $v_n$ to zero away from a neighbourhood of $\partial ^{\text{int}} \mathcal{P}$ and within a neighbourhood of $\partial Y$ in a smooth manner. By Picard–Lindelöf theorem, there exists $T\gt 0$ such that the ODE system admits a unique solution $(x,p)\in C^1(0,T; \mathcal{X})$ . Due to the convexity of the domain, we obtain $x\in C^1(0,T; C^{2,1}(\bar{Y},\mathbb{R}^d))$ . By the choice of $v_n$ , we have $\text{im}(x(s)) = \bar{Y}$ for all $s\in (0,T)$ . Since the initial conditions fulfil $x(0)=\text{id}_{\bar{Y}}\in{\text{Diff}}^{2,1}(\bar{Y})$ , the diffeomorphism property for $x(s)$ is obtained in a possibly reduced time interval $(0,T)$ . Consequently, $x\in C^1(0,T;{\text{Diff}}_\Box ^{2,1}(\bar{Y}))$ . Since the value of $\Phi$ is constant along the trajectories associated with $x$ , we finally see $h_t(\bar{\mathcal{P}}) = \left \lbrace \Phi (t,\cdot )\leq 0 \right \rbrace$ .

In order to underline the power of the above Lemma 3.10, let us reconsider the case of contracting and expanding circles as in Example 3.5. Previously, a suitable family of diffeomorphisms had to be constructed explicitly in (3.7) to apply Theorem 2. With the help of Lemma 3.10, it is now sufficient to check that the function

\begin{align*} \Phi _r(y) = r - ||y||_2 \end{align*}

has a zero-level-set of a circle of radius $r\lt 0.5$ and can be smoothed locally around $0$ to fulfil $\Phi _{r,0} \in C^4(\bar{Y})$ . The statement follows by applying Lemma 3.10 to $\Phi _{r,0}$ and $v_n \equiv 1$ . However, the presented framework also applies to much more general geometry deformations as illustrated in Figure 3 where an explicit construction of diffeomorphisms is unsuitable. It especially includes the transition to non-convex shapes as well as the treatment of disconnected components of solid.

3.5. Existence for full coupling

In geoscientific applications of our model, the evolution of the microscopic geometry of the porous medium is not known in advance but depends on the time-dependent macroscopic concentration field. Accordingly, this chapter is dedicated to the local-in-time existence of strong solutions to a model where the geometry evolution is described by the level-set equation (2.7). For clarity, we restate the model equations:

\begin{align*} \phi \partial _t c - \nabla \cdot (\mathbb{D}\nabla c) &= \sigma f(c) - \partial _t \phi c &&\quad \text{in } (0,T)\times \Omega, \\ \frac{\partial \Phi }{\partial t} + v_n |\nabla _y \Phi | &= 0, &&\quad \text{in } (0, T)\times \Omega \times Y, \nonumber \end{align*}

supplemented by boundary and initial conditions as well as diffusion cell problem (2.3), where the level-set normal interface velocity is given by

(3.18) \begin{align} v_n(t,x) = f(c(t,x)), \end{align}

i.e. depending on the solution of the macroscopic concentration. In simplification of (2.9), we consider a uniform velocity within each unit cell to facilitate the establishment of higher regularity for quantities derived from the geometry. Furthermore, we restrict to scenarios where the porosity can be used as the natural order parameter uniquely characterising the geometry state which is the typical case in dissolution/precipitation processes. Moreover, we assume linear reaction rates. By doing so, we retract to the setting investigated in [34] enabling us to use respective results stated therein. Additionally, we again consider a single evolving master unit cell $Y^*$ as in Section 3.3 and introduce spatial inhomogeneity in $\Omega$ by choosing different evolution states of $Y^*$ as the initial condition.

In the following, we approach existence of solutions $(c,\Phi )\in (\mathcal{X}_1, \Omega _x\times C^1((0,T)\times Y))$ to the fully coupled model described above, where the solution space with respect to $\Phi$ is defined as

\begin{align*} \Omega _x\times C^1((0,T)\times Y) = \left \lbrace \Phi \,:\, (0, T)\times \Omega \times Y \to \mathbb{R}\,:\, \; \forall x\in \Omega \quad \Phi (\cdot,x,\cdot ) \in C^1((0,T)\times Y) \right \rbrace. \end{align*}

To do so, we rewrite the given problem in the variables $(c,\phi )\in \mathcal{X}_1^2$ . As such, we can apply Theorem A.2 to obtain local-in-time existence in the variables $(c,\phi )$ . In order to do so, the effective parameters $\mathbb{D},\sigma$ must be expressed as functions of $\phi$ . As we will see, this relation is independent of the concentration solution $c$ and is therefore determined by the level-set initial condition alone. Given the solution $(c,\phi )$ of the transformed system, we construct the solution $(c,\Phi )$ to our model of consideration.

Let us consider the geometry evolution being preliminarily parameterised by the time $t$ of the level-set equation. Using the properties of (3.18), we can reparametrise all effective parameters as a function of $\phi$ by setting:

(3.19) \begin{align} \hat{\sigma } ={\sigma } \circ \phi ^{-1}, \quad \hat{\mathbb{D}} ={\mathbb{D}} \circ \phi ^{-1}. \end{align}

The main difficulty in proving existence for the bilaterally coupled system is establishing sufficiently high regularity in the coefficients $\hat{\sigma }, \hat{\mathbb{D}}$ as required to apply Theorem A.2.

First, we notice that the mapping $\phi \mapsto \hat{\sigma }, \hat{\mathbb{D}}$ is independent of the particular choice of a normal velocity $v_n$ in (3.18) in the set of y-independent functions which can be seen as follows: Let $\Phi _1$ solve (2.7) with $v_n \equiv 1$ . Then $\Phi (t,x) = \Phi _1(V(t),x)$ solves (2.7) with $v_n =v$ and $V$ being the primitive of $v$ . As such, solutions are solely rescaled with respect to the time variable but not with respect to space, i.e. the individual geometries relating to $\Phi _1(s,\cdot )$ remain unaffected. It is therefore sufficient to consider the case $v_n\equiv 1$ , where, due to

(3.20) \begin{align} \partial _s \phi (s) = -\sigma (s), \end{align}

$\phi$ is strictly monotonically increasing. This fact justifies the transformations introduced in (3.19). As such, the regularity of $\phi \mapsto \hat{\sigma }, \hat{\mathbb{D}}$ follows immediately from the regularity of $\phi, \sigma, \mathbb{D}$ which we investigate using the same means as in Theorem 3.

Let $\varphi \,:\,U\subset \mathbb{R}^{d-1}\to \mathbb{R}^d$ be a suitable local parameterisation of $\partial ^{\text{int}}\mathcal{P}$ . Then, $\hat{\varphi }\,:\, ({-}S,S)\times U \to \mathbb{R}^d$ defined as

(3.21) \begin{align} \hat{\varphi }(s,x) \,:\!=\, \varphi (x) + s\nu (\varphi (x)) \end{align}

is a local parameterisation of the deformed interface at times $s\in ({-}S,S)$ using the results of Section 3.4. Note that $\hat{\varphi }(s,\cdot )$ is a bijection for $s$ close to zero and in case of $\Phi _0 \in C^4(Y)$ , the map (3.21) is of class $C^3$ , cf. Theorem A.1. Using the regularity of $\hat{\varphi }$ , we infer the mapping

\begin{align*} s \mapsto \int \limits _{U} \sqrt{\text{det}\left (\nabla \hat{\varphi } (\nabla \hat{\varphi } )^T \right ) } \; dx \end{align*}

to be of class $C^2$ , which translates to $\sigma (s)$ by using a partition of unity, cf. [22]. By (3.20), we also obtain $\phi \in C^3(({-}S,S))$ . As such, we finally conclude $\hat{\mathbb{D}} \in C^1((\phi _{\text{min}},\phi _{\text{max}}))$ , $\hat{\sigma } \in C^2((\phi _{\text{min}},\phi _{\text{max}}))$ for sufficiently narrow bounds $\phi _{\text{min}}\lt \phi ^0\lt \phi _{\text{max}}$ around the initial condition. Summing up the observations of this section, we obtain the following existence result using Theorem A.2.

Theorem 4. Existence of strong solutions, full coupling, diffusive transport. Let $\Omega \subset \mathbb{R}^d,\; d\in \{2,3\}$ , be a $C^2$ -domain, $r\gt d+2$ with initial conditions $c_0\in W^{2-\frac{2}{r},r}(\Omega )$ , $c_0\geq 0$ , fulfilling $c_0=0$ on $\partial \Omega$ . Furthermore, consider an initial inclusion compactly contained in $Y^*$ given by a regular $C^4(\bar{Y}^*)$ level-set function evolving according to the level-set equation ( 2.7 ) such that solutions for $v_n \equiv \pm 1$ cover $\phi \in (\phi _{\text{min}},\phi _{\text{max}}) \subset \subset (0,1)$ . Let $\phi ^0\in W^{2,r}(\Omega )$ , $\phi ^0(x)\in (\phi _{\text{min}}+\epsilon,\phi _{\text{max}}-\epsilon )$ for some $\epsilon \gt 0$ , $x\in \Omega$ and the reaction rate given as $f(c)=c$ . Then there exists a local-in-time solution $(c,\Phi )\in (\mathcal{X}_1, \Omega _x\times C^1((0,T)\times Y))$ to the bilaterally coupled model. In case of higher regularity, i.e. $c_0\in W_r^{4-\frac{2}{r}}(\Omega )$ , $\phi ^0\in W_r^{3-\frac{2}{r}}(\Omega )$ fulfilling the compatibility condition $(\nabla c_0 \cdot \nabla \phi ^0) \equiv 0$ , $\Delta c_{0|\partial \Omega }\equiv 0$ it either holds for the maximum existence time $T$ :

(3.22) \begin{align} \lim \limits _{t\to T} \inf \limits _{x\in \Omega } \phi (t,x) = \phi _{\text{min}} \quad \quad \text{or} \quad \quad T=\infty. \end{align}

Proof. In the setup presented above, Theorem A.2 ensures the short-time existence of solutions $(c,\phi )\in (\mathcal{X}_1)^2$ to the equations

\begin{align*} \phi \partial _t c -\nabla \cdot (\hat{\mathbb{D}}(\phi )\nabla c) &= \hat{\sigma }(\phi ) c^2 - \hat{\sigma }(\phi ) c,\\ \nonumber \partial _t\phi = -\hat{\sigma }(\phi ) c, \end{align*}

with $\hat{\mathbb{D}},\; \hat{\sigma }$ as constructed in (3.19). Apparently, the second equation corresponds to the level-set function $\Phi$ fulfilling

(3.23) \begin{align} \partial _t \Phi (t,x,y) +c(t,x)|\nabla _y \Phi (t,x,y)| = 0, \quad \forall x\in \Omega, \end{align}

for $y$ sufficiently close to the zero-level-set as required. By the embedding theorems (3.11), $t\mapsto c(t,x)$ is continuous for all $x\in \Omega$ . Furthermore, a solution $\Phi _1$ related to $v_n \equiv 1$ in a neighbourhood of the zero-level-set exists and is of class $C^2$ as discussed in Section 3.4. As such, the transformed solution $\Phi (t,x,y) = \Phi _1(C(t,x),x,y)$ with $C(\cdot,x)$ being a primitive of $c(\cdot,x)$ solving (3.23) is of class $C^1((0,T)\times Y)$ for each $x\in \Omega$ . Finally, we investigate the case of high regularity initial conditions. Note that the structure function f in (3.17) is smooth by the regularity of $v_n$ , leading to a solution $x$ smooth in time. Therefore, we particularly have $\hat{\mathbb{D}} \in C^3((\phi _{\text{min}},\phi _{\text{max}}))$ due to Theorem 2. As such, the assumption follows by Corollary 1 of [34].

Remark 3.11. Global existence. The range of porosities $(\phi _{\text{min}},\phi _{\text{max}})$ the solutions in the above theorem attain depends on the radius $r$ of the tubular neighbourhood of the master unit cell’s fluid–solid interface at $s=0$ , cf. Theorem A.1. More precisely, $\phi _{\text{max}}-\phi _{\text{min}}$ is bounded by the volume of the tubular neighbourhood, beyond which the appearance of shocks prevents the continuation of solutions. Consequently, global-in-time existence cannot be expected in general due to the requirement of classical solutions to the non-linear hyperbolic level-set equation in our framework, cf. Remark 3.9. In the case of circular geometries as illustrated in Example 3.5, we can allow for $(\phi _{\text{min}},\phi _{\text{max}})\subset \subset (1-\frac{\pi }{4},1)$ , i.e. the full range between complete dissolution and clogging. Moreover, under the stronger assumptions on the smoothness of parameters and the additional compatibility conditions of the initial data the existence interval of $(c,\Phi )$ is solely limited by the well-behavedness of the underlying microscopic geometry evolution.

Remark 3.12. Uniqueness. Note that a proper level-set function uniquely determines the geometry state of the system but not vice versa. In that sense, the problem of Theorem 4 features multiple solutions $(c,\Phi )$ . However, all these solutions are equivalent as they describe the same geometry evolution. By Theorem 4.2 in [34], the solution $(c,\phi )\in \mathcal{X}_1^2$ is unique. By the one-on-one relation between geometry state and porosity, this translates to the uniqueness of the solid part within each unit cell $Y(t,x)$ . Fixing an initial value $\Phi _0$ , also the solution $(c,\Phi )$ is uniquely determined.

4. Smooth parameter dependence and existence for diffusive–advective transport

In this section, we consider an extension of the model analysed in Section 3 involving fluid flow and including advective solute transport. To do so, we follow a similar strategy as in Section 3. After fixing the setting in Section 4.1, we prove smooth dependence of the permeability tensor on the geometry using a variant of the implicit function theorem in Section 4.2. Building upon those results, the continuous dependence of Darcy velocity and pressure on the permeability field is shown. Finally, we obtain local-in-time existence results for the diffusive–advective transport case with partial one-way micro-to-macro coupling in Section 4.4. A graphical representation of the roadmap for this section is found in Figure 6.

Figure 6. Roadmap for Section 4. Main statements are given in blue boxes, linking lemmas, theorems and corollaries are highlighted in red boxes.

4.1. Setting

In the following, we consider solute transport by diffusion and advection:

(4.1) \begin{align} \phi \partial _t c +\nabla \cdot (vc)- \nabla \cdot (\mathbb{D}\nabla c) = \sigma f(c) - \partial _t \phi c\quad \text{in } (0,T)\times \Omega \end{align}

which is coupled to Darcy’s equation

(4.2) \begin{align} v &= -\mathbb{K} \nabla p &&\text{in } \Omega,\; t\in (0,T),\\ \nabla \cdot v &=\tilde{f} && \text{in } \Omega,\; t\in (0,T). \nonumber \end{align}

In comparison to the equations stated in Section 2, we allow for a general constant-in-time divergence $\tilde{f}$ . In contrast to the previous analysis, we now consider a partial coupling from the microscopic to the macroscopic scale as in Section 3.3. That is, we suppose the evolution of the underlying geometry is a-priori given by a sufficiently smooth path of diffeomorphisms. Finally, the model regarded in this chapter is closed by the two effective tensors and associated cell problems (2.3), (2.6).

4.2. Continuous dependence of permeability tensors

Using a similar strategy as in Section 3.2, we also show the continuous dependence of the permeability tensor on the geometry evolution.

In order to capture the underlying stationary Stokes equations in (2.6), we introduce the following function spaces for the velocity field $v$ and pressure field $q$ for $k\geq 1$ , cf. [12]:

\begin{align*} H^k_v(\mathcal{P}) &= \left \lbrace v \in (H_\#^k(\mathcal{P}))^d\,:\,{\text{Tr}}_{\mathcal{P}}(v) = 0 \right \rbrace,\\[3pt] H^{k-1}_q(\mathcal{P}) &= \left \lbrace q \in H_\#^{k-1}(\mathcal{P})\,:\,\int _{\mathcal{P}} q \; dy = 0 \right \rbrace, \\[3pt] H^k_{v,\Box }(\mathcal{P}) &= \left \lbrace v\in H^k_v(\mathcal{P})\,:\, \nabla \cdot v = 0 \text{ on } U \right \rbrace \end{align*}

for some fixed neighbourhood $\partial Y \subset U \subset \bar{Y}$ , cf. (3.4). In order to formulate the notion of weak solutions stated in [12], we furthermore introduce the space

\begin{align*} H^1_{v,\sigma }(\mathcal{P})\,:\!=\,\left \lbrace u\in H^1_{v}(\mathcal{P})\,:\,\nabla \cdot u =0 \right \rbrace \end{align*}

of solenoidal functions in $H^1_{v}(\mathcal{P})$ . As in Section 3.2, we are required to consider the general inhomogeneous class of Stokes equations. Therefore, the weak form of the Stokes problem with general force term $f\in (L^2(\mathcal{P}))^d$ and inhomogeneity $j\in H^1(\mathcal{P})\cap L^2_0(\mathcal{P})$ with $j=0$ in a neighbourhood of $\partial Y$ reads: Find $(u,p)\in H^1_v(\mathcal{P}) \times H^0_q(\mathcal{P})$ such that

(4.3) \begin{align} \int _{\mathcal{P}} \nabla u \,:\, \nabla v \; dy &= \int _{\mathcal{P}} f v \; dy, \nonumber \\ \int _{\mathcal{P}} \nabla u \,:\, \nabla \Psi \; dy - \int _{\mathcal{P}} p \nabla \cdot \Psi \; dy &= \int _{\mathcal{P}} f \Psi \; dy,\\ \nabla \cdot u &= j, \nonumber \end{align}

for all $(v,\Psi )\in (H^1_{v,\sigma }(\mathcal{P}) \times C^\infty _0(\mathcal{P}))$ .

Following the procedure of Section 3.2, we next implement higher regularity of solutions to the Stokes problem given a sufficiently smooth interior boundary $\partial ^{\text{int}}\mathcal{P}$ .

Lemma 4.1. Regularity Stokes. Let $\partial ^{\text{int}}\mathcal{P}$ be of class $C^2$ , $f\in (L^2(\mathcal{P}))^d$ and $j\in H^1(\mathcal{P})\cap L^2_0(\mathcal{P})$ vanishing in a neighbourhood of $\partial Y$ . Then there exists a unique solution $(u,p)$ to ( 4.3 ) in $H^2_v(\mathcal{P}) \times H^1_q(\mathcal{P})$ . Furthermore, the associated solution operator is linear and bounded.

Proof. We start considering the homogeneous case $j\equiv 0$ . By standard Hilbert space arguments, there exists a uniquely determined function $u\in H^1_{v,\sigma }(\mathcal{P})$ such that

\begin{align*} \int _{\mathcal{P}} \nabla u \,:\, \nabla v \; dy = \int _{\mathcal{P}} f v \; dy, \quad \forall v \in H^1_{v,\sigma }(\mathcal{P}). \end{align*}

Due to the smoothness of the domain and source term, there exists a unique pressure field $p\in L^2(\mathcal{P})$ with vanishing average such that

\begin{align*} \int _{\mathcal{P}} \nabla u \,:\, \nabla \Psi \; dy - \int _{\mathcal{P}} p \nabla \cdot \Psi \; dy &= \int _{\mathcal{P}} f \Psi \; dy, \quad \forall \Psi \in (C_0^\infty (\mathcal{P}))^2 \end{align*}

by Lemma IV 1.1 in [12]. Furthermore, Theorem IV 4.1 in [12] ensures the interior regularity claimed. As the boundary is of class $C^2$ , the regularity can be extended to the full domain by Theorem IV 5.1 in [12]. Next, we consider the problem for general inhomogeneities $j$ of the specified class. According to Theorem 3.4 in [15], there exists a function $\beta \in H^2_0(\mathcal{P})$ satisfying $\nabla \cdot \beta = j$ . Solving the homogeneous system for $\tilde{f}=f-\Delta \beta$ and adding $\beta$ to the velocity solution, the full statement is shown. Combining the estimates associated with the previous steps, we conclude boundedness of the solution operator.

In comparison to the case of elliptic equations in Section 3.2, we need to slightly change the setting here in order to accommodate for the additional condition on $\nabla \cdot u$ which is not invariant under pullbacks and cannot be additively compensated for without changing the solution itself. More precisely, we switch to a setting that does not require surjectivity of operators. Serving the analogous purpose as $F$ defined in equation (3.3), let us consider the following mappings:

\begin{align*} G & \,:\, H^2_{v,\Box }(\mathcal{P}) \times H^1_{q}(\mathcal{P}) \times ({-}S,S) \to L^2(\mathcal{P}) \times H^1_{\#}(\mathcal{P}),\\ & (u,p,s) \mapsto \left (h_s^*\Delta h_s^{*-1} u- h_s^*\nabla h_s^{*-1} p,\; h_s^*\nabla \cdot h_s^{*-1} u \right ),\\ \tilde{G} & \,:\, ({-}S,S) \to L^2(\mathcal{P}) \times H^1_{\#}(\mathcal{P}), \\ & (s) \mapsto (h_s^*e_1,\; 0), \end{align*}

with $h$ being a $C^1$ -path of diffeomorphisms in ${\text{Diff}}^2_{\Box }(\bar{Y})$ and $h_0=\text{id}_{\bar{Y}}$ . As the interior boundary conditions are invariant under pullbacks and already implemented in the underlying function spaces, it is sufficient to only map to the function spaces associated with the bulk data $j,\; f$ . Differentiability of this mapping is established using the same reasoning as above. Summing up our considerations, we obtain the following theorem.

Theorem 5. Smooth dependence of weak solutions to ( 2.6 ) on geometric parameter. Assume a $C^{2}$ open set $Y\setminus \bar{\mathcal{P}} \subset Y$ being compactly contained in $Y$ and $(u,p)\in H^2_{v,\Box }(\mathcal{P}) \times H^1_q(\mathcal{P})$ such that $G(u,p,0)= \tilde{G}(0)$ , i.e. $(u,p)$ is a solution to problem ( 4.3 ) with force term $f=e_1$ and heterogeneity $j \equiv 0$ . Furthermore, let $h$ be an $m$ -times differentiable path in ${\text{Diff}}_{\Box }^2(\bar{Y})$ with $h_0=\text{id}_{\bar{Y}}$ . Then there exist a neighbourhood $V$ of zero and a function $g\,:\,V\to H^2_v(\mathcal{P}) \times H^1_{q}(\mathcal{P})$ $m$ -times differentiable in $0\in V$ , $g(0)=(u,p)$ , such that

\begin{align*} G(g(s),s)= \tilde{G}(s), \quad \forall s\in V, \end{align*}

i.e. $h_s^{*-1}g(s)$ solves ( 4.3 ) on $h(\mathcal{P})$ for $s\in V$ .

Proof. By calculations similar to Section 3.2, we find $G,\; \tilde{G}$ to be differentiable in $s=0$ . Furthermore, $G(\cdot,\cdot,s)$ is a bounded linear operator for fixed $s$ . The unique solvability of

\begin{align*} G(u,p,s)= \tilde{G}(s) \end{align*}

is guaranteed for all $s\in V$ by Lemma 4.1, i.e. we can properly define the map $g(s)$ . Finally, we obtain the estimate

\begin{align*} ||G(u,p;0)||_{L^2(\mathcal{P}) \times H^1(\mathcal{P})} \geq C^{-1} ||(u,p)||_{H^2(\mathcal{P}) \times H^1(\mathcal{P})} \quad \forall (u,p) \in H^2_{v,\Box }(\mathcal{P}) \times H^1_q(\mathcal{P}) \end{align*}

using the operator norm $C\lt \infty$ of the Stokes solution operator on $\mathcal{P}$ introduced in Lemma 4.1. The assertion is a consequence of Theorem A.4. This is due to the fact that all values in the second image space of $G$ vanish close to $\partial Y$ by the restrictions on the preimage spaces.

This result immediately translates to the permeability tensor:

Corollary 4.2. Smooth dependence of $\mathbb{K}$ on geometric parameter. Under the assumptions of Theorem 5, the mapping

\begin{align*} R\,:\, V \to \mathbb{R},\quad s\mapsto \int _{h_s(\mathcal{P})} h_s^{*-1}g_u(s) \; dy \end{align*}

is $m$ -times differentiable. Therefore, the permeability tensor, cf. ( 2.5 ), depends differentiably on variations of $s$ .

4.3. Continuous dependence of Darcy velocity and pressure

The Darcy velocity field $v$ enters the first-order term of the transport equation (2.1) as a parameter. Again, in order to apply linear parabolic theory as in Section 3.3, we first need to establish sufficient regularity of $v$ which immediately rises the question of dependence on the permeability $\mathbb{K}$ . A special difficulty arises from the fact that the stationary Darcy equation (2.4) acts on time-slices of the space-time cylinder $\Omega _T$ . In this chapter, we investigate the effect of permeability tensors continuously depending on space and time on the flow field fulfilling Darcy’s equation. The following result establishes local Lipschitz continuity with respect to the $L^2$ -norm in the velocity field. Assume a typical flow-channel scenario with flux boundary conditions on $\partial \Omega _{\text{flux}}$ and an outlet $\partial \Omega _{\text{Dir}}$ with zero Dirichlet data for the pressure with $\partial \Omega = \partial \Omega _{\text{Dir}} \cup \partial \Omega _{\text{flux}}$ and Hausdorff measure $\mathcal{H}^{d-1} (\partial \Omega _{\text{Dir}}) \gt 0$ . First, we consider (2.4) as an elliptic equation for the pressure $p$ in the following weak form with a general source term $\tilde{f}\in L^2(\Omega )$ :

Find $p\in H_{\text{Dir}}^1(\Omega )\,:\!=\,\left \lbrace v\in H^1(\Omega )\,:\, v=0 \text{ on } \Omega _{\text{Dir}}\right \rbrace$ such that for all $q \in H_{\text{Dir}}^1(\Omega )$

(4.4) \begin{align} \int _{\Omega } \mathbb{K} \nabla p \cdot \nabla q \; dx - (g_{\text{flux}},q)_{L^2(\partial \Omega _{\text{flux}}) }= \int _{\Omega } \tilde{f} q \; dx. \end{align}

The associated velocity field is then given as

\begin{align*} v=-\mathbb{K} \nabla p. \end{align*}

As such, we can establish regularity of $v$ by analysing the pressure equation (4.4) and deduce the relevant properties from $p$ . Next, we follow the approach taken in [6] where continuous dependence on parameters in case of the Brinkman–Forchheimer equation was established. Considering two pressure solutions to Darcy’s equation $(p_1,p_2)$ related to a pair of coefficients $(\mathbb{K}_1, \mathbb{K}_2)$ , we test the difference of the weak formulations with $p_2-p_1$ . Using suitable estimates on the new right-hand side, the assertion is established. More precisely, the following statement holds:

Lemma 4.3. Continuous dependence of Darcy solutions on $\mathbb{K}$

Let $\Omega$ be a Lipschitz bounded and connected domain with flux boundary conditions on $\partial \Omega _{\text{flux}}$ and homogeneous pressure boundary conditions on $\partial \Omega _{\text{Dir}}$ :

\begin{align*} \mathbb{K}\nabla p_{|\partial \Omega _{\text{flux}}} \cdot \nu = g_{\text{flux}} \in L^2(\partial \Omega _{\text{flux}}), \quad p = 0 \in H^{\frac{1}{2}} (\partial \Omega _{\text{Dir}}). \end{align*}

Assume that the Dirichlet boundary part $\partial \Omega _{\text{Dir}}$ is of positive measure. Furthermore, let $\mathbb{K}_1, \mathbb{K}_2 \in \left ( L^\infty \left ({\Omega }\right ) \right )^{d \times d}$ be permeability tensor fields uniformly $\lambda \gt 0$ -coercive and uniformly bounded in the Frobenius norm a.e. by $K_{max}$ . In addition, let $\tilde{f}\in L^2(\Omega )$ . Denoting the respective solutions to Darcy’s equation by $(v_1,p_1), \; (v_2,p_2) \in L^2(\Omega )\times H^1(\Omega )$ the following estimates hold:

\begin{align*} ||p_2-p_1||_{H^1(\Omega )} \leq C ||\mathbb{K}_2-\mathbb{K}_1||_{L^\infty (\Omega )},\\ ||v_2-v_1||_{L^2(\Omega )} \leq C ||\mathbb{K}_2-\mathbb{K}_1||_{L^\infty (\Omega )}, \nonumber \end{align*}

where the constant $C$ only depends on $\lambda$ , $K_{max}$ , $\Omega$ and the Darcy data $g_{\text{flux}},\tilde{f}$ .

Proof. First of all, we note that standard elliptic theory guarantees the existence of solutions $p_1,p_2\in H_{\text{Dir}}^1(\Omega )$ to (4.4). In a first step towards continuous dependence of solutions on $\mathbb{K}$ , we show uniform boundedness of $p$ in the $H^1$ -seminorm by revisiting well-established energy methods. As $\mathbb{K}$ is symmetric positive definite, there exists a unique symmetric root $\mathbb{K}^{\frac{1}{2}}$ . As $p$ itself is an admissible test function, the weak formulation directly yields for an arbitrary $\epsilon \gt 0$ :

\begin{align*} \lambda || \nabla p ||^2_{L^2(\Omega )} &\leq || \mathbb{K}^{\frac{1}{2}} \nabla p ||^2_{L^2(\Omega )} \leq (g_{\text{flux}},p)_{L^2(\partial \Omega _{\text{flux}}) } +||\tilde{f}||_{L^2(\Omega )} ||p||_{H^1(\Omega )} \nonumber \\ &\leq C(\epsilon ) \left (||g_{\text{flux}}||^2_{L^2(\partial \Omega _{\text{flux}})}+ ||\tilde{f}||^2_{L^2(\Omega )} \right ) +\epsilon ||\nabla p||^2_{L^2(\Omega )} \end{align*}

using Young’s and Poincaré’s inequality as well as the trace theorem. The eigenvalues of $\mathbb{K}$ being bounded from below by $\lambda$ and choosing $\epsilon$ small enough, we have

\begin{align*} ||\nabla p ||_{L^2(\Omega )}&\leq C \left ( || g_{\text{flux}}||_{L^2(\partial \Omega _{\text{flux}})}+||\tilde{f}||_{L^2(\Omega )} \right ) \end{align*}

with $C$ depending on Poincaré’s constant and $\lambda$ only.

Next, assume $p_1, \; p_2$ to be pressure solutions to the Darcy problems respective to $\mathbb{K}_1,\mathbb{K}_2$ . Testing the difference of the weak formulations with $p\,:\!=\,p_2-p_1$ leads to

\begin{align*} ||\mathbb{K}_1^{\frac{1}{2}} \nabla p ||^2_{L^2(\Omega )} &= \int _\Omega (\mathbb{K}_1 -\mathbb{K}_2) \nabla p_2 \cdot \nabla p \; dx \\ & \leq \mathop{\textrm{ess sup}}\limits _{x\in \Omega }(||\mathbb{K}_1 -\mathbb{K}_2||_F) \left ( C(\epsilon ) ||\nabla p_2||^2_{L^2(\Omega )} + \epsilon ||\nabla p||^2_{L^2(\Omega )} \right ) \nonumber \\ & \leq \mathop{\textrm{ess sup}}\limits _{x\in \Omega }(||\mathbb{K}_1 -\mathbb{K}_2||_F) C(\epsilon ) ||\nabla p_2||^2_{L^2(\Omega )} + 2K_{max} \epsilon ||\nabla p||^2_{L^2(\Omega )}, \nonumber \end{align*}

where $||.||_F$ denotes the Frobenius norm of the matrices. Using the uniform boundedness in the gradient of $p_2$ established in the first step and choice of $\epsilon$ small enough yields the assertion on the pressure. Calculating

\begin{align*} ||v_1 -v_2||_{L^2(\Omega )} &= ||\mathbb{K}_1 \nabla p_1 -\mathbb{K}_2 \nabla p_2 ||_{L^2(\Omega )} \\ &\leq ||\mathbb{K}_1||_{L^\infty (\Omega )} ||\nabla p ||_{L^2(\Omega )} + ||\mathbb{K}_1-\mathbb{K}_2||_{L^\infty (\Omega )} ||\nabla p_2 ||_{L^2(\Omega )} \nonumber \\ &\leq C ||\mathbb{K}_2-\mathbb{K}_1||_{L^\infty (\Omega )},\nonumber \end{align*}

we obtain the assertion on the velocities applying the same reasoning.

Unfortunately, the last result cannot be generalised to stronger norms in $v$ by demanding higher regularity to the data, coefficient or $\Omega$ due to the mixed boundary conditions. For a counter-example in a setting of maximal smoothness see [38].

Yet, in the case of pure Dirichlet boundaries, it is straightforward to show that the above estimate indeed can be refined to hold in $L^\infty (\Omega )$ by standard elliptic theory [10]. We capture this observation in the following corollary.

Corollary 4.4. Stronger norms. Let the assumptions of Lemma 4.3 hold with $\partial \Omega _{\text{flux}}=\emptyset$ , a $C^3$ -domain $\Omega \subset \mathbb{R}^d$ , $\tilde{f} \in H^1 (\Omega )$ and $d\in \{2,3\}$ . Furthermore, let $\mathbb{K}_1,\mathbb{K}_2 \in (C^2(\bar{\Omega }))^{d\times d}$ be uniformly coercive and uniformly bounded. Then

\begin{align*} ||v_2-v_1||_{L^\infty (\Omega )} \leq C ||\mathbb{K}_2-\mathbb{K}_1||^{1-\frac{d}{4}}_{L^\infty (\Omega )}. \end{align*}

Proof. Using standard elliptic theory, the solutions $p_1,p_2$ to (4.2) are uniformly bounded in $W^{3,2}(\Omega )$ , cf. [10]. By the Gagliardo–Nirenberg interpolation theorem [28] we have for $p=p_2-p_1$

\begin{align*} ||\nabla p||_{L^\infty (\Omega )} \leq C_1 ||D^3p||^{\frac{d}{4}}_{L^2(\Omega )} ||\nabla p||_{L^2(\Omega )}^{1-\frac{d}{4}} + C_2 ||\nabla p||_{L^2(\Omega )}. \end{align*}

Using the uniform boundedness of $||D^3p||_{L^2(\Omega )}$ and Lemma 4.3, we see

\begin{align*} ||\nabla p||_{L^\infty (\Omega )} \leq C ||\mathbb{K}_2-\mathbb{K}_1||_{L^\infty (\Omega )}^{1-\frac{d}{4}} \end{align*}

which translates to the statement claimed.

Remark 4.5. Continuity on space-time cylinder. In case $\mathbb{K}\,:\,\Omega _T\to \mathbb{R}^{d\times d}$ fulfils the conditions of Corollary 4.4 for each time-slice and the additional condition $\mathbb{K} \in C\left (0,T; (L^\infty (\Omega ))^{d\times d}\right )$ holds, Darcy’s velocity $v\,:\,\Omega _T\to \mathbb{R}^d$ is a continuous function on the whole space-time cylinder.

4.4. Existence for partial coupling

Analogously to Section 3.3, we use our prior results to prove local-in-time existence of solutions to the model specified in Section 4.1. Using the prior regularity results of this section, sufficient smoothness of the Darcy velocity field is established. Again, we re-define the functions $\phi,\; \sigma,\; \mathbb{D},\; \mathbb{K}$ from mappings on $\Omega _T$ , cf. (3.1), to mappings of the real-valued order parameter without change of notation. In order to avoid additional difficulties due to degeneracy, we require

(4.5) \begin{align} \forall s\in ({-}S,S)\,:\, 0\lt \phi (s) \lt 1, \quad \sigma (s) \gt 0, \quad \mathbb{D}(s)\gt 0, \quad \mathbb{K}(s)\gt 0, \end{align}

to hold along the prescribed path $h$ of diffeomorphisms within the master unit cell $Y^*$ .

Theorem 6. Existence of strong solutions, partial coupling, advective transport. Let $\Omega \subset \mathbb{R}^d,\; d\in \{2,3\}$ be a $C^3$ -domain. Concerning the transport equation ( 4.1 ), let $r\gt d+2$ hold with initial conditions $c_0\in W^{2-\frac{2}{r},r}_r(\Omega )$ and boundary conditions $C_0 \in W^{1-\frac{1}{2r}, 2-\frac{1}{r}}(\partial \Omega _T)$ being compatible in the sense of $C_0(0,\cdot )=c_0$ on $\partial \Omega$ . Concerning Darcy’s equation ( 4.2 ), let $p_{|\partial \Omega } = 0$ and $\tilde{f}\in H^2(\Omega )$ . Furthermore, let the evolution of the pore-space geometry be given by an order parameter $s\in C^1([0,T_1);C^2(\bar{\Omega }))$ , $s(t,x)\in ({-}S+\epsilon,S-\epsilon )$ , $ \epsilon \gt 0\; \forall (t,x)\in \Omega _T$ , and a path $h \in C^3({-}S,S;{\text{Diff}}_{\Box }^{2,1}(\bar{Y}))$ of diffeomorphisms such that $h_0 = \text{id}_{\bar{Y}}$ . Assume the initial inclusion $Y^*\setminus \bar{\mathcal{P}} \subset Y^*$ to be a compactly contained $C^{2,1}$ open set. Let ( 4.5 ) hold true and the reaction rate $f$ be locally Lipschitz. Then there exists a time $0\lt T\leq T_1$ such that the system ( 4.1 ), ( 4.2 ) admits a unique solution $(c,v,p)\in \mathcal{X}_1\times C(0,T,L^{2}(\Omega )) \times C(0,T,W^{1,2}(\Omega ))$ .

Proof. The proof of this assertion follows along the lines of the proof of Theorem 3. In addition, we must now also consider the smoothness of $\mathbb{K}(s(t,x))$ as well as of the associated Darcy velocity field $v(t,x)$ . By the assumptions on the geometry alteration and Corollary 4.2, $\mathbb{K}(s(t,\cdot ))\in \left (C^2(\bar{\Omega })\right )^{d\times d}$ for every $t\in [0,T_1)$ . Due to (4.5) and Remark 4.5, we obtain $v \in C(\overline{\Omega _T})$ for some $0\lt T\leq T_1$ . By standard Sobolev embedding theorems, also $\nabla \cdot v\in C(\overline{\Omega _T})$ holds true. As such, we ensure the regularity of $v$ as required for coefficients in Theorem A.3 and may therefore proceed with the fixed-point argument analogously to Theorem 3. As mentioned before, the solutions $(v(t,\cdot ),p(t,\cdot ))$ for Darcy’s equation are uniquely determined for each time-slice $t$ . Repeating the arguments following (3.13), uniqueness for the full solution triplet is obtained.

5. Conclusion

In this research, we presented local-in-time existence results to a common class of multi-scale models for reactive transport in evolving porous media. Describing geometry alterations by diffeomorphisms, smooth dependence of the diffusion and permeability tensors on the evolving domain is proven. Being computed as affine-linear functionals of solutions to elliptic and Stokes-type PDEs, the method of transformation to a reference domain in combination with the implicit function theorem was applied for this purpose. Subsequently, we leveraged the resulting smoothness in the micro–macro coupling to short-time existence results to the partially and fully coupled system in the diffusive transport case. Similar results are presented for the partial coupling in the advective case additionally incorporating Darcy’s equation.

As the second major aspect, we proved the induction of suitable paths of diffeomorphisms by the level-set equation. Therefore, our results cover a broad range of possible underlying microscopic geometries. In particular, no explicit construction of adequate families of diffeomorphisms is needed using this methodology.

Further research is needed to derive conditions under which the short-time existence results presented in this paper can be extended to global solutions. This especially concerns the adequate treatment of topological changes in the underlying geometry which is not approachable using diffeomorphisms. As such, the methods used in this research are in particular unable to analyse the behaviour of solutions in scenarios involving clogging or the complete dissolution of the porous structure.

In addition, further work is required to extend our results to local-in-time existence for the fully coupled model including advective transport. More precisely, methods need to be refined in order to compensate for the drop of regularity between the permeability tensor field and the advective velocity field. Finally, future effort is required to generalise our findings to the multi-solute and multi-mineral case, leading to a system of possibly non-linearly coupled parabolic equations on the macroscopic scale.