1. Introduction
Multiscale techniques are known as a powerful tool to model reactive transport in porous media [Reference Molins and Knabner27]. 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 [Reference Čanić6, Reference Seigneur, Mayer and Steefel37]. For many of the aforementioned applications, the structure of the porous medium is altered. Therefore, such problems inherently belong to the class of freeboundary problems, i.e. moving boundary problems with apriori unknown interface evolution. Moving/freeboundary problems enjoy continuous mathematical concern and development [Reference Čanić6, Reference Chen, Shahgholian and Vazquez8]. An overview of methods to analyse moving boundary problems such as energy methods or tools from operator and semigroup theory is found in [Reference Prüss and Simonett33]. Multiscale models are typically derived by means of volume averaging [Reference Whitaker45] or homogenisation [Reference Hornung19] and can be extended to the situation of evolving microstructures [Reference Molins and Knabner27] in the context of reactive transport.
In twoscale models for reactive transport in evolving porous media such as derived in [Reference Bringedal, von Wolff and Pop5, Reference Van Noorden44] 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 twoway 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 twoway coupling between both scales will be referred to as full coupling. Commonly, also a simplified coupling structure is investigated in the literature, disregarding the backcoupling from the macro to the microscale, cf. Figure 1. We refer to the arising onesided coupling as a partial coupling scenario.
An additional challenge is posed by a suitable framework to capture the evolving geometries. A variety of convenient methods is outlined in [Reference Alber1] regarding the description of evolving microstructures. Commonly, levelset methods or phasefield approaches are used, in which case the macroscopic concentrations either prescribe a normal interface velocity, cf. [Reference Van Noorden44], or induce a source term to the phase field [Reference Bringedal, von Wolff and Pop5], 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 [Reference Schulz and Knabner35, Reference Schulz, Ray, Frank, Mahato and Knabner36] or the thickness of the precipitation layer as in [Reference Van Noorden43]. Furthermore, models can also handle the geometry evolution implicitly by computing porosity from conservation of mass and deducing all geometryrelated effective quantities therefrom by heuristic laws, cf. [Reference Jambhekar, Mejri, Schröder, Helmig and Shokri20, Reference Soulaine and Tchelepi41] facilitating the numerical treatment. In [Reference Gärttner, Frolkovic, Knabner and Ray14, Reference Olivares, Bringedal and Pop30] micro–macro models including fully resolved microscopic geometries have been investigated numerically. Regarding the most general setting, the associated evolution PDE (levelset or phasefield 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 (phasefield) models is relatively straightforward as the phasefield 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 sharpinterface limit, i.e. blow up for diffusive interface width tending to zero, cf. [Reference Olivares, Bringedal and Pop30] 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 [Reference Gahn, NeussRadu and Pop12, Reference Peter32] 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 twoscale convergence. Typically, the upscaling process of the transformed model is of increased complexity due to the appearance of additional factors in the highestorder terms arising from the transformation itself. Restricting to diffusion–reaction systems, [Reference Peter31] 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 porespace volume are reflected. Likewise, diffusion–advection–reaction equations are treated in [Reference Evans11] assuming an apriori 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 diffusiondriven model with partial and full coupling between the scales was shown in [Reference McLean24] 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 [Reference Ray and Schulz34], fully coupled systems with diffusiondriven transport are investigated. Furthermore, [Reference Prüss and Simonett33] considered partially coupled systems with advection including degenerating hydrodynamic parameters. Yet, a key assumption is the apriori 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 localintime solutions to the partially coupled model potentially including advective solute transport. On the other hand, we treat the scenario of full coupling and diffusiondriven transport. In that case, the macroscopic concentrations are coupled to the full levelset equation for geometry evolution. Thereby, we extend similar results known in the literature for simple and restrictive geometries like cubes or spheres, cf. [Reference Ray and Schulz34]. Moreover, we discuss sufficient conditions for the longterm 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 localintime 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 [Reference Valent42] where its derivation from a detailed porescale model by formal homogenisation arguments was performed in two spatial dimensions. A generalisation to three dimensions was derived in [Reference Ray and Schulz34] 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\}$ :
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:
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
with Kronecker delta $\delta _{ij}$ , where $\zeta _j$ are the solutions to the following elliptic problems:
with outer unit normal $\nu =\nu (t,x,y)$ . For the derivation in the context of evolving geometries, see [Reference Valent42]. In this case, the evolution is reflected in timedependent 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 [Reference Hornung18], 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:
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:
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 divergencefree velocity field in (2.4) is a common simplification as discussed in [Reference Gärttner, Frolkovic, Knabner and Ray14]. Due to the much larger timescale of geometry evolution compared to fluid flow, it is justified to disregard the flow induced by fluid displacement arising from a variable porespace volume.
The permeability tensor in (2.4) is defined as
where $(\omega _j, \pi _j)$ are the solutions to the Stokestype problems, cf. [Reference Valent42]:
Likewise, the shape of the Stokestype problem is identical to the one derived under the assumption of fixed underlying geometries in [Reference Hornung18]. Again, the evolving geometry is reflected in the timedependent domains. We refer to (2.6) as permeability cell problems. Note that both cell problems (2.3), (2.6) result in symmetric positive semidefinite tensors $\mathbb{D}$ , $\mathbb{K}$ .
Finally, the levelset method is used to describe the evolving underlying geometries [Reference Seigneur, Mayer and Steefel37]. Therefore, we assume the existence of a levelset function $\Phi _0\,:\,\Omega \times Y \to \mathbb{R}$ characterising the solid part within the unitcell 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 zerolevelset 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 levelset equation for $\Phi \,:\,(0, T)\times \Omega \times Y \to \mathbb{R}$ , cf. [Reference Seigneur, Mayer and Steefel37]:
Due to the finite propagation speed of information in (2.7), cf. [Reference Eden, Nikolopoulos and Muntean10], we refrain from posing explicit boundary conditions on $(0, T)\times \Omega \times \partial Y$ here, see also Section 3.4. The different subdomains 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 levelset function uniquely determines porosity $\phi$ and specific surface $\sigma$ as required in (2.1) via
using the ddimensional Hausdorff measure $\mathcal{H}^{d}$ . The interface velocity is coupled to the chemical reaction in a mass conserving way, e.g.
potentially using a scalar speed modification function $v_{\text{mod}}$ which allows for a varying normal interface velocity within a unitcell, cf. [Reference Galdi13]. 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 macroscale 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 levelset equation in Section 3.4. By doing so, we establish regular dependence of the required effective quantities $\phi,\sigma,\mathbb{D}$ on the levelset equation’s time variable, enabling the application of fixedpoint arguments. This ultimately leads to localintime 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.
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 righthand side, cf. [Reference Ray and Schulz34]. Therefore, we write
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 goodnatured conditions such as the smoothly bounded solid geometry being compactly contained within the unitcell $Y$ . By considering localintime 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 [Reference Geißert, Heck and Hieber16]. 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 affinelinear 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 [Reference Chen, Shahgholian and Vazquez8]. 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 Yperiodic 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:
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
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 pointwise 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. [Reference Chen, Shahgholian and Vazquez8].
Let us now suppose that $\partial ^{\text{int}}\mathcal{P}$ is $C^{k+2,1}$ regular. By the surjectivity of the higherorder trace operator continuously extending the following function, cf. [Reference Molins and Knabner27],
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 [Reference Meier25], the required interior higher regularity is established. Following the lines of Theorem 4 of [Reference Meier25], 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 pointwise almosteverywhere 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 porescale 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 pointwise 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 [Reference Geißert, Heck and Hieber16], 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 nonuniformly periodic or evolving domains in the context of porous media [Reference Evans11, Reference Henry17, Reference Olivares, Bringedal and Pop30] or shape optimisation minimising energy functionals depending on a PDE’s solution [Reference Shamir39]. In order to redefine 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.
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
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. [Reference Geißert, Heck and Hieber16, Reference Lee23]. 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.
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:
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 welldefined 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 unitcell 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:
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. [Reference Peter32], 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 oneonone 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
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
derived in [Reference Geißert, Heck and Hieber16] using the inverse transposed Jacobian matrix $(\nabla h)^{T}$ , we compute
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
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 [Reference Soulaine and Tchelepi41] 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
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
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 1parametric 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
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 convexcombination 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
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 1parametric 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 [Reference Cioranescu and Donato9, Reference Ray and Schulz34].
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
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 [Reference Soulaine and Tchelepi41].
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
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 realvalued order parameter $s$ can be replaced by its vectorvalued analogue, allowing for more sophisticated couplings and more complex geometries. A natural field of application is posed by twomineralphase solids, cf. [Reference Gärttner, Frolkovic, Knabner and Ray14], 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 localintime solutions to (3.1) under the assumption of no backcoupling from the macro to the microscale, cf. Figure 1. That is, we assume the evolution of the underlying pore geometry to be known apriori. 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 [Reference Ray and Schulz34], a similar fully coupled problem was solved under the assumption of effective parameters being parameterisable by the porosity $\phi$ in a smooth and apriori 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 [Reference Ray and Schulz34] 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 timeparameter $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 nondegenerative case, i.e. we assume
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 [Reference Ray and Schulz34]. Finally, we state the following anisotropic Sobolev spaces:
with $r\gt d+2$ which play a crucial role in the subsequent existence result. Following the major steps of [Reference Ray and Schulz34], we have
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 porespace 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:
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}^{d1}\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
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
which translates to the continuity of $\sigma$ by using a partition of unity subordinate to the domains of parameterisation, cf. [Reference Ladyzhenskaya, Solonnikov and Ural‘ceva22]. As a result of the continuity of all effective parameters with respect to $s$ , there exists a $\delta \in (0,1)$ such that:
Note that the continuity of the eigenvalues is inherited from the continuity of $\mathbb{D}$ , cf. [Reference Jambhekar, Mejri, Schröder, Helmig and Shokri20]. Following the technique presented in [Reference Ray and Schulz34], the proof is now based on Schauder’s fixedpoint theorem applied to the set
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 fixedpoint form:
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 righthand side of (3.10). Using the compact embeddings
$\mathcal{F}_1$ shows to be compact, cf. [Reference Ray and Schulz34]. 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 aboveprescribed initial and boundary conditions with coefficients defined as
and source term $f$ according to the righthand side of (3.10).
As such, for a given $\tilde{c}\in \mathcal{X}_1$ , we have
abbreviating the contributions from initial and boundary data by
Due to $f$ being locally Lipschitz, there exists a monotone function $\tilde{f}$ satisfying
Similar to the estimates established in [Reference Ray and Schulz34], we obtain using Hölder’s inequality and (3.12)
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 selfmapping 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 [Reference Ray and Schulz34] by considering two solutions $c_1,c_2\in \mathcal{X}_1$ as well as their difference $\bar{c}=c_2c_1$ . Subtracting both associated equations in fixedpoint form (3.10), we obtain
Testing (3.13) with $\bar{c}$ and estimating the righthand side with Hölder’s and Young’s inequality show
where $L$ denotes the Lipschitz constant of $f$ with respect to the compact interval
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. Levelset 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 levelset 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 levelset equation which is a much more natural setup from the viewpoint of applications. At first, we must fix the class of realvalued functions whose level sets are guaranteed to be smooth submanifolds of codimension one:
Definition 3. Regular levelset 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 levelset 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 levelset function, cf. [Reference Geißert, Heck and Hieber16] Chapter 1 or Theorem A.1 .
Remark 3.9. Levelset solution concepts. We note that viscosity solutions pose a powerful and popular weak solution concept for a wide range of nonlinear firstorder PDEs like the levelset equation, for a basic introduction to this topic, see e.g. [Reference Barles4]. 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 nonlinear levelset 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 nonlinear firstorder PDEs, we rewrite the levelset 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
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
with initial conditions
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 zerolevelset describing the position of the fluid–solid interfaces is transported along $x$ . Note that every admissible point of the characteristic ODE system is noncharacteristic, 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 localintime $C^2$ solution by standard theory [Reference Eden, Nikolopoulos and Muntean10]. 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 levelset solution. Let a normal velocity field $v_n\in C^{4,1}(\bar{Y},\mathbb{R})$ be given. Furthermore, let a regular levelset 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 localintime solution to the levelset 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
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
is welldefined 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
has a zerolevelset 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 nonconvex 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 timedependent macroscopic concentration field. Accordingly, this chapter is dedicated to the localintime existence of strong solutions to a model where the geometry evolution is described by the levelset equation (2.7). For clarity, we restate the model equations:
supplemented by boundary and initial conditions as well as diffusion cell problem (2.3), where the levelset normal interface velocity is given by
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 [Reference Ray and Schulz34] 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
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 localintime 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 levelset 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 levelset equation. Using the properties of (3.18), we can reparametrise all effective parameters as a function of $\phi$ by setting:
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 yindependent 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
$\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}^{d1}\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
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
to be of class $C^2$ , which translates to $\sigma (s)$ by using a partition of unity, cf. [Reference Ladyzhenskaya, Solonnikov and Ural‘ceva22]. 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}^*)$ levelset function evolving according to the levelset 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 localintime 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$ :
Proof. In the setup presented above, Theorem A.2 ensures the shorttime existence of solutions $(c,\phi )\in (\mathcal{X}_1)^2$ to the equations
with $\hat{\mathbb{D}},\; \hat{\sigma }$ as constructed in (3.19). Apparently, the second equation corresponds to the levelset function $\Phi$ fulfilling
for $y$ sufficiently close to the zerolevelset 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 zerolevelset 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 [Reference Ray and Schulz34].
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, globalintime existence cannot be expected in general due to the requirement of classical solutions to the nonlinear hyperbolic levelset 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 wellbehavedness of the underlying microscopic geometry evolution.
Remark 3.12. Uniqueness. Note that a proper levelset 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 [Reference Ray and Schulz34], the solution $(c,\phi )\in \mathcal{X}_1^2$ is unique. By the oneonone 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 localintime existence results for the diffusive–advective transport case with partial oneway microtomacro coupling in Section 4.4. A graphical representation of the roadmap for this section is found in Figure 6.
4.1. Setting
In the following, we consider solute transport by diffusion and advection:
which is coupled to Darcy’s equation
In comparison to the equations stated in Section 2, we allow for a general constantintime 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 apriori 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. [Reference Gahn, NeussRadu and Pop12]:
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 [Reference Gahn, NeussRadu and Pop12], we furthermore introduce the space
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
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
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
by Lemma IV 1.1 in [Reference Gahn, NeussRadu and Pop12]. Furthermore, Theorem IV 4.1 in [Reference Gahn, NeussRadu and Pop12] 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 [Reference Gahn, NeussRadu and Pop12]. Next, we consider the problem for general inhomogeneities $j$ of the specified class. According to Theorem 3.4 in [Reference Gärttner, Frolkovic, Knabner and Ray15], 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:
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
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
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
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
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 firstorder 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 timeslices of the spacetime 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 flowchannel 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}^{d1} (\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 )$
The associated velocity field is then given as
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 [Reference Čanić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_2p_1$ . Using suitable estimates on the new righthand 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}}$ :
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:
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 wellestablished 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$ :
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
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_2p_1$ leads to
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
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 counterexample in a setting of maximal smoothness see [Reference Sethian38].
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 [Reference Eden, Nikolopoulos and Muntean10]. 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
Proof. Using standard elliptic theory, the solutions $p_1,p_2$ to (4.2) are uniformly bounded in $W^{3,2}(\Omega )$ , cf. [Reference Eden, Nikolopoulos and Muntean10]. By the Gagliardo–Nirenberg interpolation theorem [Reference Necas, Simader, Necasová, Tronel and Kufner28] we have for $p=p_2p_1$
Using the uniform boundedness of $D^3p_{L^2(\Omega )}$ and Lemma 4.3, we see
which translates to the statement claimed.
Remark 4.5. Continuity on spacetime cylinder. In case $\mathbb{K}\,:\,\Omega _T\to \mathbb{R}^{d\times d}$ fulfils the conditions of Corollary 4.4 for each timeslice 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 spacetime cylinder.
4.4. Existence for partial coupling
Analogously to Section 3.3, we use our prior results to prove localintime 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 redefine the functions $\phi,\; \sigma,\; \mathbb{D},\; \mathbb{K}$ from mappings on $\Omega _T$ , cf. (3.1), to mappings of the realvalued order parameter without change of notation. In order to avoid additional difficulties due to degeneracy, we require
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 porespace 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 fixedpoint 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 timeslice $t$ . Repeating the arguments following (3.13), uniqueness for the full solution triplet is obtained.
5. Conclusion
In this research, we presented localintime existence results to a common class of multiscale 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 affinelinear functionals of solutions to elliptic and Stokestype 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 shorttime 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 levelset 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 shorttime 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 localintime 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 multisolute and multimineral case, leading to a system of possibly nonlinearly coupled parabolic equations on the macroscopic scale.
Acknowledgements
We further acknowledge the insightful discussions with Helmut Abels, University of Regensburg.
Financial support
This research was supported by the DFG Research Training Group 2339 Interfaces, Complex Structures, and Singular Limits. Nadja Ray was also supported by the DFG Research Unit 2170 MadSoil.
Competing interests
The authors have no relevant financial or nonfinancial interests to disclose.
Appendix A
Theorem A.1. (Tubular neighbourhood theorem, adapted from Theorem 1.5, [Reference Geißert, Heck and Hieber16]) Let $\Omega \subset \mathbb{R}^d$ open have a $C^{m,\alpha }$ regular boundary, $2\leq m\leq \infty$ . There exists $r\gt 0$ so that if
then $t(\cdot )\,:\,B_r(\partial \Omega ) \to ({}r,r), \; \pi (\cdot )\,:\, B_r(\partial \Omega )\to \partial \Omega$ are well defined, $\pi$ is a $C^{m1,\alpha }$ retraction onto $\partial \Omega$ ( $\pi (x)=x$ when $x\in \partial \Omega$ ) and t has the same smoothness as $\partial \Omega$ . Further
is a $C^{m1,\alpha }$ diffeomorphism with inverse
$t(\cdot )$ is the unique solution to $\nabla t(x)=1$ in $B_r(\partial \Omega )$ . The largest choice of $r$ is $r=1/\max k$ , where $k$ is the sectional curvature of the boundary in any (tangent) direction at any point of $\partial \Omega$ .
Theorem A.2. (Localintime existence of strong solutions, Theorem 4.1, [Reference Ray and Schulz34]) Let the system of PDEs be given by
Let $\Omega \subset \mathbb{R}^d$ , $d\in \{2,3\}$ , be a domain with $C^2$ smooth boundary $\partial \Omega$ , $r\gt d+2$ , $c_0\in W^{2\frac{2}{r},r}(\Omega )$ , $c_0 \geq 0$ , satisfying the compatibility condition $c_{0\partial \Omega } \equiv 0$ and let $\phi _0 \in W^{2,r}(\Omega )$ hold with $\phi _0(x)\in (\delta, 1\delta ) \subset (0,1)$ for all $x \in \Omega$ and some $\delta \in (0,\frac{1}{2})$ . Furthermore, let $D\in C^1((0,1))$ be a positive scalar function and $\sigma \in C((0,1))$ , $\tau \in C^2((0,1))$ . Then, there exists a constant $T\gt 0$ and at least one strong solution $(c,\phi ) \in \mathcal{X}_1^2$ solving ( A.1 ) with
Theorem A.3. (Parabolic regularity, specialised form of Theorem 9.1, Chapter IV [Reference Kato21]) Let a parabolic problem be given by
with the uniformly parabolic operator $\mathcal{L}$ in nondivergence form
Let $r\gt d+2$ . Suppose that the coefficients $a_{i,j}$ of the operator $\mathcal{L}$ are bounded and continuous in $\Omega _T$ , while the coefficients $a_i$ and $a$ have finite norms $a_i_{L^r(\Omega _T)}$ and $a_{L^r(\Omega _T)}$ . Furthermore, let $\Omega$ be a bounded domain of class $C^2$ . Then for any $f\in L^r(\Omega _T)$ , Dirichlet data $U_0 \in W^{1\frac{1}{2r}, 2\frac{1}{r}}_r(\partial \Omega _T)$ and $u_0 \in W^{2\frac{2}{r},r}(\Omega )$ being compatible in the sense of $U_0(0,\cdot )=u_0$ on $\partial \Omega$ , problem ( A.2 ) has a unique solution $u\in W^{1,2}_r(\Omega _T)$ satisfying the apriori estimate
Theorem A.4. (Differentiability of an implicit equation solution, Theorem 6 [Reference Shamir39]) We give us

• an open set $\mathcal{U}$ in a Banach space $U$ , $u_0 \in \mathcal{U}$ , two reflexive Banach spaces $A$ and $B$ ,

• a map $F\,:\,\mathcal{U}\times A \to B$ , such that $F(u;\cdot )\in \mathcal{L}(A;B)$ for all $u\in \mathcal{U}$ ,

• a function $m\,:\,\mathcal{U}\to A$ , and a function $f\,:\,\mathcal{U}\to B$ , such that
\begin{align*} F(u,m(u))=f(u) \quad \forall u\in \mathcal{U}. \end{align*}

1. Assume that $u \mapsto F(u;\cdot )$ is differentiable at $u_0$ into $\mathcal{L}(A;B)$ , $f$ is differentiable at $u_0$ ,
\begin{align*} F(u_0;x)_B \geq \alpha x_A \quad \forall x\in A, \quad \quad \text{for some } \alpha \gt 0. \end{align*}Then the map $u\mapsto m(u)$ is differentiable at $u_0$ . Its derivative $m^\prime (u_0;\cdot )$ is the unique solution of\begin{align*} F(u_0;m^\prime (u_0;v)) = f^\prime (u_0;v)\partial _u F(u_0;m(u_0);v)\quad \forall v \in U. \end{align*} 
2. In addition, assume that for some integer $k\geq 1$ ,
\begin{align*} u \mapsto F(u;\cdot ) \text{ and } f \text{ are } k \text{ times differentiable at } u_0. \end{align*}Then, the map $u\mapsto m(u)$ is $k$ times differentiable at $u_0$ .