# Strongly Coupled Two-scale System with Nonlinear Dispersion:
Weak Solvability and Numerical Simulation

Vishnu<EMAIL_ADDRESS>Surendra Nepal  Rainey Lyons
Michael Eden  Adrian Muntean

###### Abstract

We investigate a two-scale system featuring an upscaled parabolic dispersion-
reaction equation intimately linked to a family of elliptic cell problems. The
system is strongly coupled through a dispersion tensor, which depends on the
solutions to the cell problems, and via the cell problems themselves, where
the solution of the parabolic problem interacts nonlinearly with the drift
term. This particular mathematical structure is motivated by a rigorously
derived upscaled reaction-diffusion-convection model that describes the
evolution of a population of interacting particles pushed by a large drift
through an array of periodically placed obstacles (i.e., through a regular
porous medium).

We prove the existence of weak solutions to our system by means of an
iterative scheme, where particular care is needed to ensure the uniform
positivity of the dispersion tensor. Additionally, we use finite element-based
approximations for the same iteration scheme to perform multiple simulation
studies. Finally, we highlight how the choice of micro-geometry (building the
regular porous medium) and of the nonlinear drift coupling affects the
macroscopic dispersion of particles.

Key words: Two-scale system; Nonlinear dispersion; Weak solutions; Iterative
scheme; Simulation;
MSC2020: 35G55; 35A01; 35M30; 47J25; 65M60.

##  1 Introduction

The transport of substances through porous media typically involves a
combination of drift, diffusion, and adsorption processes, all of which take
place at the scale of the heterogenous porous structure. The interplay between
these physical processes and the geometry of the underlying heterogeneous
structure is intricate. Only in rare circumstances can one obtain useful
insight into the quantitative understanding of the effective dispersion
mechanisms (responsible for an eventual macroscopic migration of substances),
the overall storage capacity of a given heterogeneous medium, or the resulting
turbulent diffusion and anomalous forms of dissipation; see e.g. [32] for a
discussion of simple turbulent flows conveying microstructure information.
However, understanding these phenomena is crucial for a number of modern
technological applications, including the design of drug delivery systems,
remediation of groundwater contamination, design of efficient hydrogen storage
systems, creation of polymer-based morphologies to facilitate current
transport for organic solar cells; we refer the reader to [18, 26, 13] for
examples of such applications.

In this work, we are primarily interested in understanding what effects are
observable at a macroscopic level, i.e. in terms of the effective indicators,
within the interplay between fast drift and diffusion. In other words, we are
revisiting the Taylor-Aris notion of dispersion now taking place within a
specific type of heterogeneous medium222G. I. Taylor tells in his 1954 paper
[43] that shear flow smears out the concentration distribution, enhancing the
rate at which it spreads in the direction of the flow. R. Aris confirms that
the effect takes place at large Peclet numbers(see [8]).. As a first step, we
are neither looking at processes related to sorption/adsorption mechanisms nor
to surface chemical reactions, even though we are aware that such interface
processes take place in real materials and can have an effect on the
macroscopic response. We are, of course, not the first ones to look at Taylor-
Aris dispersion for spatially structured domains. The majority of the previous
approaches aimed to justify that the macroscopic transport mechanisms can be
effectively modeled via a constant dispersion tensor which encodes both the
geometric properties of the underlying microstructure (like porosity and
connectivity) and the diffusivity of the solute. Such effective models for
transport in porous media are usually derived via some type of scale analysis;
this can be the reference elementary volume (REV) approach often utilized in
the engineering community [11, 12, 45], or the more rigorous route of
mathematical homogenization [4, 6, 7]. In these works, the resulting limit
problem is structurally identical to a diffusion problem and the constant
symmetric effective dispersion can be directly calculated via cell problems.
However, in some situations, this perspective is too simplistic as the
physical reality is more complicated. For instance, earlier investigations
have shown that the resulting dispersion tensor may not be symmetric [10] and
that it can depend (linearly or nonlinearly) on the solute concentration [41].
Such effects cannot be captured in the averaging aforementioned approaches.

To this end, we are considering a dispersion-reaction equation governing the
effective solute transport and production in a porous medium where the
dispersion depends nonlinearly on the solute concentration,

$\partial_{t}u+\text{div}(-D^{*}(W)\nabla u)=f\quad\text{in}\
(0,T)\times\Omega,$ (1a) over some macroscopic domain
$\Omega\subset\mathbb{R}^{2}$ representing a porous medium. While Equation 1a
is macroscopic, the dispersion tensor $D^{*}(W)$ depends on microscopic
solutions $W=(w_{1},w_{2})$ of cell problems over the fluid part of the porous
microstruture $Y$ via
$\text{div}_{y}\left(-D\nabla_{y}w_{i}+G_{i}(u)Bw_{i}\right)=\text{div}_{y}\left(De_{i}\right)\quad\text{in}\
Y,\quad(i=1,2).$ (1b)

Here, $D$ denotes the solute diffusivity, $B$ the velocity field representing
the movement of the solvent, and $G_{i}(u)$ the drift coupling. Via the
product $G_{i}(u)Bw_{i}$, the macroscopic function $u$ nonlinearly interacts
with the cell problem solution $W$. Without this specific interaction in the
cell problems, i.e., taking the case $G_{i}(u)\equiv const.$, the dispersion
tensor assumes the standard form of linear dispersion, see [6]. Moreover, the
case $G_{i}(u)\equiv 0$ corresponds to pure diffusion [2].

The structure of our model is motivated by the homogenization limit
$\varepsilon\to 0$ of a fast drift problem of the form

$\partial_{t}u^{\varepsilon}+\mathrm{div}(-D^{\varepsilon}\nabla
u^{\varepsilon}+\frac{1}{\varepsilon}B^{\varepsilon}P(u^{\varepsilon}))=f^{\varepsilon},$
(2)

where $0<\varepsilon\ll 1$ is a small scale parameter representative of the
porous geometry which is assumed to be periodic; we refer to [40] for more
details regarding this limiting process. The $\varepsilon$-dependent problem
with nonlinear fast drift given by Equation 2 is associated with a
hydrodynamic limit of a totally asymmetric simple exclusion process (TASEP)
which describes a population of interacting particles crossing a domain with
obstacles [17, 16, 38]. In this scenario, the nonlinearity takes the form
$P(r):=r(1-r)$ which corresponds to $G_{i}(r)=1-2r$ in Equation 1b (see [40]).
Since nonlinear interactions in the diffusion/dispersion part come up quite
often in upscaled models, e.g., [31, 5, 24, 14], we study the case where the
terms $G_{i}(u)$ are general functions which satisfy certain assumptions
presented later. It is worth noting that our scientific questions are very
much in the spirit of [22], where the authors ask: what is the macroscopic
response of

$\partial_{t}u^{\varepsilon}+\mathrm{div}(-\eta\nabla
u^{\varepsilon}+u^{\varepsilon}V^{\varepsilon})=f^{\varepsilon},$ (3)

where $\eta>0$ and $V^{\varepsilon}$ is some large random field? In other
words, what precisely is the term $\langle uV\rangle$ in the corresponding
upscaled equation

$\partial_{t}\langle u\rangle+\mathrm{div}(-\eta\nabla\langle u\rangle+\langle
uV\rangle)=\langle f\rangle,$ (4)

where $\langle\cdot\rangle$ denotes some suitable average, and how can this be
computed numerically. Comparing (3) and (4) with (2) and (1a), we see that the
vector field $V^{\varepsilon}$ is in our setting what the TASEP requires, and
hence, the meaning of $\langle uV\rangle$ is for us precisely $-D^{*}(W)\nabla
u$.

System (1) is an example of a two-scale system, where we refer to the
dispersion equation (1a) as the macroscopic equation and the cell problems
(1b) as the microscopic equations. Such systems are also referred to as
distributed-microstructure models, terminology introduced by R. E. Showalter
in [42]. The presence of nonlinear coupling in the dispersion tensor has been
previously studied. For instance, in [19], the authors looked at the nonlinear
effect of colloids deposition on diffusion in porous media and established
local-in-time existence (see also [33, 35]). In the same context, in [44] a
similar model was analyzed via homogenization techniques. Additionally, a
related homogenized model for reactive flows in porous media was derived in
[5].

The simulation of such scale separated problems often come with unique
computational challenges. For instance, the nonlinear drift coupling in the
cell problems requires solving Equation 1b many times. Well-established
numerical techniques capable of handling such multiscale features are reported
in [3, 23, 1, 15]. Some of the closer works that explore the simulation of
two-scale systems similar to ours are [27, 36, 35]. In [27], the authors
address a two-scale coupled problem where the coupling is via the reaction
rate of macroscopic equation. On the other hand, authors in [36] study and
simulate a two-scale system for a phase-field model for precipitation and
dissolution in porous media via an iterative scheme and, in [35], clogging of
a porous medium via aggregation of colloidal particles was explored.

This paper is structured as follows: In Section 2, we present the mathematical
model of interest, the geometric setup, as well as the assumptions on
coefficients and data needed for the mathematical analysis. This is followed
by Section 3 where the main result, namely the existence of solutions (Theorem
2), is proven via an iteration scheme. In Section 4, we use this iteration
scheme to numerically simulate different scenarios, showcasing that the
proposed two-scale model is able to capture interesting effects at both the
microscopic and macroscopic scale. The results reported in this section are
only preliminary, as the numerical analysis and further numerical exploration
of our system are postponed for follow-up work. Finally in Section 5, we
conclude our work with a short summary of the findings and with a discussion
of possible future investigations.

##  2 Setting of the problem

Let $T>0$ represent the time horizon. We denote by
$\Omega\subset\mathbb{R}^{2}$ a non-empty open bounded domain with Lipschitz
boundary, $\partial\Omega$, and by $Y_{0}\subset(0,1)^{2}$ a compact set with
positive measure and Lipschitz boundary such that the microscopic domain
$Y:=(0,1)^{2}\backslash Y_{0}$ is connected. We let $t\in(0,T)$ represent the
time variable, $x\in\Omega$ the macroscopic space variable, $y\in Y$ the
microscopic space variable, $\Gamma_{N}=\partial Y_{0}$, and $n_{y}$ the
outward unit normal vector across the interface $\partial Y$. With the sets
$\Omega$ and $Y$ at hand, we can build a two-scale geometry. For a schematic
representation of such a geometry, see Figure 1. We take $e_{i}$, $i=1,2$, to
denote the standard unit basis vectors in $\mathbb{R}^{2}$.

$x$$\Omega$$\Gamma_{N}$$Y$ Figure 1: Typical two-scale geometry: schematic
representation of the macroscopic domain $\Omega$ and of the microscopic
domain $Y$ with internal boundary $\Gamma_{N}$.

We study the following nonlinear two-scale system:

ptProblem $P(\Omega)$ Find the unknown functions $u$ and $W:=(w_{1},w_{2})$
satisfying $\displaystyle\partial_{t}u+\mathrm{div}(-D^{*}(W)\nabla u)$
$\displaystyle=f$ in $\displaystyle\hskip 5.69046pt(0,T)\times\Omega,$ (5a)
$\displaystyle u$ $\displaystyle=0$ on $\displaystyle\hskip
8.5359pt(0,T)\times\partial\Omega,$ (5b) $\displaystyle u(0)$
$\displaystyle=g$ in $\displaystyle\hskip 8.5359pt\overline{\Omega},$ (5c)
$\displaystyle\mathrm{div}_{y}\left(-D\nabla_{y}w_{i}+G_{i}(u)Bw_{i}\right)$
$\displaystyle=\mathrm{div}_{y}(De_{i})$ in $\displaystyle\hskip 5.69046ptY,$
(5d) $\displaystyle\left(-D\nabla_{y}w_{i}+BG_{i}(u)w_{i}\right)\cdot n_{y}$
$\displaystyle=\left(De_{i}\right)\cdot n_{y}$ on $\displaystyle\hskip
5.69046pt\Gamma_{N},$ (5e) $\displaystyle w_{i}\,$ $\displaystyle\mbox{ is
$Y$--periodic}.$ (5f)

Relating back to our motivation discussed in Section 1 (cf. also [40, 17]),
$u$ represents the particle concentration on macroscopic porous domain
$\Omega$, while $W$, the solution to cell problems, relays the microscopic
information, including the shape of the domain, drift, and diffusion, to the
macroscopic scale via the dispersion tensor, $D^{*}(W)$, given by

$D^{*}(W):=\frac{1}{|Y|}\int_{Y}D(y)\left(I+\begin{bmatrix}\frac{\partial
w_{1}}{\partial y_{1}}&\frac{\partial w_{2}}{\partial y_{1}}\\\ \frac{\partial
w_{1}}{\partial y_{2}}&\frac{\partial w_{2}}{\partial
y_{2}}\end{bmatrix}\right)\,\mathrm{d}y$ (6)

where $I$ is the $2\times 2$ identity matrix and $|Y|$ is the Lebesgue measure
of the set $Y$. The periodic boundary condition (5f) and the Neumann influx
boundary condition (5e) are common in effective models derived by
homogenization (see [30]).

We assume the reaction rate $f:(0,T)\times\Omega\longrightarrow\mathbb{R}$,
the initial condition $g:\Omega\longrightarrow\mathbb{R}$, the microscopic
diffusion matrix $D:Y\rightarrow\mathbb{R}^{2\times 2}$, the microscopic
velocity field $B:Y\rightarrow\mathbb{R}^{2}$, and the nonlinear functions
$G_{1},G_{2}:\mathbb{R}\rightarrow\mathbb{R}$ are given and satisfy additional
assumptions discussed later. Note that, to calculate the macroscopic equation
and dispersion tensor, we need to, at each point $(t,x)\in(0,T)\times\Omega$,
solve the microscopic problem in $Y$.

To discuss the concept of weak solution to Problem $P(\Omega)$, we introduce
the following functional spaces:

$\displaystyle H_{\\#}^{1}(Y)$ $\displaystyle:=\\{v\in
H^{1}(Y):v\;\mbox{is}\;Y-\mbox{periodic}\\},$ $\displaystyle\mathcal{U}$
$\displaystyle:=\\{v\in L^{2}((0,T);H_{0}^{1}(\Omega)):\partial_{t}v\in
L^{2}((0,T);H^{-1}(\Omega))\\},$ $\displaystyle\text{and }\mathcal{W}$
$\displaystyle:=\left\\{v\in
H_{\\#}^{1}(Y):\int_{Y}v(y)\,\mathrm{d}y=0\right\\},$

where $H_{\\#}^{1}(Y)$ is equipped with the standard $H^{1}(Y)$ norm. We
denote by $\langle\cdot,\cdot\rangle$ the duality pairing between
$H^{-1}(\Omega)$ and $H_{0}^{1}(\Omega).$ Finally, we define the weak solution
of Problem $P(\Omega)$ in the following sense:

######  Definition 2.1.

We say that $(u,W)$ is a weak solution to Problem $P(\Omega)$ if
$u\in\mathcal{U}$ with $u(0,\cdot)=g$ and, for almost every
$(t,x)\in(0,T)\times\Omega$, $w_{i}\in\mathcal{W}$ and the following integral
equations are satisfied:

$\displaystyle\langle\partial_{t}u,\phi\rangle+\int_{\Omega}D^{*}(W)\nabla
u\cdot\nabla\phi\,\mathrm{d}x$
$\displaystyle=\int_{\Omega}f\phi\,\mathrm{d}x,$ (7a)
$\displaystyle\int_{Y}\big{(}D\nabla_{y}w_{i}-G_{i}(u(t,x))Bw_{i}\big{)}\cdot\nabla_{y}\psi\,\mathrm{d}y$
$\displaystyle=\int_{Y}\mathrm{div}_{y}(De_{i})\psi\,\mathrm{d}y-\int_{\Gamma_{N}}De_{i}\cdot
n_{y}\psi\,\mathrm{d}\sigma,$ (7b)

for all $(\phi,\psi)\in H^{1}(\Omega)\times H_{\\#}^{1}(Y)$ and
$i\in\\{1,2\\}$.

#### Assumptions.

From now on $C$ denotes a positive real number, possibly changing its value
from line to line. We consider the following restrictions on data and model
parameters:

  1. (A1)

The microscopic diffusion matrix satisfies $D\in(H^{1}_{\\#}(Y)\cap
L^{\infty}(Y))^{2\times 2}$ and there exists $\theta>0$ such that

$\theta|\eta|^{2}\leq D\eta\cdot\eta\quad\text{for all}\
\eta\in\mathbb{R}^{2}\;\mbox{and almost all }y\in Y;$

  2. (A2)

$G_{1},G_{2}:\mathbb{R}\rightarrow\mathbb{R}$ are locally Lipschitz functions,
i.e., Lipschitz on compact sets;

  3. (A3)

The microscopic drift velocity $B\in(H^{1}_{\\#}(Y)\cap L^{\infty}(Y))^{2}$
satisfies

$\begin{cases}\mathrm{div}_{y}B=0\hskip 11.38092pt\mbox{in}\hskip
11.38092ptY,\\\ B\cdot n_{y}=0\hskip 11.38092pt\mbox{on}\hskip
11.38092pt\Gamma_{N};\end{cases}$

  4. (A4)

The reaction rate satisfies $f\in L^{\infty}((0,T)\times\Omega)$ and the
initial condition $g\in L^{\infty}(\Omega)$.

Assumptions (A1)-(A4) may seem technical but can be satisfied by many
physically relevant functions. For example, Assumption (A2) can be satisfied
by taking $B$ to be the solution to an incompressible Stokes problem with no
penetration into the obstacle and periodic boundary conditions across
$\partial Y\setminus\Gamma_{N}$. This is also the approach we have taken in
our simulations, see Section 4.

We remark that in Problem $P(\Omega)$ we can also have nonhomogeneous
Dirichlet boundary conditions $u=h$ on $\partial\Omega$ for some given $h\in
H^{\frac{1}{2}}(\Omega).$ This regularity allows us to extend $h$ to
$\tilde{h}\in H^{1}(\Omega)$ and transform the original problem into

$\displaystyle\partial_{t}\tilde{u}+\mathrm{div}(-D^{*}(W)\nabla_{x}\tilde{u})$
$\displaystyle=\tilde{f}(W)$ in $\displaystyle\hskip
5.69046pt(0,T)\times\Omega,$ $\displaystyle\tilde{u}$ $\displaystyle=0$ on
$\displaystyle\hskip 8.5359pt(0,T)\times\partial\Omega,$
$\displaystyle\tilde{u}(0)$ $\displaystyle=\tilde{g}$ in $\displaystyle\hskip
8.5359pt\overline{\Omega},$
$\displaystyle\mathrm{div}_{y}\left(-D\nabla_{y}w_{i}+G_{i}(\tilde{u}+\tilde{h})Bw_{i}\right)$
$\displaystyle=\mathrm{div}_{y}(De_{i})$ in $\displaystyle\hskip 5.69046ptY,$
$\displaystyle\left(-D\nabla_{y}w_{i}+BG_{i}(\tilde{u}+\tilde{h})w_{i}\right)\cdot
n_{y}$ $\displaystyle=\left(De_{i}\right)\cdot n_{y}$ on $\displaystyle\hskip
5.69046pt\Gamma_{N},$ $\displaystyle w_{i}\,$ $\displaystyle\mbox{ is
$Y$--periodic},$

where
$\tilde{u}:=u-\tilde{h},\;\tilde{f}(W):=f-\mathrm{div}(D^{*}(W)\nabla\tilde{h})$
and $\tilde{g}:=g-\tilde{h}$.

##  3 Weak solvability of Problem $P(\Omega)$

In this section, we prove the existence of solutions to Problem $P(\Omega)$ in
the sense of Definition 2.1. There are two main difficulties in the analysis
of the involved equations: (i) the nonlinear coupling in the drift, and (ii)
the positivity of the dispersion tensor.

The main strategy of our proof follows four distinct steps:

  * Step 1:

We construct a two-scale iterative scheme based on Problem $P(\Omega)$. This
scheme is constructed by decoupling $P(\Omega)$ into a family of linear
problems, which we refer to as $P^{k}(\Omega)$, $k\in\mathbb{N}$. The notion
of weak solutions $(u^{k},W^{k})$ to Problem $P^{k}(\Omega)$ is given in
Definition 3.1.

  * Step 2:

We construct an auxiliary problem similar to (5d)–(5f) and show it is well-
posed in Lemma 1. This helps us to understand the coupling between
$G_{i}(u^{k})$ and $W^{k}$. We then show in Lemma 3 that the dispersion tensor
is bounded and positive definite, uniformly with respect to the index $k$.

  * Step 3:

With Lemma 1 and Lemma 3, we show that for each $k$ the problems
$P^{k}(\Omega)$ are well-posed, their solutions satisfy $k$-independent energy
estimates, and are also uniformly bounded; see both Theorem 1 and Lemma 4.

  * Step 4:

Finally, we prove our main result (Theorem 2) – the solution to the iterative
scheme converges in a suitable sense to a solution to Problem $P(\Omega)$.

###  3.1 Iterative scheme

We begin by introducing the iterative scheme and its corresponding weak
formulation. We set $u^{0}=g$, and, for any $k\in\mathbb{N}\cup\\{0\\}$, we
denote as $u^{k+1},w^{k}_{1},$ and $w^{k}_{2}$ the solutions to the following
decoupled system:

ptDiffusion-drift problem
$\displaystyle\mathrm{div}_{y}\left(-D\nabla_{y}w_{i}^{k}+G_{i}(u^{k})Bw_{i}^{k}\right)$
$\displaystyle=\mathrm{div}_{y}(De_{i})$ in $\displaystyle\hskip 5.69046ptY,$
(8a)
$\displaystyle\left(-D\nabla_{y}w_{i}^{k}+BG_{i}(u^{k})w_{i}^{k}\right)\cdot
n_{y}$ $\displaystyle=\left(De_{i}\right)\cdot n_{y}$ on $\displaystyle\hskip
5.69046pt\Gamma_{N},$ (8b) $\displaystyle w_{i}^{k}\,$ $\displaystyle\mbox{ is
$Y$--periodic},$ $\displaystyle i\in$ $\displaystyle\\{1,2\\}$ (8c)

ptReaction-dispersion problem
$\displaystyle\partial_{t}u^{k+1}+\mathrm{div}(-D^{*}(W^{k})\nabla_{x}u^{k+1})$
$\displaystyle=f$ in $\displaystyle\hskip 5.69046pt(0,T)\times\Omega,$ (8d)
$\displaystyle u^{k+1}(0)$ $\displaystyle=g$ in $\displaystyle\hskip
8.5359pt\overline{\Omega},$ (8e) $\displaystyle u^{k+1}$ $\displaystyle=0$ on
$\displaystyle\hskip 8.5359pt(0,T)\times\partial\Omega,$ (8f)

where the dispersion tensor $D^{*}(W^{k})$ is given by

ptEffective dispersion tensor
$D^{*}(W^{k}):=\frac{1}{|Y|}\int_{Y}D(y)\left(I+\begin{bmatrix}\frac{\partial
w^{k}_{1}}{\partial y_{1}}&\frac{\partial w^{k}_{2}}{\partial y_{1}}\\\
\frac{\partial w^{k}_{1}}{\partial y_{2}}&\frac{\partial w^{k}_{2}}{\partial
y_{2}}\end{bmatrix}\right)\,\mathrm{d}y.$ (8g)

We refer to the iterative scheme (8a)–(8f) together with (8g), as Problem
$P^{k}(\Omega)$. We define the concept of weak solutions to Problem
$P^{k}(\Omega)$ in a similar manner as for the Problem $P(\Omega)$ (see
Definition 2.1).

######  Definition 3.1.

Given $u^{k}\in\mathcal{U},$ we say that $(u^{k+1},W^{k})$ is a weak solution
to Problem $P^{k}(\Omega)$ if $u^{k+1}\in\mathcal{U}$ with
$u^{k+1}(0,\cdot)=g$ and, for almost every $(t,x)\in(0,T)\times\Omega$,
$w_{i}^{k}\in\mathcal{W}$ and the following integral equations are satisfied:

$\displaystyle\langle\partial_{t}u^{k+1},\phi\rangle+\int_{\Omega}D^{*}(W^{k})\nabla
u^{k+1}\cdot\nabla\phi\,\mathrm{d}x$
$\displaystyle=\int_{\Omega}f\phi\,\mathrm{d}x$ (9a)
$\displaystyle\int_{Y}\left(D\nabla_{y}w^{k}_{i}-G_{i}(u^{k}(t,x))Bw^{k}_{i}\right)\cdot\nabla_{y}\psi\,\mathrm{d}y$
$\displaystyle=\int_{Y}\mathrm{div}_{y}(De_{i})\psi\,\mathrm{d}y-\int_{\Gamma_{N}}De_{i}\cdot
n_{y}\psi\,\mathrm{d}\sigma,$ (9b)

for all $(\phi,\psi)\in H^{1}(\Omega)\times H_{\\#}^{1}(Y)$ and
$i\in\\{1,2\\}$.

###  3.2 Existence results for the Problem $P^{k}(\Omega)$

In this section, we study the well-posedness of the iterative scheme
(8a)–(8g). To analyze the scheme, we first focus on the elliptic equation
(8a)–(8c) by constructing an auxiliary problem similar to (5d)–(5f). We show
this auxiliary problem is well-posed, prove energy estimates, and establish
the regularity of weak solutions. Then with the help of Lemma 3, which
establishes the uniform positive definite property and boundedness of
$D^{*}(W^{k})$, we show our iterative scheme Problem $P^{k}(\Omega)$ is well-
posed.

######  Lemma 1.

Let $p\in\mathbb{R}$. Consider the following auxiliary problem: Find
$W(p,\cdot)=(w_{1}(p,\cdot),w_{2}(p,\cdot))\in\mathcal{W}^{2}$ satisfying

$\displaystyle\mathrm{div}_{y}\left(-D\nabla_{y}w_{i}+pBw_{i}\right)$
$\displaystyle=\mathrm{div}_{y}(De_{i})$ in $\displaystyle\hskip 5.69046ptY,$
(10a) $\displaystyle\left(-D\nabla_{y}w_{i}+pBw_{i}\right)\cdot n_{y}$
$\displaystyle=\left(De_{i}\right)\cdot n_{y}$ on $\displaystyle\hskip
5.69046pt\Gamma_{N},$ (10b) $\displaystyle w_{i}\,$ $\displaystyle\mbox{ is
$Y$--periodic},$ (10c)

where $i\in\\{1,2\\}$. Assume (A1)–(A2) holds, then

  1. (i)

there exists a unique weak solution $w_{i}(p,\cdot)\in\mathcal{W}$ to the
problem (10a)–(10c),

  2. (ii)

there exists a constant $C>0$ independent of $p$ such that

$\|\nabla_{y}w_{i}(p,\cdot)\|_{L^{2}(Y)}\leq C,$ (11)

  3. (iii)

there exists a constant $C>0$ such that, for all $p_{1},p_{2}\in\mathbb{R}$,

$\int_{Y}|\nabla_{y}(w_{i}(p_{1},y)-w_{i}(p_{2},y))|^{2}\,\mathrm{d}y\leq
C|p_{1}-p_{2}|^{2}.$ (12)

###### Proof.

(i) We begin by stating the weak form of (10a)–(10c):

$\int_{Y}\big{(}D\nabla_{y}w_{i}-pBw_{i}\big{)}\cdot\nabla_{y}\psi\,\mathrm{d}y=\int_{Y}\mathrm{div}_{y}(De_{i})\psi\,\mathrm{d}y-\int_{\Gamma_{N}}De_{i}\cdot
n_{y}\psi\,\mathrm{d}\sigma\quad$ (13)

for all $\psi\in H_{\\#}^{1}(Y)$ and $i\in\\{1,2\\}$. From (A1), we have

$\int_{Y}\mathrm{div}_{y}(D(y)e_{i})\,\mathrm{d}y=\int_{\Gamma_{N}}(D(y)e_{i})\cdot
n_{y}\,\mathrm{d}\sigma.$ (14)

Using identity (14) together with standard arguments involving the classical
Fredholm alternative, we have the existence and uniqueness of $w_{1}(p,\cdot)$
and $w_{2}(p,\cdot)$ in the sense of (13) (for details we refer the reader [4,
Section 4]).

(ii) First, using integration by parts, we have the following equation:

$\displaystyle\int_{Y}pBw_{i}\cdot\nabla_{y}w_{i}\,\mathrm{d}y$
$\displaystyle=\frac{p}{2}\int_{Y}B\nabla_{y}(w_{i}^{2})\,\mathrm{d}y$
$\displaystyle=\frac{p}{2}\int_{Y}w_{i}^{2}\mathrm{div}_{y}B\,\mathrm{d}y+\frac{p}{2}\int_{\partial
Y\backslash\Gamma_{N}}B\cdot n_{y}w_{i}^{2}\,\mathrm{d}\sigma$
$\displaystyle\hskip 56.9055pt+\frac{p}{2}\int_{\Gamma_{N}}B\cdot
n_{y}w_{i}^{2}\,\mathrm{d}\sigma.$ (15)

By (A2), we have that $\int_{Y}w_{i}^{2}\mathrm{div}_{y}B\,\mathrm{d}y=0$ and
$\int_{\Gamma_{N}}B\cdot n_{y}w_{i}^{2}\,\mathrm{d}\sigma=0$. Additionally, by
the periodicity of $w_{i}$ and $B$, we have $\int_{\partial
Y\backslash\Gamma_{N}}B\cdot n_{y}w_{i}^{2}\,\mathrm{d}\sigma=0$, and so we
obtain

$\int_{Y}pBw_{i}\cdot\nabla_{y}w_{i}\,\mathrm{d}y=0.$ (16)

Now, taking the test function $\psi=w_{i}$ in the weak formulation (13) and
making use of (16), we obtain

$\int_{Y}D\nabla_{y}w_{i}\cdot\nabla_{y}w_{i}\,\mathrm{d}y=\int_{Y}\mathrm{div}_{y}(De_{i})w_{i}\,\mathrm{d}y-\int_{\Gamma_{N}}\left(De_{i}\right)\cdot
n_{y}w_{i}\,\mathrm{d}\sigma.$ (17)

Making use of (A1) and integration by parts, we see that

$\theta\int_{Y}|\nabla_{y}w_{i}|^{2}\,\mathrm{d}y\leq\int_{Y}|(De_{i})\nabla_{y}w_{i}|\,\mathrm{d}y.$
(18)

Young’s inequality applied to the right-hand side of (18) together with (A1),
yields (11).

(iii) Let

$\overline{w_{i}}:=w_{i}(p_{1},y)-w_{i}(p_{2},y)\text{ and
}\overline{p}:=p_{1}-p_{2}$

and consider the weak formulations (13) with $p=p_{1},\,p_{2}$. We take the
test function $\psi=\overline{w_{i}}$ in both weak forms and subtract the
equations to arrive at

$\int_{Y}\big{(}D(y)\nabla_{y}\overline{w_{i}}-B(y)(p_{1}w_{i}(p_{1},y)-p_{2}w_{i}(p_{2},y)))\big{)}\cdot\nabla_{y}\overline{w_{i}}\,\mathrm{d}y=0.$
(19)

Adding and subtracting the term $\int_{Y}B(y)p_{2}w_{i}(p_{1},y)\,\mathrm{d}y$
to (19) and making use of assumption (A1), we obtain

$\theta\int_{Y}|\nabla_{y}\overline{w_{i}}|^{2}\,\mathrm{d}y\leq\int_{Y}B(y)\big{(}\overline{p}w_{i}(p_{1},y)+p_{2}\overline{w_{i}}\big{)}\cdot\nabla_{y}\overline{w_{i}}\,\mathrm{d}y.$
(20)

Since $w_{1}$ and $w_{2}$ are periodic, $\overline{w_{i}}$ is periodic as
well. Therefore, we may follow similar arguments which lead to (16) to obtain

$\int_{Y}B(y)p_{2}\overline{w_{i}}\cdot\nabla_{y}\overline{w_{i}}\,\mathrm{d}y=0.$
(21)

Returning to (20), using Young’s inequality, we have

$\displaystyle\theta\int_{Y}|\nabla_{y}\overline{w_{i}}|^{2}\,\mathrm{d}y\leq\int_{Y}B(y)\overline{p}w_{i}(p_{1},y)\cdot\nabla_{y}\overline{w_{i}}\,\mathrm{d}y$
$\displaystyle\leq
C\overline{p}\int_{Y}|w_{i}(p_{1},y)||\nabla_{y}\overline{w_{i}}|$
$\displaystyle\leq\frac{1}{2\theta}|\overline{p}|^{2}\|w_{i}(p_{1})\|_{L^{2}(Y)}^{2}+\frac{\theta}{2}\int_{Y}|\nabla_{y}\overline{w_{i}}|^{2}\,\mathrm{d}y.$
(22)

Then, with (11), we arrive at

$\int_{Y}|\nabla_{y}\overline{w_{i}}|^{2}\,\mathrm{d}y\leq
C|\overline{p}|^{2}.$ (23)

∎

To study the weak solvability of (8d)–(8f), we first show that the dispersion
tensor (8g) is both uniformly positive definite and uniformly bounded. To this
end, for any $p\in\mathbb{R}$, we define
$\overline{D}(p)\in\mathbb{R}^{2\times 2}$ by

$\displaystyle[\overline{D}(p)]_{i,j}:=\dfrac{1}{|Y|}\int_{Y}D(y)\left(e_{j}+\nabla_{y}w_{j}(p,y)\right)\cdot
e_{i}\,\mathrm{d}y,$ (24)

where $i,j\in\\{1,2\\}$ and $w_{j}(p,y)$ is the weak solution to (10a)–(10c).
We begin by showing that the matrix $\overline{D}(p)$ is Lipschitz continuous
with respect to $p$.

######  Lemma 2.

Assume (A1)–(A2) hold. There exists a constant $C>0$, independent of $p$ and
$q$, such that

$|[\overline{D}(p)]_{i,j}-[\overline{D}(q)]_{i,j}|\leq C|p-q|.$

###### Proof.

From definition (24), assumption (A1), and Lemma 1 $(iii)$, we have

$\displaystyle|[\overline{D}(p)]_{i,j}-[\overline{D}(q)]_{i,j}|$
$\displaystyle\leq
C\int_{Y}|D(y)|\left|\nabla_{y}(w_{i}(p,y)-w_{i}(q,y))\right|\,\mathrm{d}y$
$\displaystyle\leq
C\left(\int_{Y}|\nabla_{y}(w_{i}(p,y)-w_{i}(q,y))|^{2}\,\mathrm{d}y\right)^{\frac{1}{2}}$
$\displaystyle\leq C|p-q|.$ (25)

∎

We now show that the dispersion tensor is both uniformly bounded and uniformly
positive definite.

######  Lemma 3.

Assume (A1)–(A2) hold and let $M>0$. For $p=(p_{1},p_{2})\in[-M,M]^{2}$, let
$W_{p}:=(w_{1}(p_{1},\cdot),w_{2}(p_{2},\cdot))$ be the corresponding
solutions of the auxiliary problem (10a)–(10c). Then the macroscopic
dispersion tensor $D^{*}(W_{p})$ satisfies the following properties:

  1. (i)

There exists $\theta_{M}>0$ independent of $p$ such that

$\theta_{M}|\eta|^{2}\leq D^{*}(W_{p})\eta\cdot\eta\quad\text{for all}\
\eta\in\mathbb{R}^{2};$ (26)

  2. (ii)

There exist a constant $C>0$ independent of $p$ such that

$|[D^{*}(W_{p})]_{i,j}|\leq C,\quad(i,j\in\\{1,2\\}).$ (27)

###### Proof.

(i) From the definition of $\overline{D}(\cdot)$ (24) and the definition of
$D^{*}(\cdot)$ (6), we have

$[D^{*}(W_{p})]_{i,j}=[\overline{D}(p_{i})]_{i,j},$ (28)

and hence,

$[D^{*}(W_{p})]_{i,j}=\frac{1}{|Y|}\int_{Y}D(y)\left(e_{j}+\nabla_{y}w_{j}(p_{i},y)\right)\cdot
e_{i}\,\mathrm{d}y.$ (29)

We consider the weak formulation for the problem (10a)–(10c) (i.e., Equation
(13)) for index $j$, where we choose the test function $\psi=w_{i}$ and divide
by $|Y|$:

$\frac{1}{|Y|}\int_{Y}\big{(}D\nabla_{y}w_{j}-p_{j}Bw_{j}\big{)}\cdot\nabla_{y}w_{i}\,\mathrm{d}y=\frac{1}{|Y|}\int_{Y}\text{div}_{y}(De_{j})w_{i}\,\mathrm{d}y-\frac{1}{|Y|}\int_{\Gamma_{N}}\left(De_{j}\right)\cdot
n_{y}w_{i}\,\mathrm{d}\sigma.$ (30)

Using integration by parts on
$\int_{Y}\text{div}_{y}(De_{j})w_{i}\,\mathrm{d}y$ and noting the periodicity
of $D(\cdot)$ and $w_{i}(\cdot)$, we arrive at

$\frac{1}{|Y|}\int_{Y}\big{(}D\nabla_{y}w_{j}-p_{j}Bw_{j}+De_{j}\big{)}\cdot\nabla_{y}w_{i}\,\mathrm{d}y=0.$
(31)

Adding (31) to (29), we see that

$[D^{*}(W_{p})]_{i,j}=[A(p)]_{i,j}+[J(p)]_{i,j},$ (32)

where

$\displaystyle[A(p)]_{i,j}$
$\displaystyle:=\frac{1}{|Y|}\int_{Y}D(e_{j}+\nabla_{y}w_{j})\cdot(e_{i}+\nabla_{y}w_{i})\,\mathrm{d}y,$
(33) $\displaystyle\text{and }[J(p)]_{i,j}$
$\displaystyle:=-\frac{1}{|Y|}\int_{Y}p_{j}Bw_{j}\cdot\nabla_{y}w_{i}\,\mathrm{d}y.$
(34)

Using (A2) and integration by parts on
$\int_{Y}p_{j}Bw_{j}\cdot\nabla_{y}w_{i}\,\mathrm{d}y$, we obtain

$\displaystyle-[J(p)]_{i,j}$
$\displaystyle=-\frac{p_{j}}{|Y|}\int_{Y}\nabla_{y}(Bw_{j})w_{i}\,\mathrm{d}y+\frac{p_{j}}{|Y|}\int_{\partial
Y}B\cdot n_{y}w_{i}w_{j}\,\mathrm{d}\sigma$
$\displaystyle=-\frac{1}{|Y|}\int_{Y}p_{j}Bw_{i}\cdot\nabla_{y}w_{j}\,\mathrm{d}y$
$\displaystyle=[J(p)]_{j,i}.$ (35)

Combining (33) and (35) shows that $A(p)$ and $J(p)$ are the symmetric and
skew-symmetric parts of $D^{*}(W_{p})$. Since $J(p)$ is a $2\times 2$ skew-
symmetric matrix, we have

$J(p)\eta\cdot\eta=0\quad\text{for all}\ \eta\in\mathbb{R}^{2}.$ (36)

So, to prove (26) it is enough to establish that $A(p)$ is in fact positive
definite. Since $A(p)$ is symmetric, using the structure (33) we have for all
$\eta=(\eta_{1},\eta_{2})\in\mathbb{R}^{2}$,

$\displaystyle A(p)\eta\cdot\eta$
$\displaystyle=\eta_{1}^{2}\frac{1}{|Y|}\int_{Y}D(e_{1}+\nabla_{y}w_{1})\cdot(e_{1}+\nabla_{y}w_{1})\,\mathrm{d}y$
$\displaystyle\hskip
28.45274pt+2\eta_{1}\eta_{2}\frac{1}{|Y|}\int_{Y}D(e_{1}+\nabla_{y}w_{1})\cdot(e_{2}+\nabla_{y}w_{2})\,\mathrm{d}y$
$\displaystyle\hskip
42.67912pt+\eta_{2}^{2}\frac{1}{|Y|}\int_{Y}D(e_{2}+\nabla_{y}w_{2})\cdot(e_{2}+\nabla_{y}w_{2})\,\mathrm{d}y.$

Distributing $\eta_{1}$ and $\eta_{2}$ inside the integrals, we get

$\displaystyle A(p)\eta\cdot\eta$
$\displaystyle=\frac{1}{|Y|}\int_{Y}D(\begin{bmatrix}\eta_{1}\\\
0\end{bmatrix}+\eta_{1}\nabla_{y}w_{1})\cdot(\begin{bmatrix}\eta_{1}\\\
0\end{bmatrix}+\eta_{1}\nabla_{y}w_{1})\,\mathrm{d}y$ $\displaystyle\hskip
28.45274pt+2\frac{1}{|Y|}\int_{Y}D(\begin{bmatrix}\eta_{1}\\\
0\end{bmatrix}+\eta_{1}\nabla_{y}w_{1})\cdot(\begin{bmatrix}0\\\
\eta_{2}\end{bmatrix}+\eta_{2}\nabla_{y}w_{2})\,\mathrm{d}y$
$\displaystyle\hskip 42.67912pt+\frac{1}{|Y|}\int_{Y}D(\begin{bmatrix}0\\\
\eta_{2}\end{bmatrix}+\eta_{2}\nabla_{y}w_{2})\cdot(\begin{bmatrix}0\\\
\eta_{2}\end{bmatrix}+\eta_{2}\nabla_{y}w_{2})\,\mathrm{d}y.$ (37)

Rearranging the right side of the identity (37) and using (A1) yields

$\displaystyle A(p)\eta\cdot\eta$
$\displaystyle=\frac{1}{|Y|}\int_{Y}D(\begin{bmatrix}\eta_{1}\\\
\eta_{2}\end{bmatrix}+\sum_{i=1}^{2}\eta_{i}\nabla_{y}w_{i})\cdot(\begin{bmatrix}\eta_{1}\\\
\eta_{2}\end{bmatrix}+\sum_{i=1}^{2}\eta_{i}\nabla_{y}w_{i})\,\mathrm{d}y$
$\displaystyle\geq\theta\frac{1}{|Y|}\int_{Y}|\eta+\sum_{i=1}^{2}\eta_{i}\nabla_{y}w_{i})|^{2}\;\,\mathrm{d}y$
$\displaystyle\geq 0.$ (38)

If $A(p)\eta\cdot\eta=0$ for some $\eta\in\mathbb{R}^{2}$, then we have

$\displaystyle|\eta+\sum_{i=1}^{2}\eta_{i}\nabla_{y}w_{i})|^{2}=0\,\,\,\mbox{for
almost every }\,y\in Y.$

Consequently, this yields

$\sum_{i=1}^{2}\eta_{i}w_{i}=C-\eta\cdot y,$ (39)

for some constant $C$. Since the $w_{i}$ are periodic, identity (39) is only
possible for $\eta=0$. Hence, $A(p)\eta\cdot\eta=0$ only for $\eta=0$.
Therefore, combining this with (36), we get

$D^{*}(W_{p})\eta\cdot\eta>0.$ (40)

for all $p\in[-M,M]^{2}$ and all $\eta\neq 0$. Define the function
$F:[-M,M]^{2}\times\partial B(0,1)\rightarrow\mathbb{R}$ as

$F(p,\xi):=D^{*}(W_{p})\xi\cdot\xi\,\,\,\mbox{for all}\;\xi\in\partial
B(0,1).$

We have from Lemma 2 that the function $p_{i}\mapsto\overline{D}(p_{i})$ is
continuous. Using the definition (28), we see that the function $p\mapsto
D^{*}(W_{p})$ is also continuous. Therefore, $F$ is continuous and
$F(p,\xi)>0$ over the compact set $([-M,M]^{2}\times\partial B(0,1))$. Hence,
there exist a $\theta_{M}>0$ such that

$F(p,\xi)>\theta_{M}$

for all $(p,\xi)\in([-M,M]^{2}\times\partial B(0,1))$. So, for any $\eta\neq
0$ we have

$\displaystyle D^{*}(W_{p})\frac{\eta}{|\eta|}\cdot\frac{\eta}{|\eta|}$
$\displaystyle>\theta_{M},$ $\displaystyle\text{i.e.,
}D^{*}(W_{p})\eta\cdot\eta$ $\displaystyle>\theta_{M}|\eta|^{2}.$

(ii) The inequality (27) follows directly from (A1), $(ii)$ of Lemma 1, and
the identity

$[D^{*}(W_{p})]_{i,j}=\frac{1}{|Y|}\int_{Y}D\left(e_{j}+\nabla_{y}w_{j}\right)\cdot
e_{i}\,\mathrm{d}y$

together with the Cauchy–Schwartz inequality. ∎

######  Theorem 1 (Well-posedness of $P^{k}(\Omega)$).

Assume (A1)–(A4) hold and let $u^{0}=g$. Then, there exist a sequence
$(u^{k},W^{k})_{k\in\mathbb{N}\cup\\{0\\}}\subset\mathcal{U}\times
L^{\infty}((0,T)\times\Omega;\mathcal{W}^{2})$ such that, for each $k$,
$(u^{k+1},W^{k})$ uniquely solves the iterative scheme given by (8a)–(8g) in
the sense of Definition 3.1. Moreover, we have

$\|u^{k}\|_{L^{\infty}((0,T)\times\Omega)}\leq\|g\|_{L^{\infty}(\Omega)}+T\|f\|_{L^{\infty}((0,T)\times\Omega)}$
(41)

for all $k\in\mathbb{N}$.

###### Proof.

We proof the existence of solutions via induction. By construction, we have
$u^{0}=g$ and the linear elliptic problems (8a)–(8c) have unique solutions
$w_{i}^{0}(p,\cdot)\in\mathcal{W}$ where $p=G_{i}(u^{0}(t,x))$ at almost all
$(t,x)\in(0,T)\times\Omega$ via Lemma 1 (i). Owing to Lemma 1 (ii), we have
that $\|\nabla_{y}w_{i}(p,\cdot)\|_{L^{2}(\Omega)}\leq C$ uniformly in $p$.
This implies $\|w_{i}(p,\cdot)\|_{\mathcal{W}}\leq C$ via the
Poincaré–Wirtinger inequality and, as a consequence, also
$W^{0}=(w_{1}^{0},w_{2}^{0})\in
L^{\infty}((0,T)\times\Omega;\mathcal{W}^{2})$.

Let

$m:=\|g\|_{L^{\infty}(\Omega)}+T\|f\|_{L^{\infty}((0,T)\times\Omega)}$

and set

$M:=\max_{r\in[-m,m]}\max_{i\in\\{1,2\\}}|G_{i}(r)|$

which is well-defined as the $G_{i}$ are continuous. Since $u^{0}=g$, we have
$|G_{i}(u^{0})|\leq M$. Using Lemma 3, the dispersion tensor $D^{*}(W^{0})$
satisfies

$D^{*}(W^{0})\eta\cdot\eta\geq\theta_{M}|\eta|^{2}$

for all $\eta\in\mathbb{R}^{2}$ almost everywhere in $(0,T)\times\Omega$,
where $\theta_{M}>0$ depends only on $M$.

The problem (8d)–(8f) is a standard linear parabolic problem. Benefiting from
Assumptions (A1)–(A4) and the uniform positivity of $D^{*}(W^{0})$, the
existence of a unique weak solution $u^{1}\in\mathcal{U}$ to the problem
(8d)–(8f) with $k=0$ can be shown via standard Galerkin approximation
arguments (see, e.g., [20, Chapter 7, Theorem 3]). The estimate

$\|u^{1}\|_{L^{\infty}((0,T)\times\Omega)}\leq\|g\|_{L^{\infty}(\Omega)}+T\|f\|_{L^{\infty}((0,T)\times\Omega)},$
(42)

follows by an application of Duhamel’s principle [37, Chapter 5], for details
see [19, Lemma 10].

The solvability of the problem (8a)–(8f) for any $k\in\mathbb{N}$ follows by
induction. Indeed, assume $u^{k}\in\mathcal{U}$ satisfies $u^{k}(0,\cdot)=g$
and

$\|u^{k}\|_{L^{\infty}((0,T)\times\Omega)}\leq m$ (43)

implying

$G_{i}(u^{k}(t,x))\in[-M,M],\quad i=1,2.$ (44)

Then the existence result for $W^{k}\in
L^{\infty}((0,T)\times\Omega;\mathcal{W}^{2})$ solving (8a)–(8c) follows by
Lemma 1 in the same way as for $W^{0}$. Relying on Lemma (3) and (44), we
obtain that $D^{*}(W^{k})$ is bounded and for all $\eta\in\mathbb{R}^{2}$
satisfies

$D^{*}(W^{k})\eta\cdot\eta>\theta_{M}|\eta|^{2}$

almost everywhere in $(0,T)\times\Omega$, where $\theta_{M}>0$ does not depend
on $k$. Following the same arguments as before, we get the existence and
uniqueness of $u^{k+1}\in\mathcal{U}$. The estimate

$\|u^{k+1}\|_{L^{\infty}((0,T)\times\Omega)}\leq\|g\|_{L^{\infty}(\Omega)}+T\|f\|_{L^{\infty}((0,T)\times\Omega)},$

follows similarly to (42). ∎

###  3.3 Convergence of the iterative scheme

Now, we show that the sequence $(u^{k},W^{k})_{k\in\mathbb{N}\cup\\{0\\}}$
satisfies uniform energy estimates. These estimates are crucial in
establishing the convergence of our scheme. We also present a proposition
which will help us pass to the limit in the weak formulation.

######  Lemma 4.

Assume (A1)–(A4) hold and let $u^{k}\in\mathcal{U}$ be weak solution to the
iterative scheme (8a)–(8f). Then, there exists a constant $C>0$ independent of
$k$ such that

$\displaystyle\|u^{k}\|_{L^{\infty}((0,T);L^{2}(\Omega))}$ $\displaystyle\leq
C,$ (45a) $\displaystyle\|\nabla u^{k}\|_{L^{2}((0,T);L^{2}(\Omega))}$
$\displaystyle\leq C,$ (45b)
$\displaystyle\|\partial_{t}u^{k}\|_{L^{2}((0,T);H^{-1}(\Omega))}$
$\displaystyle\leq C.$ (45c)

###### Proof.

Using $(i)$ of Lemma 3 and (41), there exist $\theta_{M}>0$ independent of
$k$, depending on $M$ such that

$D^{*}(W^{k})\eta\cdot\eta\geq\theta_{M}|\eta|^{2}$

for all $\eta\in\mathbb{R}^{2}$. Estimates (45a)–(45c) then follow the
arguments made in [20, Chapter 7, Theorem 2]. ∎

######  Proposition 1.

Let $\xi_{k},\psi_{k},\xi,\psi\in L^{2}(\Omega)$ for all $k\in\mathbb{N}$.
Assume

$\displaystyle\xi_{k}$ $\displaystyle\rightarrow\xi\hskip 28.45274pt$ strongly
in $\displaystyle L^{2}(\Omega),$ (46) $\displaystyle\psi_{k}$
$\displaystyle\rightharpoonup\psi$ weakly in $\displaystyle L^{2}(\Omega)$
(47)

and there exist a $C>0$ independent of $k$, such that

$\|\xi_{k}\|_{L^{\infty}(\Omega)}\leq C,$ (48)

then for any $\phi\in L^{2}(\Omega)$, upto a subsequence we have

$\lim_{k\rightarrow\infty}\int_{\Omega}\xi_{k}\psi_{k}\phi\,\mathrm{d}x\rightarrow\int_{\Omega}\xi\psi\phi\,\mathrm{d}x.$
(49)

###### Proof.

Using (46) and (47), we get

$\xi_{k}\psi_{k}\;\rightharpoonup\;\xi\psi\;\;\;\mbox{weakly
in}\;\;L^{1}(\Omega).$

Hence for all $\phi\in C_{c}^{\infty}({\Omega})$, we have

$\lim_{k\rightarrow\infty}\int_{\Omega}\xi_{k}\psi_{k}\phi\,\mathrm{d}x\rightarrow\int_{\Omega}\xi\psi\phi\,\mathrm{d}x.$
(50)

Now, from (47) and (48), there exist a constant $C>0$ independent of $k$, such
that

$\|\xi_{k}\psi_{k}\|_{L^{2}(\Omega)}\leq C.$ (51)

So, by the weak compactness theorem, we get $\xi_{k}\psi_{k}$ weakly converges
to some $h\in L^{2}(\Omega)$ in $L^{2}(\Omega)$. Now, combining with (50), we
get $h=\xi\psi$ for a.e. $x\in\Omega$, which give us the required result (49).
∎

######  Theorem 2 (Solvability of $P(\Omega)$).

Assume (A1)–(A4) hold. Then, there exists a $u\in\mathcal{U}$ and a $W\in
L^{\infty}((0,T)\times\Omega;\mathcal{W}^{2})$ such that, along a subsequence,

$\displaystyle u^{k}$ $\displaystyle\rightarrow u\hskip 28.45274pt$ strongly
in $\displaystyle L^{2}((0,T)\times\Omega),$ (52a) $\displaystyle\nabla u^{k}$
$\displaystyle\rightharpoonup\nabla u\hskip 28.45274pt$ weakly in
$\displaystyle L^{2}((0,T)\times\Omega),$ (52b)
$\displaystyle\partial_{t}u^{k}$
$\displaystyle\rightharpoonup\partial_{t}u\hskip 28.45274pt$ weakly in
$\displaystyle L^{2}((0,T);H^{-1}(\Omega)),$ (52c) $\displaystyle
D^{*}(W^{k})$ $\displaystyle\rightarrow D^{*}(W)\hskip 28.45274pt$ strongly in
$\displaystyle L^{2}((0,T)\times\Omega).$ (52d) Moreover, $(u,W)$ is a weak
solution to the nonlinear parabolic-elliptic system (5a)–(5f) in the sense of
Definition 2.1 and $u\in L^{\infty}((0,T)\times\Omega)$.

###### Proof.

Using (45a)–(45c) together with Aubin-Lions Compactness Lemma (see [28] for
details or [27, Lemma 4]), we get (52a) for a subsequence. Using the weak
compactness property (see [20, Appendix D, Theorem 3]) with the energy
estimates (45b) and (45c), we get both (52b) and (52c). From (52a), it also
follows that $u^{k}\rightarrow u$ pointwise almost everywhere in
$(0,T)\times\Omega$, possibly along a further subsequence (not relabeled).
Hence, we obtain $u\in L^{\infty}((0,T)\times\Omega)$ and, more precisely,
$|u|\leq m$. Via Lemma 1, we get the existence of the weak solution
$w_{i}(p_{i},\cdot)\in\mathcal{W}$ for (5d)–(5f) where $p_{i}=G_{i}(u(t,x))$,
$i\in\\{1,2\\}$, for almost all $(t,x)\in(0,T)\times\Omega$. This then extends
to a function $W=(w_{1},w_{1})\in
L^{\infty}((0,T)\times\Omega;\mathcal{W}^{2})$ via Lemma 1 (ii).

Recalling that the functions $G_{i}$ are locally Lipschitz (Assumption (A1)),
i.e.,

$|G_{i}(u^{k})-G_{i}(u)|\leq C|u^{k-1}-u|,$ (53)

where $C>0$ is independent of $k$ as $|u^{k}|,|u|\leq m$, we can use the
estimates (25) along with (52a) to show

$\displaystyle\lim_{k\rightarrow\infty}\int_{0}^{T}\int_{\Omega}|[D^{*}(W^{k})-D^{*}(W)]_{i,j}|^{2}\,\mathrm{d}x$
$\displaystyle=\lim_{k\rightarrow\infty}\int_{0}^{T}\int_{\Omega}|[\overline{D}(G_{i}(u^{k-1})))-\overline{D}(G_{i}(u))]_{i,j}|^{2}\,\mathrm{d}x$
$\displaystyle\leq
C\lim_{k\rightarrow\infty}\int_{0}^{T}\int_{\Omega}|G_{i}(u^{k-1})-G_{i}(u)|^{2}\,\mathrm{d}x$
$\displaystyle\leq
C\lim_{k\rightarrow\infty}\int_{0}^{T}\int_{\Omega}|u^{k-1}-u|^{2}\,\mathrm{d}x$
$\displaystyle=0.$ (54)

Using (52b) and (52d) together with Proposition 1, for any $\phi\in
H_{0}^{1}(\Omega)$ we have

$\lim_{k\rightarrow\infty}\int_{\Omega}D^{*}(W^{k})\nabla
u^{k+1}\cdot\nabla\phi\,\mathrm{d}x=\int_{\Omega}D^{*}(W)\nabla
u\cdot\nabla\phi\,\mathrm{d}x.$ (55)

Combining (52c) and (55) and then passing $k\rightarrow\infty$ in (9a), we see
that $u$ satisfies (7a). Hence, we have that $(u,W)$ solves Problem
$P(\Omega)$ in the sense of Definition 2.1. ∎

We remark that the convergence of our scheme (8a)–(8f) does not guarantee the
uniqueness of solutions to Problem $P(\Omega)$. However, under the assumption
$\nabla u\in{L^{\infty}((0,T)\times\Omega)},$ uniqueness follows. Indeed, if
$u_{1},u_{2}$ solve Problem $P(\Omega)$, then by using the Lipschitz
continuity of $D^{*}(W)$ (see Lemma (3)) and standard energy estimates for
$u_{1}-u_{2}$, we can conclude that $u_{1}=u_{2}$ for almost every
$(t,x)\in(0,T)\times\Omega$.

##  4 Numerical simulation

In this section, we demonstrate that the iteration scheme (8) can be used as a
numerical method to approximate the solution of Problem $P(\Omega)$.
Additionally, we use this numerical method to illustrate the impact of
different microscopic geometries and to investigate how the effective
dispersion tensor, $D^{*}$, governs the macroscopic behaviour of the system.
To accomplish this, we make use of finite element method solvers in FEniCS
[29, 25] with Lagrange polynomials of degree one as the basis elements for
both the elliptic and parabolic problems. We begin by discretizing the
macroscopic and microscopic domains. By first fixing an initial iteration
$u^{0}$, we can then solve the linear elliptic cell problems given by
(8a)–(8c) in the microscopic domain for each point in the macroscopic domain.
Please note that this step can be perfectly parallelized as each cell problem
is independent. By calculating $D^{*}(W^{k})$ given by (8g), we can then solve
the linear parabolic problem (8d)–(8f), where we make use of an implicit Euler
discretization to handle the time dependence. Then, using this numerical
solution, we repeat these calculations for the the next iteration. More
concisely, given the numerical solution $u^{k}$, we calculate $u^{k+1}$ by
solving the linear system (8). We continue this iteration method until a
desired error tolerance between successive iterations (measured in the
$L^{2}(0,T;L^{2}(\Omega))$ norm) is reached. We outline this approach in
Algorithm 1.

To proceed with our investigation, we define two specific geometries
(henceforth named geometry 1 and 2) which have different and interesting
effects on the effective dispersion tensor. For geometry 1, we choose
$Y=(0,1)^{2}\setminus\mathcal{B}_{0.25}((0.5,0.5))$ where
$\mathcal{B}_{r}((y_{1},y_{2}))$ denotes the closed disk with radius $r$ and
center $(y_{1},y_{2})$ (for a visualization, see the top row of Figure 2). For
geometry 2, we let $Y=(0,1)^{2}\setminus(\mathcal{R}_{1}\cup\mathcal{R}_{2})$
with rectangles $\mathcal{R}_{1}:=[0.1,0.9]\times[0.1,0.2]$ and
$\mathcal{R}_{2}:=[0.1,0.9]\times[0.8,0.9]$(for a visualization, see the
bottom row of Figure 2).

When choosing model ingredients, one difficulty is ensuring that the
microscopic drift, $B(y)$, is both interesting, physically relevant, and
satisfies the assumptions (A2) for each geometry. To accomplish this, we take
the velocity field as the solution to the following Stokes problem:

ptStokes problem for the drift velocity $\displaystyle-\mu\Delta B+\nabla p$
$\displaystyle=F(y)\quad$ $\displaystyle\text{in}\;\;Y,$ (56a)
$\displaystyle\text{div}\,B$ $\displaystyle=0\;\;$
$\displaystyle\text{in}\;\;Y,$ (56b) $\displaystyle B$ $\displaystyle=0\;\;$
$\displaystyle\text{on}\;\;\Gamma_{N},$ (56c) $\displaystyle y$
$\displaystyle\mapsto B(y)$ $\displaystyle\text{is}\;Y\text{-periodic},$ (56d)

where $Y$ is either geometry 1 or geometry 2. To ensure stability while
solving (56), we use the Taylor-Hood elements in FEniCS, where second-order
polynomials are used as basis function for the velocity field and first-order
polynomials are used for the pressure. We refer the reader, for instance, to
[9] for more information on the use of stable finite elements for the
computation of Stokes flow.

For all simulations to follow, we choose $\mu=0.01$ and

$F:Y\to\mathbb{R}^{2}\mbox{ given by }F(y)=(10\sin(2\pi y_{1})\sin(2\pi
y_{2}),10\sin(2\pi y_{1})\cos(2\pi y_{2}))$

for $y=(y_{1},y_{2})\in Y$. The velocity vector $B(y)=(B_{1}(y),B_{2}(y))$
corresponding to geometry 1 and geometry 2 is shown in Figure 2.

Figure 2:  Top: solution to the Stokes problem, $B_{1}(y)$ (left) and
$B_{2}(y)$ (right), in geometry 1. Bottom: solution to the Stokes problem,
$B_{1}(y)$ (left) and $B_{2}(y)$ (right), in geometry 2. Algorithm 1 Procedure
to compute the weak solution to (8a)–(8f).

1:Discretize the space micro domain $Y$ and macro domain $\Omega$

2:Discretize the time domain $[0,T]$ with step size $\Delta t$

3:Solve the Stokes problem (56a)-(56d) to get $B(y)$

4:Set initial iteration guess $u^{0}$

5:Choose data $f,D,G_{i},g$

6:Set tolerance value $\epsilon$

7:Set the maximum number of iterations, Maxiter.

8:Initialize iteration and time. i.e. $iter=0$, $t=0$

9:while $iter<Maxiter$ do

10: Set $u_{\text{old}}=g$

11: Set $u=[u_{\text{old}}]$

12: for each time discrete node on time domain do

13: for each node on macroscopic grid do

14: Solve for $(w_{1},w_{2})$ using
$(G_{1}(u_{\text{old}}),G_{2}(u_{\text{old}}))$

15: end for

16: Compute $D^{*}$ from $(w_{1},w_{2})$

17: Solve for $u_{\text{new}}$ using $D^{*}$

18: append $u$ with $u_{\text{new}}$

19: $u_{\text{old}}\leftarrow u_{\text{new}}$

20: end for

21: if $\|u-u^{0}\|<\epsilon$ then

22: Stop

23: end if

24: $u^{0}\leftarrow u$

25:end while

###  4.1 Dispersion tensor

We now study the behaviour of the effective dispersion tensor with respect to
the nonlinear drift interactions in both geometries. To begin with, we divide
the interval $[-10,10]$ into $101$ equidistant nodes denoted $p_{i}$. We then
solve the auxiliary cell problem given by (10a)-(10c) for each $p_{i}$. To
enforce the fact that the cell solution has zero average, we use the Lagrange
multiplier method in a similar manner as, for instance, in [21]. With the
solutions to the cell problems, we compute the entries of the effective
dispersion tensor using equation (6). These results are then interpolated and
the components of $D^{*}$ are plotted as a function of $p\in[-10,10]$ in
Figures 3 and 4.

#### Case 1: Fast diffusion

We first consider the case where the microscopic diffusion is comparable in
magnitude to the velocity field $B(y)$. To this end, we choose the diffusion
matrix

$\displaystyle D(y)=\begin{pmatrix}2+\sin(\pi y_{1})\sin(\pi y_{2})&0\\\
0&2+\sin(\pi y_{1})\end{pmatrix},\quad y=(y_{1},y_{2})\in Y,$

which satisfies $1\leq\det D(y)\leq 9$ everywhere in $Y$.

Figure 3: Comparison of the components of $D^{*}(W)$ for the fast diffusion
case with both geometry 1 (solid red line) and geometry 2 (dashed green line).
Here, we plot the components $D^{*}_{11}$ (top left), $D^{*}_{12}$ (top
right), $D^{*}_{21}$ (bottom left), and $D^{*}_{22}$ (bottom right) as
functions of the parameter $p$.

Looking at Figure 3, we notice that the main diagonal entries $D^{*}_{11}$ and
$D^{*}_{22}$ are symmetric with respect to the vertical axis while the off-
diagonal entries $D^{*}_{12}$ and $D^{*}_{21}$ are symmetric about the origin.
However, we note that the off-diagonal entries are close to $0$ and have
therefore only a small impact on the macroscopic dispersion. This is a
consequence of both geometries being symmetric: in non-symmetric setups, e.g.,
angled rectangles or ellipses, the off-diagonal entries would play a bigger
role. The particular case $p=0$ corresponds to “no drift” and, in both
geometries, it is the value with the slowest dispersion.

We can also see that while $D^{*}_{11}$ and $D^{*}_{22}$ are similar (both
quantitatively and qualitatively) in the first geometry, they differ quite a
bit in the second geometry. Focusing on geometry 2, the dispersion in $y_{1}$
direction appears to be much faster than in the $y_{2}$ direction as can be
seen by $D^{*}_{11}$ being almost twice as large as $D^{*}_{22}$. This is to
be expected since the geometric setup with the two rectangles impedes vertical
flows much more than lateral flows. More interestingly, $D^{*}_{11}$ is close
to constant in $p$ while, on the other hand, $D^{*}_{22}$ is much more dynamic
with changes of up to $50\%$ relative to the lowest value for $p=0$. This
shows that the drift interaction can play a large role in facilitating
vertical flows for geometry 2.

#### Case 2: Slow diffusion

Now, we consider the slow diffusion case where the microscopic diffusion is
small relative to the velocity field. We compute the components of the
effective dispersion tensor and compare the results in both geometries in
Figure 4. For this case, we choose the diffusion tensor to be

$\displaystyle D(y):=\begin{pmatrix}0.05+\frac{1}{50}\sin(\pi y_{1})\sin(\pi
y_{2})&0\\\ 0&0.05+\frac{1}{50}\sin(\pi y_{1})\end{pmatrix},\quad
y=(y_{1},y_{2})\in Y,$

which satisfies $0.009\leq\det D(y)\leq 0.0049$ everywhere in $Y$.

Figure 4: Comparison of the components of $D^{*}(W)$ for the slow diffusion
case with both geometry 1 (solid red line) and geometry 2 (dashed green line).
Here, we plot the components $D^{*}_{11}$ (top left), $D^{*}_{12}$ (top
right), $D^{*}_{21}$ (bottom left), and $D^{*}_{22}$ (bottom right) as
functions of the parameter $p$.

Comparing Figure 4 with Figure 3, we can observe very similar trends
qualitatively although the values are, of course, much lower due to the slow
diffusion. At first glance, it appears that there are cusps forming at the
critical point $p=0$ in $D^{*}_{11}$ for geometry 1 and $D^{*}_{22}$ for both
geometries. However, a closer look shows that the transition is in fact
smooth, see Figure 5. These observed “cusps” do however point to the fact that
the slow diffusion case is much more volatile with respect to small changes
close to the zero drift case $p=0$ when compared to the fast diffusion case.
Also, the relative changes in the dispersion values are higher than in the
case of fast diffusion, e.g., in geometry 2 the value $D^{*}_{22}$ experiences
changes of more than 100$\%$ relative to the lowest value for $p=0$.

Figure 5: Comparison of two components of $D^{*}(W)$, $D^{*}_{11}$ (left) and
$D^{*}_{22}$ (right), in the neighbourhood of $p=0$ for the slow diffusion
case.

###  4.2 Macroscopic solution

The goal of this section is to illustrate how the behavior of the dispersion
tensor presented above affects the solution to the macroscopic equation. In
this section, we fix $\Omega=(0,1)\times(0,2)$ and the initial profile of
concentration

$\displaystyle
g(x_{1},x_{2})=\begin{cases}\exp{(-10((x_{1}-0.5)^{2}+(x_{2}-0.5)^{2}))},\;\;\text{if}\;\;(x_{1},x_{2})\in\mathcal{B}_{0.25}((0.5,0.5)),\\\
0,\;\;\text{otherwise}.\end{cases}$ (57)

We also fix a constant source term $f:\Omega\to\mathbb{R}$ given by

$\displaystyle
f(x_{1},x_{2})=\begin{cases}1000,\;\;\;\text{if}\;\;\;(x_{1},x_{2})\in\mathcal{B}_{0.25}((0.5,0.5)),\\\
0,\;\;\;\text{otherwise}.\end{cases}$ (58)

The iteration scheme requires an initial guess $u^{0}$ for which we choose
$u^{0}(t,x)=g$. For the microscopic ingredients, we take the same diffusion
matrix as in the fast diffusion case above for all simulations. The velocity
field, $B(y)$, is treated in the same manner as before.

As explained previously, we iterate the FEniCS solvers until the maximum
number of iterations is achieved or the error
$e_{k}:=||u^{k+1}-u^{k}||_{L^{2}(0,T,L^{2}(\Omega))}<\mbox{tol}$, where tol is
a small tolerance value prescribed a priori and $k\in\mathbb{N}$ is the
iteration index. For all simulations, we choose ${\rm tol}=10^{-7}$ and use
$50$ space nodes in both directions in the macroscopic domain $\Omega$.
Additionally, we take $G_{1}(u)=G_{2}(u)=G(u)$ for these simulations.

#### Effect of the microscopic geometries

In Figure 6, heat maps of the concentration profile $u$ at time $T=2$ with
$G(u)=1-2u$ using both microscopic geometries are displayed. Comparing the
first and second plots in Figure 6, similar dispersive behavior of the
concentration profiles is exhibited in both geometries. However, the
concentration profile corresponding to geometry 1 disperses faster than that
corresponding to geometry 2. This is expected as, due to the long horizontal
rectangular obstacles in geometry 2, the flow is hindered in the vertical
direction. The third plot in Figure 6 shows the difference between both
solutions corresponding to the two different microscopic geometries. Looking
at this difference, we observe that the concentration in the neighborhood of
the point $(0.5,0.5)$ is higher for geometry 2 than for geometry 1. To
quantify this observed behavior, we calculate the total amount of the mass
$M(t)$ of concentration $u$ over time on the upper subdomain,
$(0,1)\times(1,2)$, using the following formula:

$\displaystyle
M(t):=\int_{1}^{2}\int_{0}^{1}u(t,x_{1},x_{2})\,\mathrm{d}x_{1}\,\mathrm{d}x_{2}.$
(59)

We then plot the indicator $M(t)$ over time $t$ for both choices of geometry
in the left plot of Figure 8.

Figure 6:  Macro-solution with $G(u)=1-2u$ corresponding to microscopic
geometry 1 (left), geometry 2 (middle), and the difference between these two
solutions (right) at $T=2$.

#### Effect of the nonlinear drift

We now investigate how the nonlinear drift interactions affect the macroscopic
dispersion. Here, all parameters are kept the same as in the previous
numerical experiment, with the exception of the new nonlinearity
$\tilde{G}(u)=1/(0.0001+|1-2u|)$ which we compare with the previous
$G(u)=1-2u$. The particular form of $\tilde{G}$ is chosen to simulate the
effect of high particle concentration impeding the flow. This can be seen by
looking at Figure 3, where low values of $|G(u)|$ or $\tilde{G}(u)$ correspond
to slightly slower dispersion than high values.

Figure 7: Macro-solution corresponding to geometry 2 with $G(u)=1-2u$ (left),
$\tilde{G}(u)=1/(0.0001+|1-2u|)$ (middle), and the difference between these
two solutions (right) at $T=2$.

The macroscopic solution for both nonlinearities corresponding to geometry 2
at $T=2$ is plotted in Figure 7. Similar to before, the concentration diffuses
slower with the nonlinearity, $\tilde{G}(u)$, than with $G(u)$. Comparing the
first and second plots in Figure 7, we observe that the solution profile is
more concentrated around the source in the case of $\tilde{G}(u)$ (middle)
than $G(u)$ (left). These effects can easily be seen in the right-most plot in
Figure 7 where we calculate the difference between the concentration profiles.
Similar observations can be made in the case of geometry 1, however, with
these choices of ingredients, the difference is less drastic, and so, we omit
these plots. As before, we examine the evolution of the mass in the upper
subdomain, $M(t)$, and plot this in Figure 8 on the right.

Figure 8: Comparison of mass $M(t)$ for geometry 1 (solid red line) and
geometry 2 (dashed green line) with $G(u)=1-2u$ (left) and
$\tilde{G}(u)=1/(0.0001+|1-2u|)$ (right).

##  5 Conclusion and outlook

Understanding the effect microscopic fast drifts have on the structure of the
macroscopic dispersion is generally a difficult task, even in the case of
simple microstructures underlying a regular porous medium. We have shown in
this paper that solutions to our model (Problem $P(\Omega)$) not only exist,
but can also capture subtle two-scale effects of potential interest for
applications. We will explore in follow-up works the numerical analysis of our
setting as well as more numerical simulations of relevant practical case
studies.

By constructing a convergent iterative scheme to Problem $P(\Omega)$, we have
shown the existence of weak solutions, as well as the well-posedness of linear
iterative approximations to the original two-scale system. We have implemented
the constructed iterative scheme using finite element methods and numerically
illustrated some potential effects the microscopic geometry and the drift
coupling have on the macroscopic dispersion tensor. Such a posteriori effects
can be seen either by examining the different profiles of the macroscopic
solutions, or by following the time and space distribution of values in the
entries of the effective dispersion tensor. We strongly believe that the
working methodology proposed here can be used to approach a larger class of
two-scale problems, covering among others the settings discussed, for
instance, in [31, 5, 24, 39].

While the iteration scheme is a useful analytical tool, it is computationally
expensive as the microscopic problem is solved for every time and space node
in each iteration. However, given suitable assumptions about the choice of
nonlinearity $G$, the dispersion tensor can be precomputed by solving the
auxiliary problem (10a)–(10c) over a range of values for the parameter $p$ in
a similar way to the computations in Section 4.1. Additionally, this
precomputing can be parallelized, as the microscopic cell problems are
independent, and save a considerable amount of computation time. In light of
this, we plan to study how this parallelization and precomputing technique
affects the numerical results in terms of both the computation time and the
order of accuracy. Of course, this requires us to additionally study the fully
discrete error analysis of the Galerkin approximation (cf. e.g. [27, 34]) in
order to determine the order of convergence in space and time.

Within the frame of this manuscript, we considered that the macroscopic
problem is defined in a two-dimensional space since our motivation originated
from a hydrodynamic limit performed for an interacting particles system
(TASEP) posed in two dimensions. So, it is natural to take
$\Omega\subset\mathbb{R}^{2}$. However, without much additional difficulty, we
can extend our problem and the corresponding results to dimension $3$ and
higher. The only notable change is that, if Problem $P(\Omega)$ is posed in
$\Omega\times Y\subset\mathbb{R}^{n}\times\mathbb{R}^{d}$, then there will be
$d$ systems of nonlinear elliptic cell problems associated with the dispersion
tensor. Additionally, the mathematical analysis work is not restricted to
Dirichlet boundary conditions. Other choices of boundary conditions (e.g.
homogeneous Neumann) are allowed as well. Furthermore, the macroscopic
equation does not need to be scalar. A system of macroscopic equations coupled
correctly with a family of elliptic cell problems can be handled as well in a
similar fashion.

It is worth noting that we are currently unable to prove the uniqueness of
solutions to Problem $P(\Omega)$ without assuming higher regularity. The main
technical issue seems to be connected to the type of the nonlinearity arising
in the structure of the dispersion tensor. Difficulties in proving the
uniqueness of weak solutions to two-scale systems were also mentioned in [19].
There, the issue is a quite similar type of nonlinearity in the dispersion
tensor. We are currently working to develop techniques tailored to prove
uniqueness of solutions to this class of multiscale systems.

## Acknowledgments

We thank H. Notsu (Kanazawa) for fruitful discussions during his visit to
Karlstad. The work of V.R., S.N., and A.M. is partially supported by the
Swedish Research Council’s project “Homogenization and dimension reduction of
thin heterogeneous layers” (grant nr. VR 2018-03648). The research activity of
M.E. is funded by the European Union’s Horizon 2022 research and innovation
program under the Marie Skłodowska-Curie fellowship project MATT (project nr.
101061956). R.L. and A.M. are grateful to Carl Tryggers Stiftelse for their
financial support via the grant CTS 21:1656.

## References

  * [1] A. Abdulle, E Weinan, B. Engquist, and E. Vanden-Eijnden, _The heterogeneous multiscale method_ , Acta Numerica 21 (2012), 1–87.
  * [2] G. Allaire, _Homogenization and two scale convergence_ , SIAM J. Math. Anal. 23 (1992), no. 6, 1482–1518.
  * [3] G. Allaire and R. Brizzi, _A multiscale finite element method for numerical homogenization_ , Multiscale Modeling & Simulation 4 (2005), no. 3, 790–812.
  * [4] G. Allaire, R. Brizzi, A. Mikelić, and A. Piatnitski, _Two-scale expansion with drift approach to the Taylor dispersion for reactive transport through porous media_ , Chemical Engineering Science 65 (2010), no. 7, 2292–2300.
  * [5] G. Allaire and H. Hutridurga, _Upscaling nonlinear adsorption in periodic porous media – homogenization approach_ , Applicable Analysis 95 (2016), no. 10, 2126–2161.
  * [6] G. Allaire and A. Raphael, _Homogenization of a convection–diffusion model with reaction in a porous medium_ , Comptes Rendus Mathematique 344 (2007), no. 8, 523–528.
  * [7] B. Amaziane, A. Bourgeat, and M. Jurak, _Effective macrodiffusion in solute transport through heterogeneous porous media_ , Multiscale Modeling & Simulation 5 (2006), no. 1, 184–204.
  * [8] R. Aris, _On the dispersion of a solute in a fluid flowing through a tube_ , Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 235 (1956), no. 1200, 67–77.
  * [9] D. N. Arnold, F. Brezzi, and M. Fortin, _A stable finite element for the Stokes equations_ , Calcolo 21 (1984), no. 4, 337–344.
  * [10] J.-L. Auriault, C. Moyne, and H. P. Amaral Souto, _On the asymmetry of the dispersion tensor in porous media_ , Transport in Porous Media 85 (2010), no. 3, 771–783.
  * [11] J. Bear, _Dynamics of Fluids in Porous Media_ , Dover Publications, 1988\.
  * [12] J. Bear and L. G. Fel, _A phenomenological approach to modeling transport in porous media_ , Transport in Porous Media 92 (2012), 649–665.
  * [13] K. Bhattacharya, G. Friesecke, and R. D. James, _The mathematics of microstructure and the design of new materials_ , Proceedings of the National Academy of Sciences 96 (1999), no. 15, 8332–8333.
  * [14] C. Bringedal, L. von Wolff, and I. S. Pop, _Phase field modeling of precipitation and dissolution processes in porous media: Upscaling and numerical experiments_ , Multiscale Modeling & Simulation 18 (2020), no. 2, 1076–1112.
  * [15] C. Le Bris, F. Legoll, and F. Madiot, _Multiscale finite element methods for advection-dominated problems in perforated domains_ , Multiscale Modeling & Simulation 17 (2019), no. 2, 773–825.
  * [16] E. N. M. Cirillo, I. de Bonis, A. Muntean, and O. Richardson, _Upscaling the interplay between diffusion and polynomial drifts through a composite thin strip with periodic microstructure_ , Meccanica 55 (2020), 2159–2178.
  * [17] E. N. M. Cirillo, O. Krehel, A. Muntean, R. van Santen, and A. Sengar, _Residence time estimates for asymmetric simple exclusion dynamics on strips_ , Physica A: Statistical Mechanics and its Applications 442 (2016), 436–457.
  * [18] P. Donovan, Y. Chehreghanianzabi, M. Rathinam, and S. P. Zustiak, _Homogenization theory for the prediction of obstructed solute diffusivity in macromolecular solutions_ , PLOS ONE 11 (2016), no. 1, 1–26.
  * [19] M. Eden, C. Nikolopoulos, and A. Muntean, _A multiscale quasilinear system for colloids deposition in porous media: Weak solvability and numerical simulation of a near-clogging scenario_ , Nonlinear Analysis: Real World Applications 63 (2022), 103408.
  * [20] L. C. Evans, _Partial Differential Equations_ , vol. 19, American Mathematical Society, 2010.
  * [21] L. Formaggia, J.-F. Gerbeau, F. Nobile, and A. Quarteroni, _Numerical treatment of defective boundary conditions for the Navier–Stokes equations_ , SIAM Journal on Numerical Analysis 40 (2002), no. 1, 376–401.
  * [22] T. Goudon and F. Poupaud, _Homogenization of transport equations: A simple PDE approach to the Kubo formula_ , Bulletin des Sciences Mathématiques 131 (2007), no. 1, 72–88.
  * [23] P. Henning, M. Ohlberger, and B. Schweizer, _An adaptive multiscale finite element method_ , Multiscale Modeling & Simulation 12 (2014), no. 3, 1078–1107.
  * [24] E. Ijioma and A. Muntean, _Fast drift effects in the averaging of a filtration combustion system: A periodic homogenization approach_ , Quarterly of Applied Mathematics 77 (2019), no. 1, 71–104.
  * [25] H. P. Langtangen and A. Logg, _Solving PDEs in Python_ , Springer International Publishing, 2016.
  * [26] W. Li, L. Jiang, W. Jiang, Y. Wu, X. Guo, Z. Li, H. Yuan, and M. Luo, _Recent advances of boron nitride nanosheets in hydrogen storage application_ , Journal of Materials Research and Technology 26 (2023), 2028–2042.
  * [27] M. Lind, A. Muntean, and O. Richardson, _A semidiscrete Galerkin scheme for a coupled two-scale elliptic–parabolic system: well-posedness and convergence approximation rates_ , BIT Numerical Mathematics 60 (2020), no. 4, 999–1031.
  * [28] J-L. Lions, _Quelques Méthodes de Résolution des Problèmes aux Limites Non-Linéaires_ , Dunod, 1969.
  * [29] A. Logg, K.-A. Mardal, and G. Wells, _Automated Solution of Differential Equations by the Finite Element Method: The FEniCS Book_ , vol. 84, Springer Science & Business Media, 2012.
  * [30] D. Lukkassen, G. Nguetseng, and P. Wall, _Two-scale convergence_ , International journal of pure and applied mathematics 2 (2002), no. 1, 35–86.
  * [31] E. Marušić-Paloka and A. L. Piatnitski, _Homogenization of a nonlinear convection-diffusion equation with rapidly oscillating coefficients and strong convection_ , Journal of the London Mathematical Society 72 (2005), 391–409.
  * [32] D.W. McLaughlin, G.C. Papanicolaou, and O.R. Pironneau, _Convection of microstructure and related problems_ , SIAM Journal on Applied Mathematics 45 (1985), no. 5, 780–797.
  * [33] A. Muntean and C. Nikolopoulos, _Colloidal transport in locally periodic evolving porous media—an upscaling exercise_ , SIAM Journal on Applied Mathematics 80 (2020), no. 1, 448–475.
  * [34] S. Nepal, Y. Wondmagegne, and A. Muntean, _Analysis of a fully discrete approximation to a moving-boundary problem describing rubber exposed to diffusants_ , Applied Mathematics and Computation 442 (2023), 127733.
  * [35] C. Nikolopoulos, M. Eden, and A. Muntean, _Multiscale simulation of colloids ingressing porous layers with evolving internal structure: A computational study_ , GEM-International Journal on Geomathematics 14 (2023), no. 1, 1.
  * [36] M. B. Olivares, C. Bringedal, and I. S. Pop, _A two-scale iterative scheme for a phase-field model for precipitation and dissolution in porous media_ , Applied Mathematics and Computation 396 (2021), 125933.
  * [37] R. Precup, _Linear and Semilinear Partial Differential Equations: An Introduction_ , Berlin, Boston: De Gruyter, 2013.
  * [38] V. Raveendran, E.N.M. Cirillo, I. de Bonis, and A. Muntean, _Scaling effects on the periodic homogenization of a reaction-diffusion-convection problem posed in homogeneous domains connected by a thin composite layer_ , Quarterly of Applied Mathematics 80 (2022), 157–200.
  * [39] V. Raveendran, E.N.M Cirillo, and A. Muntean, _Upscaling of a reaction-diffusion-convection problem with exploding non-linear drift_ , Quarterly of Applied Mathematics 80 (2022), 641–667.
  * [40] V. Raveendran, I. de Bonis, E. N.M. Cirillo, and A. Muntean, _Homogenization of a reaction-diffusion problem with large nonlinear drift and Robin boundary data_ , arXiv preprint arXiv:2307.04567 (2023), 1–37.
  * [41] R.J. Schotting, H. Moser, and S.M. Hassanizadeh, _High-concentration-gradient dispersion in porous media: experiments, analysis and approximations_ , Advances in Water Resources 22 (1999), no. 7, 665–680.
  * [42] R. E. Showalter, _Distributed microstructure models of porous media_ , Flow in Porous Media: Proceedings of the Oberwolfach Conference, June 21–27, 1992, Springer, 1993, pp. 155–163.
  * [43] G. I. Taylor, _Conditions under which dispersion of a solute in a stream of solvent can be used to measure molecular diffusion_ , Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 225 (1954), no. 1163, 473–477.
  * [44] D. Wiedemann and M. A. Peter, _Homogenisation of local colloid evolution induced by reaction and diffusion_ , Nonlinear Analysis 227 (2023), 113168\.
  * [45] B. D. Wood, F. Cherblanc, M. Quintard, and S. Whitaker, _Volume averaging for determining the effective dispersion tensor: Closure using periodic unit cells and comparison with ensemble averaging_ , Water Resources Research 39 (2003), no. 8.