# A domain decomposition solution of the Stokes-Darcy system in 3D based on
boundary integrals

Svetlana Tlupova Department of Mathematics, Farmingdale State College, SUNY,
Farmingdale, NY 11735, USA<EMAIL_ADDRESS>
###### Abstract

A framework is developed for a robust and highly accurate numerical solution
of the coupled Stokes-Darcy system in three dimensions. The domain
decomposition method is based on a Dirichlet-Neumann type splitting of the
interface conditions and solving separate Stokes and Darcy problems
iteratively. Second kind boundary integral equations are formulated for each
problem. The integral equations use a smoothing of the kernels that achieves
high accuracy on the boundary, and a straightforward quadrature to discretize
the integrals. Numerical results demonstrate the convergence, accuracy, and
dependence on parameter values of the iterative solution for a problem of
viscous flow around a porous sphere with a known analytical solution, as well
as more general surfaces.

Keywords: Stokes-Darcy, domain decomposition methods, boundary integral
equations, porous sphere.

##  1 Introduction

Many important environmental and biological phenomena that occur in our daily
lives involve fluids partly flowing freely and partly filtrating through a
porous medium. Such processes are of interest in physiology when studying
cancer growth and filtration of blood through arterial vessel walls. The
free/porous flow plays an important role in numerous industrial applications
of air and oil filter design, oil exploration, or chemical reactor
simulations. There is also considerable environmental interest in studies of
groundwater remediation, geologic CO2 sequestration, and bacterial biofilms.

Due to its wide applicability, the modeling of free/porous flow has recently
received considerable attention, and mathematical and numerical analyses have
been done. The flow is often modeled by the incompressible Stokes equations in
the free domain and the Darcy equations in the porous domain, coupled across
the interface through suitable conditions. The numerical solutions of these
flows have been almost exclusively based on the finite element method [53, 38,
20, 15, 27, 14, 11, 17, 2, 21, 22, 24, 26, 28, 37, 41, 42]; a singularity
method was used to obtain the solution in [23, 48], and a MAC scheme in [52].
There are also recently developed partitioned methods for time-dependent
Stokes-Biot problems, which utilize Robin-type coupling conditions at the
interface [10, 1].

Despite the significant efforts in recent years to design numerical methods of
coupled free/porous flows, these problems remain a challenge today, owing to
the multiphysics nature of the system and the large differences in the
physical parameters in the governing equations. The severe limitations of the
direct solution necessitate the development of techniques that would enable
one to split the system and solve the smaller Stokes and Darcy problems
separately while exchanging information through the boundary conditions only.
These techniques, called domain decomposition methods (DDM), have been used to
solve various problems [47]. DDMs are based on dividing the computational
domain into several subdomains with or without overlap, and formulating the
so-called transmission conditions along the subdomain interfaces. The coupling
across the interface is then replaced by an iterative procedure, where the
smaller problems are solved independently at each iteration. For the coupled
flow problem, similar techniques have been developed and investigated in the
context of the finite element method, e.g. [53, 12, 57, 47, 21, 22, 24]. In
particular, we are interested in the Dirichlet-Neumann iterative methods [47,
21], where the interface conditions are split between the Stokes and Darcy
problems.

In this paper we develop a simple sequential iterative algorithm inspired by
the Dirichlet-Neumann methods to decouple the free and porous problems, and
formulate highly accurate second-kind boundary integral equations to solve for
the fluid quantities. Boundary integral equation methods (BIEM) are extremely
powerful in solving many differential equations [3, 16, 35, 30], and their
importance is well recognized. They have been extensively applied to solve
Stokes problems [45, 18, 19, 30, 25], as well as porous media flows [39, 49,
40, 46]. In two dimensions, a Robin-Robin DDM was designed for the Stokes-
Darcy system and boundary integral equations were used to solve the local
problems [9]. Non-local operators based on the integral formulations were used
in the Robin transmission conditions, leading to a significant improvement in
convergence of the DDM. This method is more efficient in overcoming the
challenges of the direct solution of [55, 8]. However, to the best of our
knowledge, no constructive and efficient way of solving the coupled
free/porous system in 3D has been developed using BIEM. Attempting a direct
solution of the coupled Stokes-Darcy system with BIEM will lead to a
formulation that is inevitably of mixed-kind, due to the incompatibility of
the Stokes and Darcy differential operators. Once the system is split through
the DDM procedure however, suitable formulations can be applied to each
problem. The significant advantages of BIEM are the reduced dimensionality as
the boundary alone needs to be discretized, high achievable accuracy at points
on and off the boundary, the availability of well-conditioned formulations
where the efficiency of the solution is maintained with grid refinement, and
efficient solutions for problems where the free domain is large or infinite,
with boundaries that have complicated geometry, as well as moving boundaries.

To address the singularities that develop in the integrands when solving the
integral equations, we use a regularization technique developed in [4, 5] for
Laplace’s equation, and recently developed for the Stokes equations in [54].
The smoothing can be chosen to have high order on the boundary, so that the
accuracy in solving the integral equations generally approaches $O(h^{5})$,
where $h$ is the grid size in the coordinate planes. This approach is quite
simple in that it does not require any special quadrature for evaluating the
integrals on the boundary. To discretize the integrals a quadrature rule for
closed surfaces from [59, 6] is applied. This quadrature rule works well for
general surfaces without requiring coordinate charts, using projections on
coordinate planes instead. As the weights for the quadrature points are found
from a partition of unity on the unit sphere, they do not depend on the
particular surface. With this approach, the overall formulation is simple to
implement, and the data structure needed to describe the boundary is minimal.

The paper is organized as follows. In Section 2, we state the governing
equations and suitable interface conditions for a Stokes-Darcy system. In
Section 3, we outline a Dirichlet-Neumann type iterative procedure to split
and sequentially solve the Darcy and Stokes problems. A convergence analysis
of this algorithm for a spherical geometry based on spherical harmonics
representations is conducted in Section 4. In Section 5, a boundary integral
formulation is developed for both Darcy and Stokes equations, including a
description of the treatment of singularities in the integrands in Section 5.3
and the quadrature in Section 5.4. The complete algorithm is presented in
Section 6 and the discrete iteration operator is computed in Section 6.1.
Numerical results are presented in Section 7, where we use a known analytical
solution for a viscous flow around a porous sphere as a benchmark problem, as
well as test the algorithm on other geometries.

##  2 Problem formulation

To describe the model equations, we denote the free flow quantities and the
porous domain quantities by subscripts $S$ and $D$, respectively. Figure 1
shows an example two-dimensional setting, where a free fluid flow in a bounded
region $\Omega_{S}$ and a flow in a porous region $\Omega_{D}$ are coupled
across the common interface $\Sigma$.

Figure 1: The domain with the free fluid region $\Omega_{S}$ and the porous
medium region $\Omega_{D}$, with the common interface $\Sigma$.

In the region of free fluid flow $\Omega_{S}$, we assume the low Reynolds
number regime modeled by the incompressible Stokes equations,

$\textrm{In}\ \Omega_{S}:\left\\{\begin{array}[]{l}\nabla
p_{S}-\mu\Delta\mathbf{u}_{S}=0,\\\\[2.0pt]
\nabla\cdot\mathbf{u}_{S}=0,\end{array}\right.$ (1)

where $\mu$ is the fluid viscosity, $p_{S}$ and $\mathbf{u}_{S}$ are the fluid
pressure and velocity, respectively. Darcy’s equations govern the flow in the
porous medium,

$\textrm{In}\ \Omega_{D}:\left\\{\begin{array}[]{l}\nabla p_{D}+\mu
K^{-1}\mathbf{u}_{D}=0,\\\\[2.0pt]
\nabla\cdot\mathbf{u}_{D}=0.\end{array}\right.$ (2)

Here $\mathbf{u}_{D}$ and $p_{D}$ are respectively the (averaged) fluid
velocity and the hydrostatic pressure, and $K$ is a symmetric and positive
definite permeability tensor that depends on the microstructure of the porous
medium. The medium is called homogeneous if the permeability $K(\mathbf{x})$
is constant throughout the domain, and isotropic when permeability is a
scalar, $K=\kappa I$. In this work, we treat porous media that are homogeneous
and isotropic.

The interface between these two systems of partial differential equations
serves as a replacement of the boundary layer through which the velocity
changes rapidly. A lot of attention has been devoted to the formulation and
analysis of appropriate coupling conditions [7, 50, 31, 32, 43, 44, 51]. The
incompressibility condition leads to continuity of normal velocity components,

$\mathbf{u}_{S}\cdot\mathbf{n}_{S}=-\mathbf{u}_{D}\cdot\mathbf{n}_{D},$ (3)

where $\mathbf{n}_{S}$ and $\mathbf{n}_{D}$ are the unit normal vectors that
point out of the regions $\Omega_{D}$ and $\Omega_{S}$ respectively, so that
$\mathbf{n}_{D}=-\mathbf{n}_{S}$ on $\Sigma$ (see Fig. 1). In addition, a
suitable condition on the tangential velocity is the slip condition of
Beavers-Joseph-Saffman [7, 50],

$\bm{\tau}_{S}\cdot\bm{\sigma}_{S}\cdot\mathbf{n}_{S}=\frac{\gamma}{\sqrt{\kappa}}\mathbf{u}_{S}\cdot\bm{\tau}_{S},$
(4)

where
$\bm{\sigma}_{S}=-p_{S}I+\mu[\nabla\mathbf{u}_{S}+(\nabla\mathbf{u}_{S})^{T}]$
is the Stokes stress tensor, $\bm{\tau}_{S}$ is the tangent vector to the
interface, and $\gamma$ is a dimensionless quantity that depends on the
structure of the porous material. Well-posedness of the mathematical
formulation with this condition has been analyzed for steady flows [11, 21,
24, 37] as well as unsteady flows [13]. The final condition represents the
balance of normal forces across the interface,

$\mathbf{n}_{S}\cdot\bm{\sigma}_{S}\cdot\mathbf{n}_{S}=-p_{D},$ (5)

which allows the pressure to be discontinuous across the interface.

A fundamental difficulty in the development of robust computational methods
for the free/porous flow system (1)-(5) is due to the fact that the PDEs that
govern the fluid flow in the free region and the flow in the porous medium
have different properties. The differential operators in the two subdomains
are incompatible, inhibiting the simultaneous solution of the equations by a
direct method. In addition, the equations include parameters, the viscosity
$\mu$ and permeability $\kappa$, that could differ by several orders of
magnitude in practice.

##  3 The Dirichlet-Neumann iterative method

The severe limitations of the direct solution necessitate the development of
techniques that would enable one to split the system and solve the smaller
Stokes and Darcy problems separately while exchanging information through the
boundary conditions only. We develop a numerical scheme based on the
sequential Dirichlet-Neumann (D-N) iterative method [47, 21]. The outline of
the algorithm is as follows. First, we define
$q=-\mathbf{u}_{S}\cdot\mathbf{n}_{S}$ as the iteration variable on the
interface $\Sigma$, and choose an initial guess $q^{(0)}=0$. Then, for
$k=1,2,3,...$ until convergence, perform the following steps:

  1. 1.

Solve the Darcy problem with the interface condition (3),

$\left\\{\begin{array}[]{ll}\nabla
p^{(k)}_{D}+\mu\kappa^{-1}\mathbf{u}^{(k)}_{D}=0&\textrm{in}\
\Omega_{D},\\\\[2.0pt] \nabla\cdot\mathbf{u}^{(k)}_{D}=0&\textrm{in}\
\Omega_{D},\\\\[2.0pt]
\mathbf{u}^{(k)}_{D}\cdot\mathbf{n}_{D}=q^{(k-1)}&\textrm{on}\
\Sigma.\end{array}\right.$ (6)

The interface condition is a Neumann condition for $p_{D}$.

  2. 2.

Using $p^{(k)}_{D}$ from Step 1, solve the Stokes problem with the interface
conditions (4)-(5),

$\left\\{\begin{array}[]{ll}\nabla
p^{(k)}_{S}-\mu\Delta\mathbf{u}^{(k)}_{S}=\mathbf{0}&\textrm{in}\
\Omega_{S},\\\\[2.0pt] \nabla\cdot\mathbf{u}^{(k)}_{S}=0&\textrm{in}\
\Omega_{S},\\\\[2.0pt]
\bm{\tau}_{S}\cdot\bm{\sigma}^{(k)}_{S}\cdot\mathbf{n}_{S}=\frac{\gamma}{\sqrt{\kappa}}\mathbf{u}^{(k-1)}_{S}\cdot\bm{\tau}_{S}&\textrm{on}\
\Sigma,\\\\[2.0pt]
\mathbf{n}_{S}\cdot\bm{\sigma}^{(k)}_{S}\cdot\mathbf{n}_{S}=-p^{(k)}_{D}&\textrm{on}\
\Sigma.\end{array}\right.$ (7)

  3. 3.

Update the iteration variable

$q^{(k)}=(1-\theta)\,q^{(k-1)}-\theta\,\mathbf{u}^{(k)}_{S}\cdot\mathbf{n}_{S}$
(8)

using $\mathbf{u}^{(k)}_{S}$ from Step 2, where $\theta\in(0,1)$ is a
relaxation parameter.

The convergence criterion is taken as the relative error between $q^{(k+1)}$
and $q^{(k)}$, which we call the DDM residual, falling below a prescribed
tolerance. The solutions of the now separate Darcy (6) and Stokes (7) problems
are based on the boundary integral formulations which we describe in Section
5.

##  4 Convergence analysis

We analyze the convergence properties of the Dirichlet-Neumann iterative
scheme using the spherical harmonics representation of both problems, assuming
as the interface a sphere of radius $R$ centered at the origin, with the Darcy
problem inside and the Stokes problem outside the sphere.

###  4.1 Darcy problem inside a sphere

We consider the problem

$\left\\{\begin{array}[]{ll}\tilde{\kappa}\,\nabla p+\mathbf{u}=0&\textrm{in}\
\Omega,\\\\[2.0pt] \nabla\cdot\mathbf{u}=0&\textrm{in}\ \Omega,\\\\[2.0pt]
\mathbf{u}\cdot\mathbf{n}=q&\textrm{on}\ \Sigma,\end{array}\right.$ (9)

where $\Omega$ is the sphere of radius $R$ centered at the origin, $\Sigma$ is
its boundary, $\mathbf{n}$ is the outward unit normal vector, and
$\tilde{\kappa}=\kappa/\mu$.

In spherical coordinates $(r,\theta,\phi)$, the general solution to Laplace’s
equation $\Delta p=0$ in $\Omega$ can be written as

$p=\sum_{n=0}^{\infty}p_{n},$ (10)

with

$p_{n}=\sum_{m=-n}^{n}a_{n}^{m}\,r^{n}\,Y_{n}^{m}(\theta,\phi),$ (11)

where $p_{n}$ is a solid harmonic of order $n$, $Y_{n}^{m}$ is the spherical
harmonic function of degree $n$ and order $m$, and $a_{n}^{m}\in\mathbb{C}$
are constants.

With $\mathbf{n}=\mathbf{x}/r$ and

$\nabla p=\sum_{n=0}^{\infty}\sum_{m=-n}^{n}\left[\frac{\partial}{\partial
r}(a_{n}^{m}\,r^{n})\,\frac{\mathbf{x}}{r}\,Y_{n}^{m}(\theta,\phi)+a_{n}^{m}\,r^{n}\,\nabla
Y_{n}^{m}(\theta,\phi)\right],$ (12)

for the normal velocity we get

$\mathbf{u}\cdot\mathbf{n}=-\tilde{\kappa}\,\nabla
p\cdot\mathbf{n}=-\tilde{\kappa}\sum_{n=0}^{\infty}\sum_{m=-n}^{n}\left[n\,a_{n}^{m}\,r^{n-1}\,Y_{n}^{m}(\theta,\phi)\right].$
(13)

Using a similar representation for $q$ on $\Sigma$,

$q=\sum_{n=0}^{\infty}\sum_{m=-n}^{n}q_{n}^{m}\,Y_{n}^{m}(\theta,\phi),$ (14)

and imposing the Neumann boundary condition in (9) on $r=R$, we find for the
coefficients $a_{n}^{m}$,

$-\tilde{\kappa}\,n\,a_{n}^{m}\,R^{n-1}=q_{n}^{m}.$ (15)

###  4.2 Stokes problem outside a sphere

We now consider the solution of the Stokes problem with a no-slip boundary
condition for simplicity,

$\left\\{\begin{array}[]{ll}\nabla p-\mu\,\Delta\mathbf{u}=0&\textrm{in}\
\Omega^{c},\\\\[2.0pt] \nabla\cdot\mathbf{u}=0&\textrm{in}\
\Omega^{c},\\\\[2.0pt] \mathbf{u}\cdot\bm{\tau}=0&\textrm{on}\
\Sigma,\\\\[2.0pt]
\mathbf{n}\cdot\bm{\sigma}\cdot\mathbf{n}=-p_{D}&\textrm{on}\
\Sigma,\end{array}\right.$ (16)

where $p_{D}$ is the Darcy pressure on $\Sigma$ and $\Omega^{c}$ is the
exterior of the sphere of radius $R$.

The solution to the Stokes equations in spherical harmonics is referred to as
Lamb’s general solution [36]. Since the Stokes pressure satisfies Laplace’s
equation, it can be represented similar to (10)-(11), but for the exterior
solution we keep only the negative harmonics,

$p=\sum_{n=0}^{\infty}p_{-n-1},$ (17)

with

$p_{-n-1}=\sum_{m=-n}^{n}b_{n}^{m}\,r^{-n-1}\,Y_{n}^{m}(\theta,\phi).$ (18)

The disturbance velocity field is represented by (see [34])

$\mathbf{u}=\sum_{n=1}^{\infty}\left[-\frac{(n-2)\,r^{2}\,\nabla
p_{-n-1}}{2\mu\,n(2n-1)}+\frac{(n+1)\,\mathbf{x}\,p_{-n-1}}{\mu\,n(2n-1)}+\nabla\phi_{-n-1}+\nabla\times(\mathbf{x}\,\chi_{-n-1})\right],$
(19)

where

$\phi_{-n-1}=\sum_{m=-n}^{n}c_{n}^{m}\,r^{-n-1}\,Y_{n}^{m}(\theta,\phi),$ (20)
$\chi_{-n-1}=\sum_{m=-n}^{n}d_{n}^{m}\,r^{-n-1}\,Y_{n}^{m}(\theta,\phi).$ (21)

The first two terms in (19) correspond to a particular solution, while the
rest corresponds to the homogeneous solution, constructed from a potential
(the $\nabla\phi$ term) and a toroidal field (the
$\nabla\times(\mathbf{x}\,\chi)$ term) [34].

The stress vector acting across the sphere can be expressed in general as (see
[29] equation (3-2.37) or [34] exercise 4.6, p.104),

$\displaystyle\bm{\sigma}_{n}=\bm{\sigma}\cdot\mathbf{n}=\frac{1}{r}\sum_{n=-\infty}^{\infty}\left[\frac{n(n+2)\,r^{2}\,\nabla
p_{n}}{(n+1)(2n+3)}-\frac{(2n^{2}+4n+3)\,\mathbf{x}\,p_{n}}{(n+1)(2n+3)}\right]$
$\displaystyle+\frac{\mu}{r}\sum_{n=-\infty}^{\infty}\left[2(n-1)\nabla\phi_{n}+(n-1)\nabla\times(\mathbf{x}\,\chi_{n})\right].$
(22)

For the exterior problem, we again only keep the negative harmonics and
replace $n$ in (4.2) with $-n-1$, to obtain

$\displaystyle\bm{\sigma}_{n}=\bm{\sigma}\cdot\mathbf{n}=\frac{1}{r}\sum_{n=1}^{\infty}\left[\frac{(n-1)(n+1)\,r^{2}\,\nabla
p_{-n-1}}{n(2n-1)}-\frac{(2n^{2}+1)\,\mathbf{x}\,p_{-n-1}}{n(2n-1)}\right]$
$\displaystyle-\frac{\mu}{r}\,\sum_{n=1}^{\infty}\left[2(n+2)\nabla\phi_{-n-1}+(n+2)\nabla\times(\mathbf{x}\,\chi_{-n-1})\right].$
(23)

We now turn to the boundary conditions. The no-slip condition on $r=R$ implies
(see [34] p.89) that $\chi$ is not needed and

$\sum_{n=1}^{\infty}\left[-\frac{n(n+1)\,R}{2\mu\,(2n-1)}\,p_{-n-1}|_{r=R}+\frac{(n+1)(n+2)}{R}\phi_{-n-1}|_{r=R}\right]=0,$
(24)

which leads to

$\phi_{-n-1}|_{r=R}=\frac{n\,R^{2}}{2\mu\,(n+2)(2n-1)}p_{-n-1}|_{r=R}.$ (25)

Note that for the radial component of velocity we have from (19),

$u_{r}=\sum_{n=1}^{\infty}\left[-\frac{(n-2)\,r^{2}}{2\mu\,n(2n-1)}\,\frac{\partial
p_{-n-1}}{\partial
r}+\frac{(n+1)\,r\,p_{-n-1}}{\mu\,n(2n-1)}+\frac{\partial\phi_{-n-1}}{\partial
r}\right].$ (26)

Using Euler’s theorem for homogeneous polynomials,

$r\,\frac{\partial h_{n}}{\partial r}=n\,h_{n},$ (27)

where $h_{n}$ is any solid harmonic of order $n$, the radial velocity in (26)
simplifies to

$u_{r}=\sum_{n=1}^{\infty}\left[\frac{(n+1)\,r\,p_{-n-1}}{2\mu\,(2n-1)}-\frac{(n+1)}{r}\phi_{-n-1}\right].$
(28)

Similarly, for the radial component of the stress vector $\bm{\sigma}_{n}$,

$\displaystyle\sigma_{nr}=\frac{1}{r}\sum_{n=1}^{\infty}\left[\frac{(n-1)(n+1)\,r^{2}}{n(2n-1)}\,\frac{\partial
p_{-n-1}}{\partial r}-\frac{(2n^{2}+1)\,r\,p_{-n-1}}{n(2n-1)}\right]$
$\displaystyle-\frac{\mu}{r}\,\sum_{n=1}^{\infty}\left[2(n+2)\frac{\partial\phi_{-n-1}}{\partial
r}\right],$ (29)

which simplifies again using (27),

$\sigma_{nr}=\sum_{n=1}^{\infty}\left[-\frac{(n^{2}+3n-1)}{(2n-1)}p_{-n-1}+\frac{2\mu}{r^{2}}(n+1)(n+2)\phi_{-n-1}\right].$
(30)

On the boundary $r=R$ using (25), this becomes

$\sigma_{nr}|_{r=R}=\sum_{n=1}^{\infty}\left[-p_{-n-1}\right]_{r=R},$ (31)

and therefore enforcing the boundary condition on the normal stress

$\bm{\sigma}_{n}\cdot\mathbf{n}=\sigma_{nr}|_{r=R}=-p_{D}$ (32)

leads to

$\sum_{n=1}^{\infty}\sum_{m=-n}^{n}b_{n}^{m}\,R^{-n-1}\,Y_{n}^{m}(\theta,\phi)=\sum_{n=1}^{\infty}\sum_{m=-n}^{n}a_{n}^{m}\,R^{n}\,Y_{n}^{m}(\theta,\phi).$
(33)

Using (15), we write

$b_{n}^{m}=-\frac{R^{n+2}}{\tilde{\kappa}\,n}q_{n}^{m}.$ (34)

###  4.3 Spectrum of the Dirichlet-Neumann procedure

We now diagonalize the iteration operator of the Dirichlet-Neumann method
defined on the sphere. To derive the spectrum of the iteration operator, note
that the iteration variable $q$ is updated at each iteration according to (8),
so we need the normal component of Stokes velocity on the sphere. We evaluate
the radial velocity (28) on the boundary $r=R$ and use (25) to get for the
normal Stokes velocity,

$u_{r}|_{r=R}=\sum_{n=1}^{\infty}\left[\frac{(n+1)\,R}{\mu\,(n+2)(2n-1)}p_{-n-1}|_{r=R}\right].$
(35)

Using (34), this becomes

$u_{r}|_{r=R}=-\frac{R^{2}}{\kappa}\sum_{n=1}^{\infty}\left[\frac{(n+1)}{n(n+2)(2n-1)}\sum_{m=-n}^{n}q_{n}^{m}\,Y_{n}^{m}(\theta,\phi)\right].$
(36)

For the DDM update we get

$q^{(k)}=(1-\theta)q^{(k-1)}+\theta u_{r}^{(k)}|_{r=R},$ (37)

which we write in the form of the diagonalized iteration operator
$q^{(k)}=\mathcal{A}q^{(k-1)}$ as

$q_{n}^{m(k)}=\mathcal{A}_{n}q_{n}^{m(k-1)},$ (38)

where

$\mathcal{A}_{n}=(1-\theta)-\theta\frac{R^{2}}{\kappa}\frac{(n+1)}{n(n+2)(2n-1)},\qquad
n=1,2,3,...$ (39)

Figure 2 shows these coefficients for $n=1,...,50$ for different values of
permeability $\kappa$ as well as different $\theta$. We can observe that for
$\kappa=1$, a larger value of $\theta$ leads to smaller coefficients
$\mathcal{A}_{n}$ and therefore the iterations are expected to converge
faster, as we will show in the numerical tests. For smaller values of $\kappa$
however, the coefficients $\mathcal{A}_{n}$ get close to 1 as $n$ increases,
independent of the value of $\theta$. This can hinder the convergence rate of
the algorithm when using a successive approximation technique as in (8). We
will show in our numerical tests that this indeed is what happens, and switch
to using a more efficient method, and one that does not require all
eigenvalues to be less than 1, such as GMRES to help mitigate this issue.
Interesting approaches that will be investigated in future work are (i)
different decoupling, e.g. Robin-Robin transmission conditions, and (ii) using
in place of $\theta$ non-local operators based on the boundary integral
formulation to speed up the convergence rates. Similar techniques were shown
to dramatically speed up the iterations in two dimensions [9].

Figure 2: Harmonic coefficients $\mathcal{A}_{n}$ from (39) of the iteration
operator on the unit sphere, for permeability values
$\kappa=1,10^{-2},10^{-4}$.

##  5 Boundary integral formulations

There are two main approaches to formulating boundary integral equations,
usually referred to as the ’direct’ and ’indirect’ approaches. The direct
approach is based on the Green’s representation formula and is the approach we
take in this work. While the analytic properties of the integral operators are
the same with both approaches, our primary interest in the direct approach is
that the resulting equations are for the physically relevant Dirichlet or
Neumann data.

###  5.1 Darcy problem

The Green’s representation formula,

$p_{D}(\mathbf{y})=-\int_{\partial\Omega_{D}}\frac{\partial
p_{D}(\mathbf{x})}{\partial
n(\mathbf{x})}\,G(\mathbf{y,x})\,dS(\mathbf{x})+\int_{\partial\Omega_{D}}p_{D}(\mathbf{x})\,\frac{\partial
G(\mathbf{y,x})}{\partial
n(\mathbf{x})}\,dS(\mathbf{x}),\qquad\mathbf{y}\in\Omega_{D},$ (40)

represents any harmonic function $p_{D}$ in terms of its Cauchy data, that is,
its boundary values and its normal derivative on the boundary. The first
integral in (40) is the single-layer potential and the second integral is the
double-layer potential, where the kernels are
$G(\mathbf{y,x})=-1/4\pi|\mathbf{y}-\mathbf{x}|$ and its derivative in the
outward normal $\frac{\partial G(\mathbf{y,x})}{\partial
n(\mathbf{x})}=\nabla_{\mathbf{x}}G(\mathbf{y},\mathbf{x})\cdot\mathbf{n}(\mathbf{x})=-(\mathbf{y}-\mathbf{x})\cdot\mathbf{n}(\mathbf{x})/4\pi|\mathbf{y}-\mathbf{x}|^{3}$.
The double layer potential undergoes a discontinuity across the boundary, and
the jump relations

$\int_{\partial\Omega_{D}}\frac{\partial G(\mathbf{y,x})}{\partial
n(\mathbf{x})}\,dS(\mathbf{x})=\left\\{\begin{array}[]{ll}1,&\mathbf{y}\in\Omega_{D}\\\
1/2,&\mathbf{y}\in\partial\Omega_{D}\\\
0,&\mathbf{y}\in\mathbb{R}^{3}\setminus\bar{\Omega}_{D}\end{array}\right.$
(41)

are used to extend (40) to the boundary,

$\frac{1}{2}p_{D}(\mathbf{y})=-\int_{\partial\Omega_{D}}\frac{\partial
p_{D}(\mathbf{x})}{\partial
n(\mathbf{x})}\,G(\mathbf{y,x})\,dS(\mathbf{x})+\int_{\partial\Omega_{D}}p_{D}(\mathbf{x})\,\frac{\partial
G(\mathbf{y,x})}{\partial
n(\mathbf{x})}\,dS(\mathbf{x}),\qquad\mathbf{y}\in\partial\Omega_{D}.$ (42)

We write (42) as a second-kind integral equation assuming known Neumann data
$\partial p_{D}/\partial n$ on the boundary,

$\frac{1}{2}p_{D}=Hp_{D}+b_{D},\qquad\qquad b_{D}=-L\frac{\partial
p_{D}}{\partial n},$ (43)

where by $L$ and $H$ we denote the integral operators corresponding to the
single- and double-layer potentials, respectively. The integral equation (43)
can be solved by successive approximations, otherwise known as a Neumann’s
iteration scheme (see, for example, [35] or [30], Sec. 5.6.7),

$p_{D}^{(n+1)}=\frac{1}{2}p_{D}^{(n)}+Hp_{D}^{(n)}+b_{D},\qquad n=0,1,2,...$
(44)

###  5.2 Stokes problem

We follow a similar ’direct’ approach for the Stokes problem and use the
following integral representation for the velocity,

$\mathbf{u}_{S}(\mathbf{y})=\mathbf{u}^{\infty}-\frac{1}{8\pi}\int_{\partial\Omega_{S}}S(\mathbf{y,x})\,\mathbf{f}_{S}(\mathbf{x})\,dS(\mathbf{x})+\frac{1}{8\pi}\int_{\partial\Omega_{S}}\mathbf{u}_{S}(\mathbf{x})\cdot
T(\mathbf{y,x})\cdot\mathbf{n}_{S}(\mathbf{x})\,dS(\mathbf{x}),\quad\mathbf{y}\in\Omega_{S},$
(45)

where $\mathbf{n}_{S}$ is the unit outward normal vector to the Stokes domain
boundary, $\mathbf{f}_{S}=\bm{\sigma}_{S}\cdot\mathbf{n}_{S}$ is the surface
force, or traction, on the boundary $\partial\Omega_{S}$, $\bm{\sigma}_{S}$ is
the stress tensor, and

$S_{ij}(\mathbf{y,x})=\frac{\delta_{ij}}{|\mathbf{\hat{y}}|}+\frac{\hat{y}_{i}\
\hat{y}_{j}}{|\mathbf{\hat{y}}|^{3}},\qquad\qquad
T_{ijk}(\mathbf{y,x})=-6\frac{\hat{y}_{i}\hat{y}_{j}\hat{y}_{k}}{|\mathbf{\hat{y}}|^{5}},$
(46)

are the Stokeslet and stresslet fundamental solutions of the Stokes equations,
where $\mathbf{\hat{y}}=\mathbf{y}-\mathbf{x}$, $\delta_{ij}$ is the Kronecker
delta, and $i,j,k=1,2,3$ are Cartesian coordinates. A background flow
$\mathbf{u}^{\infty}$ is added for external flows.

The representation (45) is based on the Lorentz reciprocal theorem, and again,
the first integral is called the single-layer potential and the second
integral is the double-layer potential of Stokes flow. The double-layer
undergoes a discontinuity across the boundary much like (41), and we again use
jump conditions to continuously extend the representation (45) to the
boundary:

$\frac{1}{2}\mathbf{u}_{S}(\mathbf{y})=\mathbf{u}^{\infty}-\frac{1}{8\pi}\int_{\partial\Omega_{S}}S(\mathbf{y,x})\,\mathbf{f}_{S}(\mathbf{x})\,dS(\mathbf{x})+\frac{1}{8\pi}\int_{\partial\Omega_{S}}\mathbf{u}_{S}(\mathbf{x})\cdot
T(\mathbf{y,x})\cdot\mathbf{n}_{S}(\mathbf{x})\,dS(\mathbf{x}),\quad\mathbf{y}\in\partial\Omega_{S}.$
(47)

We write (47) as

$\frac{1}{2}\mathbf{u}_{S}=K\mathbf{u}_{S}+b_{S},\qquad\qquad
b_{S}=-M\mathbf{f}_{S}+\mathbf{u}^{\infty},$ (48)

where by $M$ and $K$ we denote the integral operators corresponding to the
Stokes single- and double-layer potentials, respectively. Assuming a known
surface force, equation (48) is a second-kind integral equation for
$\mathbf{u}_{S}$, and we solve it by successive approximations

$\mathbf{u}_{S}^{(n+1)}=\frac{1}{2}\mathbf{u}_{S}^{(n)}+K\mathbf{u}_{S}^{(n)}+b_{S},\qquad
n=0,1,2,...$ (49)

It is a well-known challenge in the numerical solution of the integral
equations (44) and (49) that singularities develop in the integrands
$H,L,K,M$, and we address this issue by using a high accuracy regularization
technique which we describe next.

###  5.3 High-order regularization of kernels on the surface

The integrands in (42) and (47) exhibit singularities for
$\mathbf{y}=\mathbf{x}$. To address this issue, we first use subtraction in
the double layer integrals to reduce the singularities,

$\displaystyle(Hp_{D})(\mathbf{y})$
$\displaystyle=\int_{\partial\Omega_{D}}[p_{D}(\mathbf{x})-p_{D}(\mathbf{y})]\,\frac{\partial
G(\mathbf{y,x})}{\partial
n(\mathbf{x})}\,dS(\mathbf{x})+\frac{1}{2}p_{D}(\mathbf{y}),\quad\mathbf{y}\in\partial\Omega_{D},$
(50a) $\displaystyle(K\mathbf{u}_{S})(\mathbf{y})$
$\displaystyle=-\frac{1}{8\pi}\int_{\partial\Omega_{S}}[\mathbf{u}_{S}(\mathbf{x})-\mathbf{u}_{S}(\mathbf{y})]\cdot
T(\mathbf{y,x})\cdot\mathbf{n}(\mathbf{x})dS(\mathbf{x})-\frac{1}{2}\mathbf{u}_{S}(\mathbf{y}),\quad\mathbf{y}\in\partial\Omega_{S}.$
(50b)

Next, we apply the regularization method of [4, 5] for Laplace’s equation, and
recently developed for the Stokes equations in [54]. The approach is to
replace $1/r^{p}$, where $r=|\mathbf{\hat{y}}|$, by a smooth version
$s(r/\delta)/r^{p}$, where $s(r/\delta)$ is a radial smoothing function and
$\delta>0$ is a smoothing, or regularization, parameter. As we are interested
in evaluating the integrals on the boundary, the smoothing factor $s$ can be
chosen to have high order, so that the accuracy in solving the integral
equations approaches $O(h^{5})$, where $h$ is the grid size in the coordinate
planes. The derivation of these smoothing functions can be found in [4, 5,
54]. For completeness, we state the regularization formulas for both the Darcy
and Stokes problems. The kernels in the Darcy representation (42) are replaced
with their regularized versions,

$G^{\delta}(\mathbf{y,x})=-\frac{s_{1}(r/\delta)}{4\pi r},\qquad\qquad\nabla
G^{\delta}(\mathbf{y,x})=\nabla
G^{\delta}(\mathbf{y,x})s_{2}(r/\delta)=\frac{\mathbf{\hat{y}}}{4\pi
r^{3}}s_{2}(r/\delta).$ (51)

The regularized single and double layer Stokes kernels in the representation
(47) are

$S^{\delta}_{ij}(\mathbf{y,x})=\delta_{ij}\frac{s_{1}(r/\delta)}{r}+\hat{y}_{i}\
\hat{y}_{j}\frac{s_{3}(r/\delta)}{r^{3}},\qquad\qquad
T^{\delta}_{ijk}(\mathbf{y,x})=-6\hat{y}_{i}\hat{y}_{j}\hat{y}_{k}\frac{s_{4}(r/\delta)}{r^{5}}.$
(52)

The radial smoothing functions $s_{i}(r/\delta)$ are chosen to yield high
order accuracy when the integrals are evaluated on the surface, with the error
approaching $O(h^{5})$ in discretization size $h$ as shown in [4, 5, 54],

$\displaystyle s_{1}(r)$
$\displaystyle=\operatorname{erf}(r)-\frac{2}{3\sqrt{\pi}}r(2r^{2}-5)e^{-r^{2}},$
(53a) $\displaystyle s_{2}(r)$
$\displaystyle=\operatorname{erf}(r)+\frac{2}{3\sqrt{\pi}}r(2r^{2}-3)e^{-r^{2}},$
(53b) $\displaystyle s_{3}(r)$
$\displaystyle=\operatorname{erf}(r)-\frac{2}{3\sqrt{\pi}}r(4r^{4}-14r^{2}+3)e^{-r^{2}},$
(53c) $\displaystyle s_{4}(r)$
$\displaystyle=\operatorname{erf}(r)-\frac{2}{9\sqrt{\pi}}r(8r^{6}-36r^{4}+6r^{2}+9)e^{-r^{2}},$
(53d)

where $\operatorname{erf}$ is the error function. This approach is quite
simple in that it does not require any special quadrature for evaluating the
integrals on the boundary.

###  5.4 Quadrature

With the integrands smoothed out, to discretize the integrals a quadrature
rule for closed surfaces from [59, 6] is applied, and we briefly describe it
here. First, an angle $\theta$ is chosen and a partition of unity on the unit
sphere is defined, consisting of functions $\psi_{1},\psi_{2},\psi_{3}$ with
$\Sigma_{i}\psi_{i}\equiv 1$ such that $\psi_{i}(\mathbf{n})=0$ if
$|\mathbf{n}\cdot\mathbf{e}_{i}|\leq\cos{\theta}$, where $\mathbf{e}_{i}$ is
the $i$th coordinate vector. Here we use $\theta=70^{o}$. For mesh size $h$, a
set $R_{3}$ of quadrature points consists of points $\bf{x}$ on the surface of
the form $(j_{1}h,j_{2}h,x_{3})$ such that
$|\mathbf{n(x})\cdot\mathbf{e}_{3}|\geq\cos{\theta}$, where
$\mathbf{n}(\mathbf{x})$ is the unit normal at $\mathbf{x}$, see Fig. 3. Sets
$R_{1}$ and $R_{2}$ are defined similarly. For a function $f$ on the surface
the integral is computed as

$\int_{S}f(\mathbf{x})\,dS(\mathbf{x})\approx\sum_{i=1}^{3}\sum_{\mathbf{x}\in
R_{i}}\frac{\psi_{i}(\mathbf{n}(\mathbf{x}))\,f(\mathbf{x})}{|\mathbf{n}(\mathbf{x})\cdot\mathbf{e}_{i}|}\,h^{2}.$
(54)

The partition of unity functions $\psi_{i}$ are constructed from the
$C^{\infty}$ bump function $b(r)=e^{r^{2}/(r^{2}-1)}$ for $|r|<1$ and zero
otherwise. The quadrature is effectively reduced to the trapezoidal rule
without boundary. Thus for regular integrands the quadrature has arbitrarily
high order accuracy, limited only by the degree of smoothness of the integrand
and surface. The points in $R_{i}$ can be found by a line search since they
are well separated; see [59] and [6].

Suppose the surface is given by $\phi(x_{1},x_{2},x_{3})=0$, with $\phi>0$
outside, and the normal vector $\nabla\phi$ is predominantly in the $x_{3}$
direction. Then the parameterization is given by $x_{3}=z(x_{1},x_{2})$, where
$z$ is the vertical coordinate on the surface. Derivatives of this
parameterization are computed implicitly by $z_{i}\equiv\partial z/\partial
x_{i}=-\phi_{i}/\phi_{3}$, the tangent vectors are defined as
$\mathbf{T}_{1}=(1,0,z_{1})$ and $\mathbf{T}_{2}=(0,1,z_{2})$. The outward
normal is $\mathbf{n}=\nabla\phi/|\nabla\phi|$ or
$\mathbf{n}=\pm(-z_{1},-z_{2},1)/\sqrt{1+z_{1}^{2}+z_{2}^{2}}$, where ”+” is
used if $x_{3}>z(x_{1},x_{2})$ outside and ”-” is used otherwise.

In Figure 3, we show the set $R_{3}$ of quadrature points, given by the top
and bottom coordinate grids $(j_{1}h,j_{2}h,x_{3})$, for the following three
surfaces:

$\displaystyle\phi(x_{1},x_{2},x_{3})$
$\displaystyle=x_{1}^{2}+x_{2}^{2}+x_{3}^{2}-1,$ (55a)
$\displaystyle\phi(x_{1},x_{2},x_{3})$
$\displaystyle=\frac{x_{1}^{2}}{a^{2}}+\frac{x_{2}^{2}}{b^{2}}+\frac{x_{3}^{2}}{c^{2}}-1,$
(55b) $\displaystyle\phi(x_{1},x_{2},x_{3})$
$\displaystyle=\sum_{k=1}^{4}\exp(-|{\bf x}-{\bf x}_{k}|^{2}/r^{2})-c.$ (55c)

For the ellipsoid (55b) we set $a=1,b=0.6,c=0.4$, and for the ’four-atom
molecule’ surface (55c) we use centers
$\mathbf{x}_{1}=(\sqrt{3}/3,0,-\sqrt{6}/12)$,
$\mathbf{x}_{2,3}=(-\sqrt{3}/6,\pm.5,-\sqrt{6}/12)$,
$\mathbf{x}_{4}=(0,0,\sqrt{6}/4)$ and $r=.5$, $c=.6$, as in [6].

Figure 3: Quadrature points generated by $\psi_{3}$ and $h=1/16$ on the
surface of the unit sphere (55a) (top left), the ellipsoid with semiaxes
$a=1,b=0.6,c=0.4$ (55b) (top right), and the ’molecule’ (55c) (bottom row).
The quadrature points are at the intersections of the lines.

This quadrature rule works well for general surfaces without requiring
coordinate charts, using projections on coordinate planes instead. As the
weights for the quadrature points are found from a partition of unity on the
unit sphere, they do not depend on the particular surface. The quadrature
points can be found efficiently if, for example, the surface is given
analytically or numerically as the level set of a function. With this
approach, the formulation is very simple to implement, the data structure and
information needed to describe the boundary are minimal, and there are no
parameters to fine-tune except the regularization parameter $\delta$ (which is
usually taken as a multiple of grid size $h$, see the discussion in [54]).

##  6 Algorithm

Now we summarize the algorithm for solving the Stokes-Darcy system. First, we
note that

$\frac{\partial p_{D}}{\partial
n}=-\frac{\mu}{\kappa}\mathbf{u}_{D}\cdot\mathbf{n}_{D},$ (56)

and by combining the interface conditions (4) and (5),

$\mathbf{f}_{S}=\bm{\sigma}_{S}\cdot\mathbf{n}_{S}=-p_{D}\mathbf{n}_{S}-\frac{\gamma}{\sqrt{\kappa}}[\mathbf{u}_{S}-(\mathbf{u}_{S}\cdot\mathbf{n}_{S})\mathbf{n}_{S}],$
(57)

so we rewrite the Darcy (43) and Stokes (48) integral equations as

$\displaystyle\frac{1}{2}p_{D}$ $\displaystyle=Hp_{D}+b_{D},\qquad\qquad
b_{D}=\frac{\mu}{\kappa}L(\mathbf{u}_{D}\cdot\mathbf{n}_{D}),$ (58a)
$\displaystyle\frac{1}{2}\mathbf{u}_{S}$
$\displaystyle=K\mathbf{u}_{S}+b_{S},\qquad\qquad
b_{S}=-M\mathbf{f}_{S}+\mathbf{u}^{\infty},$ (58b)

We outline the algorithm of the iterative Stokes and Darcy solutions using
successive approximations approach, keeping in mind that the successive
approximations are eventually replaced by GMRES for a more efficient solution.

  1. 1.

Let $k=0$, $\mathbf{u}_{S}^{(0)}=\mathbf{0}$ and
$q^{(0)}=-\mathbf{u}_{S}^{(0)}\cdot\mathbf{n}_{S}=0$.

  2. 2.

For $k=1,2,...$ until convergence, do:

    1. (a)

Compute the Laplace single layer potential $b_{D}$ in (58a) using
$\mathbf{u}_{D}\cdot\mathbf{n}_{D}=q^{(k-1)}$.

    2. (b)

Solve the Darcy problem (58a) for $p_{D}^{(k)}$ via (i) successive
approximations (44) or (ii) GMRES.

    3. (c)

Compute $\mathbf{f}_{S}$ as in (57), using $p_{D}^{(k)}$ and
$\mathbf{u}^{(k-1)}_{S}$.

    4. (d)

Compute the Stokes single layer potential $b_{S}$ in (58b).

    5. (e)

Solve the Stokes problem (58b) for $\mathbf{u}_{S}^{(k)}$ via (i) successive
approximations (49) or (ii) GMRES.

    6. (f)

Set
$q^{(k)}=(1-\theta)q^{(k-1)}-\theta\mathbf{u}_{S}^{(k)}\cdot\mathbf{n}_{S}$.

Here $\theta\in(0,1)$ is a relaxation parameter. We assume convergence when
the relative error between $q^{(k)}$ and $q^{(k-1)}$, which we call the
residual,

$\left(\frac{\sum_{i=1}^{N}\left(q^{(k)}_{i}-q^{(k-1)}_{i}\right)^{2}}{\sum_{i=1}^{N}\left(q^{(k)}_{i}\right)^{2}}\right)^{1/2},$
(59)

falls below a prescribed tolerance.

###  6.1 Discrete iteration operator

To analyze the convergence properties of DDM in the discrete solution, we look
at the iteration matrix obtained from the integral formulation and its
eigenvalues. For simplicity, we assume $\mathbf{u}^{\infty}=0$ and the slip
coefficient $\gamma=0$. With the iteration variable defined as
$q=-\mathbf{u}_{S}\cdot\mathbf{n}_{S}$, each Dirichlet-Neumann iteration
includes solving the following two problems:

$\frac{1}{2}p_{D}^{(k)}=\tilde{H}p_{D}^{(k)}+\frac{\mu}{\kappa}\tilde{L}q^{(k-1)},$
(60)
$\frac{1}{2}\mathbf{u}_{S}^{(k)}=\tilde{K}\mathbf{u}_{S}^{(k)}+\tilde{M}_{n}p_{D}^{(k)},$
(61)

where $\tilde{H},\tilde{K}$ are the matrices corresponding to the double layer
potentials for Laplace and Stokes, respectively, and $\tilde{L},\tilde{M}_{n}$
are the matrices corresponding to the single layer for Laplace and Stokes (the
Stokes single layer $\tilde{M}_{n}$ including an extra multiplication by the
normal). The ’tilde’ symbol denotes the discrete version of the operators,
with the kernels smoothed out as in Section 5.3 and the quadrature in Section
5.4 applied. All these are dense matrices.

We can then write the discrete representation of the iteration operator as
follows:

$\displaystyle q^{(k)}$ $\displaystyle=$
$\displaystyle(1-\theta)q^{(k-1)}-\theta\mathbf{u}_{S}^{(k)}\cdot\mathbf{n}_{S}$
(62) $\displaystyle=$
$\displaystyle(1-\theta)q^{(k-1)}-\theta\left[(I/2-\tilde{K})^{-1}\tilde{M}_{n}p_{D}^{(k)}\right]\cdot\mathbf{n}_{S}$
$\displaystyle=$
$\displaystyle(1-\theta)q^{(k-1)}-\theta\left[(I/2-\tilde{K})^{-1}\tilde{M}_{n}\right]\cdot\mathbf{n}_{S}\,p_{D}^{(k)}$
$\displaystyle=$
$\displaystyle(1-\theta)q^{(k-1)}-\theta\left[(I/2-\tilde{K})^{-1}\tilde{M}_{n}\right]\cdot\mathbf{n}_{S}\,(I/2-\tilde{H})^{-1}\,\frac{\mu}{\kappa}\,\tilde{L}q^{(k-1)}.$

We write (62) as

$q^{(k)}=\tilde{\mathcal{A}}q^{(k-1)},$ (63)

with the iteration matrix

$\tilde{\mathcal{A}}=I-\theta\left\\{I+\frac{\mu}{\kappa}\big{[}(I/2-\tilde{K})^{-1}\tilde{M}_{n}\big{]}\cdot\mathbf{n}_{S}\
\big{[}(I/2-\tilde{H})^{-1}\tilde{L}\big{]}\right\\}.$ (64)

The operator in the first square brackets corresponds to solving the Stokes
problem for $\mathbf{u}_{S}$, while the operator in the second square brackets
corresponds to solving the Darcy problem for $p_{D}$.

##  7 Numerical Results

Here we analyze the algorithm in terms of accuracy, convergence, and
dependence on physical parameter values for several surfaces. The code was
written in C++, compiled with the clang compiler, and all computations were
performed on a MacBook Pro 2.3 GHz and 32 GB RAM, running the Big Sur OS.

###  7.1 BIEM accuracy

We first demonstrate the high achievable accuracy of the integral formulations
on the boundary. We solve the Darcy and Stokes integral equations
independently of each other, using successive approximations (44) and (49),
each with a known exact solution. For the Darcy problem test, we used
$\kappa=1$ and the harmonic function $p_{D}=e^{x}\sin y$ on the sphere of
radius $R=0.8$. For the Stokes problem test, we used the solution similar to
Sec. 4.4. of [54]. On the outside of a sphere of radius 0.8, we assume the
velocity given by a point force singularity of strength $(1,0,0)$, placed at
$(0.2,0,0)$. We then take the solution on the inside of the sphere to be 0, so
that the strengths in the Stokes single and double layer potentials are the
jumps in velocity and surface force across the surface.

Figure 4: Convergence in grid size $h$ of separate Darcy and Stokes problems
solved on the surface of a unit sphere by successive approximations (44) and
(49), with $\delta=3h$.

Figure 4 shows the accuracy of this solution, where error is defined as the
$L_{2}$ norm of the difference in the computed and exact solutions, $p_{D}$
for Darcy and $\mathbf{u}_{S}$ for Stokes. We used a smoothing parameter
$\delta=3h$ here and throughout this paper, where $h$ is the grid size in the
coordinate planes. A larger choice of $\delta$ is needed to ensure the
regularization error is dominant over the discretization error, so that the
total error approaches $O(h^{5})$ (see [54]). As expected, the convergence
rate observed in Fig. 4 is $O(h^{5})$. Other surfaces, such as a thin
ellipsoid, are generally expected to give somewhat larger errors due to the
larger curvature and varied spacing [54].

###  7.2 Discrete iteration operator

Here we validate the numerical implementation of the Dirichlet-Neumann
iterations for the spherical geometry. Similar to the setup used in Section 4,
consider the Darcy problem inside the unit sphere centered at the origin,
coupled with the Stokes problem outside the sphere. Parameters are set to
$\mu=1$, $\gamma=0$, and the unit sphere was discretized using $h=1/16$
leading to $N=4302$ quadrature points. In Figure 5, the eigenvalues of the
discrete iteration matrix $\tilde{\mathcal{A}}$ in (64) are computed and the
first 50 are shown, to compare with the harmonic coefficients
$\mathcal{A}_{n}$ of the iteration operator in Figure 2.

  * •

$\kappa=1$. In the left figure, the eigenvalues are shown for permeability
$\kappa=1$ and varying relaxation parameter $\theta$ values. Again we see that
a larger value of $\theta\in(0,1)$ leads to smaller eigenvalues, where most
cluster around the value $1-\theta$. We then expect the successive
approximations to converge faster for larger $\theta$.

  * •

$\kappa\ll 1$. Next, keeping $\mu=1$ fixed, we compute the eigenvalues for
smaller values of permeability, $\kappa=10^{-2}$ (center figure), and
$\kappa=10^{-4}$ (right figure). For these, we set the relaxation parameter
$\theta=c\kappa$, and vary the value of coefficient $c=0.5,0.75,0.9$ as
before. In all cases, the vast majority of the eigenvalues cluster around 1
(we saw this with the harmonic coefficients as well). As we will show, this is
going to severely inhibit convergence speeds of DDM by successive
approximations, and we will perform the DDM iterations using GMRES instead.

Figure 5: Eigenvalues of the iteration matrix $\tilde{\mathcal{A}}$ in (64) on
the unit sphere for $h=1/16$, for permeability values
$\kappa=1,10^{-2},10^{-4}$.

###  7.3 Benchmark problem: Stokes flow past a porous sphere

In this section we test our algorithm using the problem of uniform flow of a
viscous fluid past a porous spherical shell. Assuming that the fluid in the
permeable sphere $0\leq r\leq R$ obeys Darcy’s law and the fluid outside the
sphere satisfies the Stokes equations, the solution was obtained by Joseph and
Tao [33]. The matching conditions enforced on the sphere are continuity of
pressure and normal velocity, as well as the no-slip condition for the
tangential velocity of the exterior flow. Considering a stationary porous
sphere of radius $R$ in a uniform streaming viscous fluid
$\mathbf{U}=(0,0,U)$, the solution is given in spherical coordinates
$(r,\theta,\phi)$,

$\displaystyle p^{D}$ $\displaystyle=$
$\displaystyle-3\mu\,U\,r\,\cos\theta/(2R^{2}+\kappa),$ (65) $\displaystyle
u^{D}_{r}$ $\displaystyle=$ $\displaystyle
3U\,\kappa\,\cos\theta/(2R^{2}+\kappa),$ (66) $\displaystyle u^{D}_{\theta}$
$\displaystyle=$ $\displaystyle-3U\,\kappa\,\sin\theta/(2R^{2}+\kappa),$ (67)

for the Darcy problem, and

$\displaystyle p^{S}$ $\displaystyle=$
$\displaystyle-3\mu\,R\,U\,\cos\theta/r^{2}(2+\kappa/R^{2}),$ (68)
$\displaystyle u^{S}_{r}$ $\displaystyle=$
$\displaystyle\left\\{\frac{-3RU}{(2+\kappa/R^{2})\,r}\left(1-\frac{R^{2}}{3r^{2}}\left(1+\frac{2\kappa}{R^{2}}\right)\right)+U\right\\}\cos\theta,$
(69) $\displaystyle u^{S}_{\theta}$ $\displaystyle=$
$\displaystyle\left\\{\frac{3RU}{(2+\kappa/R^{2})\,2\,r}\left(1+\frac{R^{2}}{3r^{2}}\left(1+\frac{2\kappa}{R^{2}}\right)\right)-U\right\\}\sin\theta,$
(70)

for the Stokes problem, where $u_{\phi}=0$ due to axial symmetry. The drag on
the porous sphere, defined as the hydrodynamic force exerted on the sphere,

$\mathbf{D}=\int_{\partial\Omega}\mathbf{f}\,dS(\mathbf{x}),$ (71)

where $\mathbf{f}=\bm{\sigma}_{S}\cdot\mathbf{n}$,
$\mathbf{n}=\mathbf{e}_{r}$, was found to equal

$\mathbf{D}=\frac{6\pi\,\mu\,R\,U}{1+\kappa/(2R^{2})}.$ (72)

The usual solution and drag for the flow past a solid sphere are recovered
from these formulas by setting $\kappa=0$.

To test our algorithm against this solution, we compute the Stokes stress
vector as $\mathbf{f}=f_{r}\mathbf{e}_{r}+f_{\theta}\mathbf{e}_{\theta}$ with

$f_{r}=-p_{D}+2\mu\frac{\partial u_{r}}{\partial r},$ (73)
$f_{\theta}=\mu\left(r\frac{\partial}{\partial
r}\left(\frac{u_{\theta}}{r}\right)+\frac{1}{r}\frac{\partial
u_{r}}{\partial\theta}\right),$ (74)

being the normal and tangential components. Note that we replaced $p_{S}$ with
$p_{D}$ in (73) since the exact solution was derived assuming continuous
pressure $p_{D}=p_{S}$ on the boundary. To compute the drag numerically, we
integrate (71) using our usual quadrature, and compare to the exact value
given by (72).

In our tests, we use a unit sphere $R=1$, $U=1,\mu=1$ and no-slip $\gamma=0$.
We use a tolerance of $10^{-9}$ for the residual (59) in the Dirichlet-Neumann
iterations, and for the solution of the integral equations, we use successive
approximations (44) and (49), as well as GMRES, both with tolerance $10^{-9}$.
Table 1 shows the results for $\kappa=1$. For different grid sizes $h$,
displayed are: the number of quadrature points $N$ on the sphere surface, the
number of Dirichlet-Neumann iterations via successive approximations (D-N SA)
for two values of $\theta=0.5$ and $\theta=0.75$, the number of successive
approximations (local SA) needed to solve the Darcy and Stokes problems at
each D-N step to reach the set tolerance, the number of GMRES iterations
(local GMRES) for each problem at each D-N step for the same tolerance, and
finally the errors in the Darcy and Stokes solutions computed in $L_{2}$ norm
as

$||p^{D}_{h}-p^{D}||_{2}=\left(\sum_{i=1}^{N}\left(p_{D}(\mathbf{x}_{i})-p_{D}^{ex}(\mathbf{x}_{i})\right)^{2}\Big{/}N\right)^{1/2},$
(75)
$||\mathbf{u}^{S}_{h}-\mathbf{u}^{S}||_{2}=\left(\sum_{i=1}^{N}\left|\mathbf{u}_{S}(\mathbf{x}_{i})-\mathbf{u}_{S}^{ex}(\mathbf{x}_{i})\right|^{2}\Big{/}N\right)^{1/2},$
(76)

where $p_{D}^{ex}$ and $\mathbf{u}_{S}^{ex}$ are the exact values and
$|\mathbf{u}|$ is the Euclidean norm. For a more detailed view, Figure 6 shows
the DDM residual and the errors in drag, $p_{D}$ and $\mathbf{u}_{S}$, by
comparing the numerical solution to the exact solution at every D-N iteration
for grid sizes $h=1/16$ and $h=1/32$ and varying $\theta$.

We make the following observations:

  * •

D-N convergence and dependence on the relaxation parameter $\theta$. For
$\theta=0.5$, about 19 D-N iterations are required to reach the tolerance,
independent of the grid size, while for $\theta=0.75$ only 7 iterations are
needed for $h=1/16$, and this number reduces quickly to 2 iterations for
smaller grid sizes. Figure 6 (top left) shows the gradual decrease in the
residual for $\theta=0.5$ and a more rapid decrease for $\theta=0.75$. This
faster convergence can be attributed to the smaller spectral radius of the
iteration operator, as was observed for $\kappa=1$ in Figures 2 and 5. This is
remarkable given that with smaller $h$ the local Stokes and Darcy problems
become larger.

  * •

Solution errors. Further, we see that the errors in drag, $p_{D}$ and
$\mathbf{u}_{S}$ in Figure 6 all level off at the BIEM accuracy shown in Table
1 (again, much faster with $\theta=0.75$). From the error values in Table 1 we
see that the convergence rate is $O(h^{5})$ or higher with grid refinement.

  * •

Local problems: SA vs. GMRES. Table 1 shows the number of successive
approximations (local SA) performed at each D-N iteration, as well as the
number of GMRES iterations (local GMRES) as an alternative. For SA, a range is
shown because the number decreases slowly with each D-N iteration. It is clear
that successive approximations do not converge quickly, and GMRES is far
superior in solving the local Stokes and Darcy problems at each iteration.
Both methods require evaluating the double layer potential at each iteration,
and since GMRES performs far fewer double layer evaluations, we use GMRES in
all our remaining tests.

  * •

Local problems: CPU time. Both SA and GMRES require evaluating the double
layer potential at each iteration, which in discretized form is a dense linear
system. The Darcy system is $N\times N$, while the Stokes system is $3N\times
3N$, where $N$ is the number of quadrature points. The computational cost is
then $O(N^{2})$ for each problem. The actual CPU times for each double layer
integral evaluation are shown in Figure 7. Note that the discretized integrals
were implemented to take advantage of the target-source symmetry in all
kernels, reducing the computational time in half. If very large systems are
considered, the computational time could be reduced to $O(N\log N)$ or $O(N)$
using a fast summation algorithm such as a treecode [58] or a fast multipole
method [56].

  * •

Dependence on $\kappa$. We next proceed to study the convergence of D-N
iterations for different values of permeability. For $\kappa=1$, we see a
strong dependence on the values of $\theta$ when successive approximations are
used, because the spectral radius of the iteration operator is smaller for
larger $\theta$. Table 2 reiterates some of these results and also shows the
results for $\kappa=10^{-2}$ and $\kappa=10^{-4}$. For small values of
$\kappa$ such as these, we saw that the spectral radius of the iteration
operator is very close to 1, and this indeed inhibits convergence. The number
of successive approximations needed for D-N convergence (D-N by SA in Table 2)
exceeds 100 because the residual decreases very slowly. This is not improved
by using a larger $\theta$ as it did for $\kappa=1$ (in all cases, we set
$\theta=c\kappa$ with $c\in(0,1)$). Because of this, we perform the D-N
iterations using GMRES, which converges much faster. The D-N by GMRES column
in Table 2 shows the number of iterations for the same tolerance $10^{-9}$ for
the residual. The number of iterations is 7 for $\kappa=10^{-2}$ and 8 for
$\kappa=10^{-4}$ for the coarsest grid $h=1/16$. As the grid size is reduced
and the solution becomes more accurate, the number of D-N iterations by GMRES
reduces. For the solution errors, which are also shown in the left graph in
Figure 8, we observe convergence of the order $O(h^{5})$.

$h$ | $N$ | D-N SA | local SA | local GMRES | $||p^{D}_{h}-p^{D}||_{2}$ | $||\mathbf{u}^{S}_{h}-\mathbf{u}^{S}||_{2}$
---|---|---|---|---|---|---
$\theta=0.5$ | $\theta=0.75$ | Darcy | Stokes | Darcy | Stokes
1/16 | 4302 | 19 | 7 | 37-2 | 26-2 | 5 | 8 | 3.450e-05 | 1.053e-04
1/32 | 17070 | 19 | 5 | 37-1 | 26-1 | 4 | 7 | 1.442e-06 | 5.525e-06
1/64 | 68166 | - | 2 | - | - | 3 | 5 | 2.605e-08 | 1.716e-07
1/128 | 272718 | - | 1 | - | - | 2 | 3 | 3.236e-10 | 4.735e-09
Table 1: Flow past a porous sphere with $\kappa=1$, DDM residual tolerance =
$10^{-9}$, local solution tolerance = $10^{-9}$.

Figure 6: Flow past a porous sphere with $\kappa=1$, tolerance $10^{-9}$ for
residual/local solutions. Top left: DDM residual, top right: error in drag,
bottom left: error in $p_{D}$, bottom right: error in $\mathbf{u}_{S}$.

Figure 7: Flow past a porous sphere with $\kappa=1$. CPU time in seconds on each evaluation of the double layer potential (one successive approximation or one GMRES iteration). $\kappa$ | $h$ | D-N by SA | D-N by GMRES | local GMRES | $||p^{D}_{h}-p^{D}||_{2}$ | $||\mathbf{u}^{S}_{h}-\mathbf{u}^{S}||_{2}$
---|---|---|---|---|---|---
$\theta=0.5\kappa$ | $\theta=0.75\kappa$ | $\theta=0.5\kappa$ | $\theta=0.75\kappa$ | Darcy | Stokes
$1$ | 1/16 | 19 | 7 | 4 | 4 | 5 | 8 | 3.450e-05 | 1.053e-04
| 1/32 | 19 | 5 | 4 | 4 | 6 | 6 | 1.442e-06 | 5.525e-06
| 1/64 | - | 2 | 3 | 3 | 6 | 6 | 2.605e-08 | 1.716e-07
$10^{-2}$ | 1/16 | $>100$ | $>100$ | 7 | 7 | 6 | 9 | 9.484e-05 | 4.500e-05
| 1/32 | $>100$ | $>100$ | 6 | 7 | 7 | 7 | 3.498e-06 | 2.576e-06
| 1/64 | - | - | 4 | 4 | 8 | 7 | 3.084e-08 | 8.546e-08
$10^{-4}$ | 1/16 | $>100$ | $>100$ | 8 | 8 | 8 | 11 | 5.474e-04 | 4.243e-05
| 1/32 | $>100$ | $>100$ | 2 | 2 | 8 | 9 | 3.490e-06 | 2.554e-06
| 1/64 | - | - | 1 | 1 | 5 | 5 | 3.899e-08 | 8.451e-08
Table 2: Flow past a porous sphere with varying $\kappa$, tolerance $10^{-9}$
for residual/local solutions.

###  7.4 Thin porous ellipsoid

Next we test the method on a geometry that has a larger curvature, a thin
ellipsoid defined in (55b) with $a=1,b=0.6,c=0.4$. In this case the exact
solution is not known, and we compute the error empirically as
$||p_{h}-p_{h/2}||_{2}$ and $||\mathbf{u}_{h}-\mathbf{u}_{h/2}||_{2}$, with
similar definitions of the $L_{2}$ norm as in (75)-(76). Table 3 shows the
results for permeability values $\kappa=1,10^{-2},10^{-4},10^{-7}$. As we saw
in the sphere example, the number of DDM iterations increases slightly when we
go from $\kappa=1$ to $\kappa<1$. Unlike the sphere, however, as the grid is
refined, the number of DDM iterations does not decrease but rather increases
slightly. Furthermore, for the small value of the permeability
$\kappa=10^{-4}$, the DDM iterations do not converge to the prescribed
tolerance $10^{-9}$ quickly. In the number of iterations shown in the table
(denoted by *), the relative residual reaches a tolerance of about $10^{-6}$,
with the consecutive residuals differing by less than $10^{-9}$. As expected
for the thin ellipsoid, the solution accuracy is lower than what we observed
with the sphere.

$\kappa$ | $h$ | $N$ | D-N by GMRES | local GMRES | $||p^{D}_{h}-p^{D}_{h/2}||_{2}$ | $||\mathbf{u}^{S}_{h}-\mathbf{u}^{S}_{h/2}||_{2}$
---|---|---|---|---|---|---
$\theta=0.5\kappa$ | $\theta=0.75\kappa$ | Darcy | Stokes
$1$ | 1/16 | 1742 | 5 | 5 | 7 | 11 | 3.731e-04 | 6.052e-04
| 1/32 | 6902 | 5 | 5 | 6 | 10 | 3.486e-05 | 4.399e-05
| 1/64 | 27566 | 5 | 5 | 5 | 9 | 3.300e-06 | 3.589e-06
| 1/128 | 110250 | 5 | 5 | 5 | 8 | - | -
$10^{-2}$ | 1/16 | 1742 | 9 | 9 | 9 | 14 | 2.378e-02 | 5.744e-03
| 1/32 | 6902 | 12 | 12 | 9 | 14 | 1.772e-03 | 1.797e-03
| 1/64 | 27566 | 12 | 12 | 8 | 11 | 8.347e-05 | 8.423e-05
| 1/128 | 110250 | 12 | 12 | 8 | 11 | - | -
$10^{-4}$ | 1/16 | 1742 | 9* | 8* | 11 | 16 | 1.689e-01 | 1.873e-02
| 1/32 | 6902 | 11* | 11* | 11 | 17 | 2.913e-02 | 2.914e-02
| 1/64 | 27566 | 13* | 13* | 11 | 17 | 1.435e-03 | 1.436e-03
| 1/128 | 110250 | 13* | 13* | 11 | 18 | - | -
$10^{-7}$ | 1/16 | 1742 | 6 | 6 | 14 | 20 | 1.900e-01 | 1.865e-02
| 1/32 | 6902 | 8 | 8 | 11 | 20 | 3.406e-02 | 3.408e-02
| 1/64 | 27566 | 10 | 10 | 11 | 18 | 1.448e-03 | 1.449e-03
| 1/128 | 110250 | 10 | 10 | 12 | 18 | - | -
Table 3: Flow past a porous ellipsoid with varying $\kappa$, tolerance
$10^{-9}$ for residual/local solutions.

###  7.5 Porous molecule

The four-atom molecule surface used in these tests is defined in (55c). Again,
we compute the errors empirically similar to the ellipsoid. The results for
different permeability values and grid sizes are shown in Table 4. The method
performs well for this surface, converging to the tolerance in just a few
iterations, with the speed of convergence improving when the discretization of
the surface is refined. For $h=1/128$ for example, which results in $N=272718$
quadrature points, only 2 DDM iterations are needed for $\kappa=1$ and
$\kappa=10^{-2}$, while the $\kappa=10^{-4}$ and $\kappa=10^{-7}$ cases only
require one iteration for the tolerance of $10^{-9}$. The solution accuracy is
higher than for the ellipsoid and similar to the sphere accuracy, with the
order of convergence as predicted $O(h^{5})$ \- this is shown in Figure 8.

$\kappa$ | $h$ | $N$ | D-N by GMRES | local GMRES | $||p^{D}_{h}-p^{D}_{h/2}||_{2}$ | $||\mathbf{u}^{S}_{h}-\mathbf{u}^{S}_{h/2}||_{2}$
---|---|---|---|---|---|---
$\theta=0.5\kappa$ | $\theta=0.75\kappa$ | Darcy | Stokes
$1$ | 1/16 | 4302 | 4 | 4 | 5 | 7 | 1.793e-05 | 1.221e-05
| 1/32 | 17070 | 4 | 4 | 6 | 6 | 8.947e-07 | 9.358e-07
| 1/64 | 68166 | 3 | 3 | 6 | 6 | 1.892e-08 | 1.974e-08
| 1/128 | 272718 | 2 | 2 | 6 | 6 | - | -
$10^{-2}$ | 1/16 | 4302 | 8 | 8 | 6 | 9 | 1.689e-04 | 3.228e-05
| 1/32 | 17070 | 7 | 7 | 7 | 7 | 6.102e-06 | 6.174e-06
| 1/64 | 68166 | 4 | 5 | 7 | 7 | 5.877e-08 | 6.263e-08
| 1/128 | 272718 | 2 | 2 | 7 | 7 | - | -
$10^{-4}$ | 1/16 | 4302 | 8 | 11 | 8 | 10 | 1.052e-03 | 2.158e-05
| 1/32 | 17070 | 2 | 8 | 9 | 9 | 5.883e-06 | 5.981e-06
| 1/64 | 68166 | 1 | 1 | 5 | 5 | 7.676e-08 | 8.000e-08
| 1/128 | 272718 | 1 | 1 | 3 | 3 | - | -
$10^{-7}$ | 1/16 | 4302 | 1 | 1 | 8 | 12 | 1.330e-04 | 4.356e-05
| 1/32 | 17070 | 1 | 1 | 8 | 10 | 5.710e-06 | 5.797e-06
| 1/64 | 68166 | 1 | 1 | 7 | 8 | 7.679e-08 | 8.002e-08
| 1/128 | 272718 | 1 | 1 | 3 | 3 | - | -
Table 4: Flow past a porous molecule with varying $\kappa$, tolerance
$10^{-9}$ for residual/local solutions.

Figure 8: Solution errors for three surfaces: unit sphere, thin ellipsoid,
four-atom molecule.

##  8 Conclusions and future work

We have presented a domain decomposition method of a sequential Dirichlet-
Neumann type to solve the Stokes-Darcy system of PDEs. The method is based on
splitting the interface conditions so that first, the Darcy problem is solved
with one of the conditions, and then the Stokes problem is solved with the
other two. The information is exchanged at each step of this iterative
procedure. The local Stokes and Darcy problems are solved by second kind
boundary integral formulations. A regularization technique specifically
designed to achieve high accuracy on the boundary is applied, and the
quadrature method of [59, 6] is used, which does not use coordinate charts and
is simple to implement. Numerical results validate the algorithm against a
problem with a known analytical solution of viscous flow past a porous sphere,
as well as demonstrate the applicability of the method to more general
geometries.

Our convergence analysis of the method shows strong dependence of the spectral
radius of the iteration operator on the permeability $\kappa$. For small
values of $\kappa$, the eigenvalues cluster close to 1, inhibiting convergence
of successive approximations. While we demonstrate that using a Krylov
subspace method such as GMRES mitigates this issue, more robust DDMs that are
less dependent on the physical parameters are desirable and will be
investigated. It was shown in the two-dimensional case [9] that a Robin-Robin
type domain decomposition method with non-local operators (based on the
boundary integral formulations) can be chosen in the de-coupling of
transmission conditions to dramatically improve the convergence properties of
the solution. The dependence of convergence rates on the physical parameters
of viscosity $\mu$ and permeability $\kappa$ can be reduced as well. Similar
approaches will be investigated in the three-dimensional case.

In many relevant applications, the porous medium has a heterogeneous
structure, and the BIEM is not a suitable choice to model the flow there. This
work can be extended to model such porous media, by coupling the BIEM in the
domain of free flow with more appropriate techniques such as the finite
element method (FEM) in the porous domain with variable properties. With this
approach, the FEM can handle the porous domain where the heterogeneities
occur, while the BIEM will efficiently deal with the Stokes domain, which can
be large or even infinite.

To limit the scope of this paper, only the steady state problem was considered
with a single interface. In future work, time-dependent problems and multi-
surface cases will be considered. In such a simulation, several surfaces can
be close to each other, and it is well-known that the integrands become nearly
singular and therefore numerically non-trivial to integrate with sufficient
accuracy. The correction methods of [4, 5, 54] will be used to improve the
accuracy in the nearly singular case. Another challenge in such a simulation
is the high computational cost of the direct evaluation of the resulting dense
linear systems, which is $O(N^{2})$, where $N$ is the number of unknowns. To
be able to simulate large enough systems, the computational cost will be
reduced using fast summation techniques such as the kernel-independent
treecode [58] or a fast multipole method [56].

## Acknowledgments

The author gratefully acknowledges support by the National Science Foundation
under grant DMS-2012371.

## References

  * [1] I. Ambartsumyan, E. Khattatov, T. Nguyen, and I. Yotov.  Flow and transport in fractured poroelastic media.  Int. J. Geomath., 10(1):1–34, 2019.
  * [2] T. Arbogast and D.S. Brunson.  A computational method for approximating a Darcy-Stokes system governing a vuggy porous medium.  Comput Geosci, 11:207–218, 2007.
  * [3] K. E. Atkinson.  The numerical solution of boundary integral equations.  In I. Duff and G. Watson, editors, The State of the Art in Numerical Analysis, pages 223–259. Clarendon Press, 1997.
  * [4] J. T. Beale.  A convergent boundary integral method for three-dimensional water waves.  Math. Comp., 70:977–1029, 2001.
  * [5] J. T. Beale.  A grid-based boundary integral method for elliptic problems in three dimensions.  SIAM J. Numer. Anal., 42:599–620, 2004.
  * [6] J. T. Beale, W. Ying, and J. R. Wilson.  A simple method for computing singular or nearly singular integrals on closed surfaces.  Commun. Comput. Phys., 20(3):733–753, 2016.
  * [7] G. S. Beavers and D. D. Joseph.  Boundary conditions at a naturally permeable wall.  J. Fluid Mech., 30:197–207, 1967.
  * [8] Y. Boubendir and S. Tlupova.  Stokes-Darcy boundary integral solutions using preconditioners.  J. Comput. Phys., 228(23):8627–8641, 2009.
  * [9] Y. Boubendir and S. Tlupova.  Domain decomposition methods for solving Stokes-Darcy problems with boundary integrals.  SIAM J. Sci. Comput., 35:B82–B106, 2013.
  * [10] M. Bukač, I. Yotov, R. Zakerzadeh, and P. Zunino.  Partitioning strategies for the interaction of a fluid with a poroelastic material based on a Nitsche’s coupling approach.  Comp. Meth. Appl. Mech. and Engin., 292(1):138–170, 2015.
  * [11] E. Burman and P. Hansbo.  A unified stabilized method for Stokes’ and Darcy’s equations.  J. Comput. Appl. Math., 198:35–51, 2007.
  * [12] Y. Cao, M. Gunzburger, X. He, and X. Wang.  Parallel, non-iterative, multi-physics domain decomposition methods for time-dependent Stokes-Darcy systems.  Math. Comp., 83(288):1617–1644, 2014.
  * [13] Y. Cao, M. Gunzburger, F. Hua, and X. Wang.  Coupled Stokes-Darcy model with Beavers-Joseph interface boundary conditions.  Comm. Math. Sci., 8(1):1–25, 2010.
  * [14] W. Chen, M. Gunzburger, F. Hua, and X. Wang.  A parallel Robin-Robin domain decomposition method for the Stokes-Darcy system.  SIAM J. Numer. Anal., 49(3):1064–1084, 2011.
  * [15] P. Chidyagwai, S. Ladenheim, and D.B. Szyld.  Constraint preconditioning for the coupled Stokes-Darcy system.  SIAM J. Sci. Comput., 38(2):A668–A690, 2016.
  * [16] D. Colton and R. Kress.  Inverse Acoustic and Electromagnetic Scattering Theory.  Springer-Verlag, Berlin, 2nd edition, 1998.
  * [17] M.R. Correa and A.F.D. Loula.  A unified mixed formulation naturally coupling Stokes and Darcy flows.  Comp. Meth. Appl. Mech. and Engin., 198:2710–2722, 2009.
  * [18] R. Cortez.  The method of regularized Stokeslets.  SIAM J. Sci. Comput., 23(4):1204–1225, 2001.
  * [19] R. Cortez, L. Fauci, and A. Medovikov.  The method of regularized Stokeslets in three dimensions: Analysis, validation, and application to helical swimming.  Phys. Fluids, 17:1–14, 2005.
  * [20] M. Discacciati and L. Gerardo-Giorda.  Optimized Schwarz methods for the Stokes-Darcy coupling.  IMA J. Numer. Anal., 38:1959–1983, 2018.
  * [21] M. Discacciati, E. Miglio, and A. Quarteroni.  Mathematical and numerical models for coupling surface and groundwater flows.  Appl. Numer. Math., 43:57–74, 2002.
  * [22] M. Discacciati, A. Quarteroni, and A. Valli.  Robin-Robin domain decomposition methods for the Stokes-Darcy coupling.  SIAM J. Numer. Anal., 45(3):1246–1268, 2007.
  * [23] L. Elasmi and F. Feuillebois.  Green function for a Stokes flow near a porous slab.  Z. Angew. Math. Mech., 81(11):743–752, 2001.
  * [24] J. Galvis and M. Sarkis.  Balancing domain decomposition methods for mortar coupling Stokes-Darcy systems.  Lecture notes in Computational Science and Engineering, 2006.
  * [25] L. Greengard, M. C. Kropinski, and A. Mayo.  Integral equation methods for Stokes flow and isotropic elasticity in the plane.  J. Comp. Phys., 125:403–414, 1996.
  * [26] J. K. Guest and J. H. Prevost.  Topology optimization of creeping fluid flows using a Darcy-Stokes finite element.  Int. J. Numer. Meth. Engng., 66:461–484, 2006.
  * [27] N. Hanspal, V. Nassehi, and A. Kulkarni.  Three-dimensional finite element modelling of coupled free/porous flows: applications to industrial and environmental flows.  Int. J. Numer. Meth. Fluids, 71:1382–1421, 2013.
  * [28] N. Hanspal, A. Waghode, V. Nassehi, and R. Wakeman.  Numerical analysis of coupled Stokes/Darcy flows in industrial filtrations.  Transport in Porous Media, 64:73–101, 2006.
  * [29] J. Happel and H. Brenner.  Low Reynolds number hydrodynamics.  Springer, 1st edition, 1983.
  * [30] G. C. Hsiao and W. L. Wendland.  Boundary Integral Equations.  Springer-Verlag, Berlin, 1st edition, 2008.
  * [31] W. Jager and A. Mikelić.  On the boundary conditions at the contact interface between a porous medium and a free fluid.  Ann. Scuola Norm. Sup. Pisa Cl. Sci., 23:403–465, 1996.
  * [32] W. Jager and A. Mikelić.  On the interface boundary condition of Beavers, Joseph, and Saffman.  SIAM J. Appl. Math., 60:1111–1127, 2000.
  * [33] D. D. Joseph and L. N. Tao.  The effect of permeability on the slow motion of a porous sphere in a viscous liquid.  Z. Angew. Math. Mech., 44(8-9):361–364, 1964.
  * [34] S. Kim and S. J. Karrila.  Microhydrodynamics: Principles and Selected Applications.  Dover, 2005.
  * [35] R. Kress.  Linear Integral Equations.  Springer-Verlag, New York, 2nd edition, 1999.
  * [36] H. Lamb.  Hydrodynamics.  Cambridge UP, 6th edition, 1932.
  * [37] W. J. Layton, F. Schieweck, and I. Yotov.  Coupling fluid flow with porous media flow.  SIAM J. Numer. Anal., 40(6):2195–2218, 2003.
  * [38] R. Li, Y. Gao, J. Li, and Z. Chen.  A weak Galerkin finite element method for a coupled Stokes-Darcy problem on general meshes.  J. Comp. and Appl. Math., 334:111–127, 2018.
  * [39] J. A. Liggett and P. L-F. Liu.  The Boundary Integral Equation Method for Porous Media Flow.  George Allen and Unwin, London, 1983.
  * [40] M. F. Lough, S. H. Lee, and J. Kamath.  An efficient boundary integral formulation for flow though fractured porous media.  J. Comput. Phys., 143:462–483, 1998.
  * [41] K. A. Mardal, X.-C. Tai, and R. Winther.  A robust finite element method for Darcy-Stokes flow.  SIAM J. Numer. Anal., 40:1605–1631, 2002.
  * [42] A. Masud and T. J. R. Hughes.  A stabilized mixed finite element method for Darcy flow.  Comput. Methods Appl. Mech. Eng., 191:4341–4370, 2002.
  * [43] J. A. Ochoa-Tapia and S. Whitaker.  Momentum transfer at the boundary between a porous medium and a homogeneous fluid - I. Theoretical development.  Int. J. Heat Mass Transfer, 38:2635–2646, 1995.
  * [44] J. A. Ochoa-Tapia and S. Whitaker.  Momentum transfer at the boundary between a porous medium and a homogeneous fluid - II. Comparison with experiment.  Int. J. Heat Mass Transfer, 38:2647–2655, 1995.
  * [45] C. Pozrikidis.  Boundary Integral and Singularity Methods for Linearized Viscous Flow.  Cambridge University Press, Cambridge, 1992.
  * [46] C. Pozrikidis and D. A. Farrow.  A model for fluid flow in solid tumors.  Annals Biomed. Eng., 31:181–194, 2003.
  * [47] A. Quarteroni and A. Valli.  Domain Decomposition Methods for Partial Differential Equations.  Oxford University Press, 1999.
  * [48] G. P. Raja Sekhar and T. Amaranath.  Stokes flow inside a porous spherical shell.  Z. angew. Math. Phys., 51:481–490, 2000.
  * [49] J. Rungamornrat and M. F. Wheeler.  Weekly singular integral equations for Darcy’s flow in anisotropic porous media.  Eng. Analysis with Boundary Elements, 30:237–246, 2006.
  * [50] P. G. Saffman.  On the boundary condition at the surface of a porous medium.  Stud. Appl. Math., L(2):93–101, 1971.
  * [51] M. Sahraoui and M. Kaviany.  Slip and no-slip velocity boundary conditions at interface of porous, plain media.  Int. J. Heat Mass Transfer, 35:927–943, 1992.
  * [52] M.-C. Shiue, K. C. Ong, and M.-C. Lai.  Convergence of the MAC scheme for the Stokes-Darcy coupling problem.  J. Sci. Comput., 76(2):1216–1251, 2018.
  * [53] Y. Sun, W. Sun, and H. Zheng.  Domain decomposition method for the fully-mixed Stokes-Darcy coupled problem.  Comput. Methods Appl. Mech. Engrg., 374:113578, 2021.
  * [54] S. Tlupova and J.T. Beale.  Regularized single and double layer integrals in 3D Stokes flow.  J. Comput. Phys., 386(1):568–584, 2019.
  * [55] S. Tlupova and R. Cortez.  Boundary integral solutions of coupled Stokes and Darcy flows.  J. Comput. Phys., 228(1):158–179, 2009.
  * [56] A.-K. Tornberg and L. Greengard.  A fast multipole method for the three-dimensional Stokes equations.  J. Comput. Phys., 227(1):1613–1619, 2008.
  * [57] D. Vassilev, C. Wang, and I. Yotov.  Domain decomposition for coupled Stokes and Darcy flows.  Comput. Methods Appl. Mech. Engrg., 268:264–283, 2014.
  * [58] L. Wang, R. Krasny, and S. Tlupova.  A kernel-independent treecode based on barycentric Lagrange interpolation.  Commun. Comput. Phys., 28(4):1415–1436, 2020.
  * [59] J. R. Wilson.  On computing smooth, singular and nearly singular integrals on implicitly defined surfaces.  PhD thesis, Duke University, 2010.


11institutetext: Arizona State University, Tempe AZ, USA
11email: {bjiang14, psheth5, baoxin.li<EMAIL_ADDRESS>
# CoVaxNet: An Online-Offline Data Repository for COVID-19 Vaccine Hesitancy
Research

Bohan Jiang(✉) Paras Sheth  Baoxin Li  Huan Liu

###### Abstract

Despite the astonishing success of COVID-19 vaccines against the virus, a
substantial proportion of the population is still hesitant to be vaccinated,
undermining governmental efforts to control the virus. To address this
problem, we need to understand the different factors giving rise to such a
behavior, including social media discourses, news media propaganda, government
responses, demographic and socioeconomic statuses, and COVID-19 statistics,
etc. However, existing datasets fail to cover all these aspects, making it
difficult to form a complete picture in inferencing about the problem of
vaccine hesitancy. In this paper, we construct a multi-source, multi-modal,
and multi-feature online-offline data repository
CoVaxNet111https://github.com/jiangbohan/CoVaxNet. We provide descriptive
analyses and insights to illustrate critical patterns in CoVaxNet. Moreover,
we propose a novel approach for connecting online and offline data so as to
facilitate the inference tasks that exploit complementary information sources.

###### Keywords:

Repository Online-Offline Data COVID-19 Vaccine.

##  1 Introduction

The COVID-19 pandemic has killed over six million people and infected 536
million globally as of mid-June, 2022222https://covid19.who.int/. It has been
pointed out that the FDA-authorized COVID-19 vaccines are highly effective at
protecting against severe illness and reinfection[12]. However, according to a
recent KFF survey333https://www.kff.org/coronavirus-covid-19/dashboard/, 32%
of people in the U.S. showed hesitancy to receive at least one dose of the
COVID-19 vaccine. SARS-CoV-2 is constantly mutating to highly contagious
variants and subvariants, causing the number of cumulative cases of infection
continues to grow in most countries. Therefore, the effort to end the pandemic
will be severely hindered if a substantial proportion of the population still
show vaccine hesitancy.

Mitigating the ongoing COVID-19 vaccine hesitancy presents unique challenges.
First, as defined in [15], vaccine hesitancy is a continuum between full
acceptance (pro-vaccine) and outright refusal of all vaccines (anti-vaccine),
and the COVID-19 vaccine hesitancy is driven by a complex set of dynamic
context-specific factors, including social media influence, public stance and
sentiment, structural inequality, risk of the disease, and trust toward
government, etc. All these make the characterization of COVID-19 vaccine
hesitancy a non-trivial task. Second, the research community lacks
comprehensive COVID-19 vaccine datasets containing adequate information to
support effective analysis of the factors mentioned above and their impact on
COVID-19 vaccine hesitancy. Given the complexity of the problem, a data
repository capturing online-offline COVID-19 vaccine information would be
critical to the analysis.

In this paper, we present a multi-source, multi-modal, and multi-feature
online-offline data repository CoVaxNet. As shwon in Figure 1, it contains two
online datasets: (i) a social media dataset and (ii) a fact-checking dataset;
and four offline datasets: (i) COVID-19 statistics, (ii) U.S. Census Bureau
data, (iii) government responses, and (iv) local news reports. To the best of
our knowledge, we are the first to build such a diverse online-offline
COVID-19 vaccine repository. We aim to continuously update this data
repository with new sources and features, as well as maintain completeness.
The main contributions of this work are:

  * $\bullet$

We construct a multi-source, multi-modal, and multi-feature online-offline
repository to facilitate COVID-19 vaccine hesitancy related research such as
(i) detecting COVID-19 vaccine misinformation and stance, (ii) exploring the
effect of online activities on offline outcomes and vice-versa, and (iii)
finding correlations between structural inequality and vaccine hesitancy;

  * $\bullet$

We provide descriptive analyzes of the repository from different perspectives.
We illustrate insights of various features and characteristics in each
dataset, which covers textual, visual, spatio-temporal, and network
information; and

  * $\bullet$

We propose a new online and offline data connection approach with examples
from CoVaxNet and demonstrate baseline performances for stance detection.

Social Media Data Fact-checking Data Online Data Tweets User Profile Text,
Photo, Video Following, Follwer Username, Photo, Bio Location, Creation Date
Favourite, Listed Membership Age, Gender, Race Fact-checking Report Malicious
Media Originated URL Claim, Label Fact-checking Agency Offline DataCOVID-19
Statistics Cases, Deaths Recoveries, Hospitalizations Vaccinations U.S Census
Bureau Data  Population Demographics Vaccine Hesitancy Survey Government
Responses Nationwide/Statewide Policies  URL Local News Title, Content
Pubisher URL CoVaxNet

Figure 1: An overview of CoVaxNet. This figure shows the online-offline data
collected.

##  2 Related Work

Previous studies have investigated vaccine hesitancy on vaccine-preventable
diseases such as HPV, H1N1, and Flu [16, 17]. However, vaccine hesitancy can
vary on different diseases. The SAGE Working Group on Vaccine Hesitancy (WG)
suggests that a complex set of behavioral and social factors should be
considered to construct vaccine hesitancy determinants. [15].

To better understand the COVID-19 vaccine hesitancy, previous works have been
conducted on online data. [6] released the first English Twitter dataset and
an online dashboard about COVID-19 vaccines. [7] presented a large-scale
multilingual Twitter dataset about COVID-19 vaccines. [19] collected 1.8
million anti-vaccine tweets. They investigated the online implicit communities
on topic network using the Louvain algorithm. [13] adopt a human-guided
machine learning framework on social media data to analyze public opinions on
COVID-19 vaccines. [11] proposed a sentiment analysis framework to understand
the public attitude and concerns on mainstream social media toward COVID-19
vaccines.

Researchers also analyzed the COVID-19 vaccine hesitancy on offline data. The
Oxford COVID-19 Government Response Tracker (OxCGRT) involved over 200
volunteers from the Oxford community to collect publicly available global
information on 21 indicators of government response, including data on
vaccination policies [10]. [18] studied the COVID-19 vaccine hesitancy by
analyzing the association between vaccination rates and Social Vulnerability
Index (SVI) from CDC on a county level. [20] examined the geospatial trends in
COVID-19 incidence of the neighborhood-level U.S. Census data. [9] used human
mobility data from Google to analyze the relationship of vaccination rates to
mobility and new COVID-19 cases.

##  3 Construction: What is in CoVaxNet?

###  3.1 Data in the Digital World (Online Data)

Tweets: Inspired by [7], we utilize a snowball sampling method in our crawling
strategy. First, we start with collecting publicly available tweets with the
two most relevant keywords, covid and vaccine [6], as the initial seeds using
Tweepy API444https://www.tweepy.org/. Next, we run n-gram and tf-idf on the
collected tweets to identify the most representative and frequent keywords and
hashtags for pro-COVID-19-vaccine and anti-COVID-19-vaccine stances. Then, we
collect more tweets with newly added keywords and hashtags to expand the list
of seeds. We repeat this process until no new seeds can be extracted. We also
carefully analyze the sampled tweets for each keyword and hashtag to check
whether they are labeled correctly. For example, #Antivaxxers is initially
considered as an anti-vaccine hashtag. However, we find that the sampled
tweets are primarily posted by pro-vaxxers to criticize anti-vaxxers.
Therefore, we eventually categorize #Antivaxxers as one of the pro-vaccine
search keywords. Finally, we remove keywords and hashtags that are not closely
related to the COVID-19 vaccine.

As shown in Table 1, the shortlisted keywords and hashtags are well-classified
to collect COVID-19 pro-vaccine and anti-vaccine tweets. A complete list of
them is available in our data repository, which includes 25 pro-vaccine seeds
and 28 anti-vaccine seeds. In agreement with Twitter’s Terms of Service, we
only release the tweets IDs which can be used to retrieve more objects via
APIs.

Table 1: Keywords and hashtags for collecting pro- and anti-COVID-19-vaccine tweets Stance |  Keywords & Hashtags
---|---
Pro-vaccine |  get vaccinated, vaccine mandate, vaccination work, fully vaccinated, vaccines save lives, vaccinated for covid, #GetVaccinatedOrGetCovid, #Antivacinidiots, #Antivaxxers, #GetYourBooster, … (25)
Anti-vaccine |  no vaccine, f**k vaccines, no forced vaccines, no vaccine mandates, #NoVaccineForMe, #CovidVaccineIsPoison, #NoVaccinePassports, #StopVaccination, #VaccineSideEffects, … (28)

Twitter User Profiles: We collect metadata from Twitter user profiles, as
listed in Figure 1. Due to the lack of geotagged tweets, we leverage the self-
reported locations in the users’ profiles to tag non-geotagged tweets for data
enrichment. Observing that 70.3% of the collected tweets contain either tweet-
level or user-level geo-location information. Furthermore, we use GeoPy [1] to
convert the unstructured self-reported locations to structured geographical
places with coordinates. We also employ state-of-the-art computer vision and
natural language processing tools to infer the demographics of social media
users. We use OpenCV [4] and DeepFace [22] to detect human faces and predict
age, gender, and race on the profile images. For textual features such as
usernames and profile descriptions/bio, by applying Nameparser [3] and
Ethnicolr [21], we extract human names to infer race and gender.

Fact-checking Reports and Malicious URLs: Fact-checking reports can illustrate
how COVID-19-vaccine-related mis- and disinformation, rumors, and conspiracy
theories spread online. In this dataset, we collect 4,263 COVID-19 vaccine-
related fact-checking reports from Poynter555https://www.poynter.org/ifcn-
covid-19-misinformation/, where the journalists carefully check the
authenticity of news pieces and social media posts. Poynter provides an
international fact-checking network (IFCN) alliance collaborating with 100
fact-checking agencies such as FactCheck.org, PolitiFact, and Science Feedback
in more than 70 countries and 40 languages. The fact-checkers labeled the
articles and social media posts as false, partially false, misleading, or no
evidence according to the reports. We also obtain 813 low credibility sources
from Iffy+666https://iffy.news/iffy-plus/ dataset, which merges lists of sites
that regularly publish mis- and disinformation, as identified by major fact-
checking and journalism organizations.

###  3.2 Data in the Physical World (Offline Data)

COVID-19 Statistics: We adopt the COVID-19 data repository operated by the
Center for Systems Science and Engineering (CSSE) at Johns Hopkins University
[8]. This data repository obtained statewide, nationwide, and global streaming
COVID-19 statistics from reliable sources, including daily confirmed cases,
deaths, recoveries, hospitalizations, and vaccinations since January, 2020.

U.S. Census Bureau Data: We comprise the historical state-level and 5-digit
zip code-level census records from the U.S. Census
Bureau777https://data.census.gov/cedsci/ such as the population’s age, gender,
race, income, health insurance coverage, and employment status. Moreover, the
Household Pulse Survey(HPS), which focused on estimating how people’s lives
have been impacted by the COVID-19 pandemic, is available in our data
repository. It is worth noting that it consists of survey questions such as
“reasons for children ages 5-17 and adults not receiving or planning to
receive a COVID-19 vaccine”, which is valuable for COVID-19 vaccine hesitancy
research. Besides, the COVID-19 Vaccination Tracker from the Census Bureau is
also included for a high-level data visualization.

Government Responses: During the COVID-19 vaccine campaign, government
policies and responses play an essential role. We obtain the COVID-19-vaccine-
related federal and state policies from
Ballotpedia888https://ballotpedia.org/Documenting_America%27s_Path_to_Recovery.
It contains spatio-temporal information about COVID-19 vaccine authorization,
distribution, and legislation.

Local News: We use the GoogleNews [2] library to collect real-world events and
local news pieces from Google News, which aggregates news articles from more
than 4.5k sites in English and 50k news sources worldwide [5]. In this study,
after filtering out low-credibility sources, we acquired around 4.8k news
reports with keyword COVID-19 vaccine from Jan 2020 to Jan 2022.

##  4 Insight: What is in the Data?

###  4.1 Social Media Data

Table 2: Twitter data statistics Features Stances |  Pro-vaccine |  Anti-vaccine
---|---|---
Time period |  01.2020 - 01.2022 |  01.2020 - 01.2022
# Tweets |  1,495,991 |  335,229
# Retweets, Replies |  317,683/433,465 |  101,059/73,369
# Users |  480,327 |  87,405
# Tweets with URL |  233,673 |  48,635
# Tweets with image |  283,365 |  43,127
# Tweets with video/GIF |  310,699 |  55,550
Avg # tweets per user |  3.11 |  3.84
Avg # followings per user |  1,354 |  1,473
Avg # followers per user |  4,590 |  3,392
Avg # favourites per user |  1,374 |  1,506
Avg # lists per user |  36.2 |  28.1
Avg accounts age (years) |  7.91 |  6.55
Avg # replies per tweet |  1.08 |  1.42
Avg # retweets per tweet |  2.73 |  6.53
Avg # likes per tweet |  13.04 |  21.09
Avg # words per tweet |  25.67 |  23.67

(a) Provax Tweets Word Cloud

(b) Antivax Tweets Word Cloud

Figure 2: Word clouds of pro- and anti-vaccine tweets. The font size scaled to
the frequency.

The general statistics on our Twitter dataset are presented in Table 2. The
Twitter dataset contains 1,831,220 tweets, 366,276 ($\sim$20$\%$) of which
have both textual and visual (images or videos) information for multi-modal
studies. Figure 2 reveals the textual characteristics within pro- and anti-
COVID-19-vaccine tweets. In the pro-vaccine group, a sufficiently large number
of tweets encourage people to get the COVID-19 vaccines and wear masks to save
lives. Thus, we can observe that the words frequently appear in this group
include fully, vaccinate, wear, mask, save. The frequent words from the anti-
COVID-19-vaccine tweets, such as people, stop, freedom, government, passport,
show concern about excessive government control and people’s freedom during
the pandemic.

###  4.2 Fact-checking and Local News Data

Figure 3 shows the number of weekly fact-checking reports about COVID-19
vaccines from Poynter and the distribution of active fact-checking websites.
We can observe that the number of COVID-19 fact-checking articles published is
significantly increased from December 2020 to March 2020, highlighted with a
pink rectangle. Note that the emergency use authorization (EUA) approved the
first COVID-19 vaccine from Pfizer and BioNTech on December 11, 2020, which
potentially led to a massive spread of misinformation about vaccine safety,
effectiveness, and adverse reactions. Furthermore, this dataset contains fact-
checking websites from various countries, which provides an potential for
multilingual misinformation studies on COVID-19 vaccines. Figure 4 presents
the word cloud and top news publishers of the local news dataset. We can
observe that news authors’ frequently used topics and vocabularies are more
neutral and objective than those of Twitter users. A possible explanation is
that the news articles collected are from multiple sources, including cross-
domain for-profit and non-profit news outlets, multidisciplinary science
journals, and international health magazines.

(a) Timeline

(b) Top Fact-checking Websites

Figure 3: Number of weekly fact-checking reports and top publishers.

(a) Word Cloud of News Articles

(b) Top News Publisher

Figure 4: Insights of local news articles.

###  4.3 Offline Statistics

In the repository, CSV files are well-documented with sufficient offline
statistics. The JHU COVID-19 dashboard[8] provides easy-access, interactive,
and comprehensive data visualization of spatio-temporal COVID-19 statistics.
Note that COVID-19 statistics and Census data have shown to be correlated with
vaccine hesitancy in a quantitative study[14]. Similar findings from surveys
can be good validations. For example, Kaiser Family Foundation (KFF) compares
public attitudes towards COVID-19 vaccines among different demographic and
socioeconomic statuses through a monthly survey. We can observe that: (i)
political lean is correlated with vaccine attitudes $-$ Republicans tend to
refuse the COVID-19 vaccine more than Democrats; (ii) The vaccine attitudes
are similar among different races across time; (iii) People in rural areas are
more likely to be anti-vaxxers than those in suburban and urban areas; (iv)
People with a college degree are more willing to vaccinate than those with
less than a college degree; and (v) The percentage of pro-vaxxers is
increasing while anti-vaxxers’ remains stable and low over time among all
demographic and socioeconomic groups.

##  5 Online and Offline Data Connection and Application

The COVID-19 vaccine campaign gave rise to a surplus of online and offline
data expressing people’s views and different policies imposed by the
government to promote the vaccine. These two kinds of data share a strong
relationship. For example, offline COVID-19 statistics can alter the vaccine
stance of the online crowd [13]. Thus, online and offline data can complement
each other to optimize the information capture. However, connecting online and
offline data is challenging because they reside in different sources,
modalities, and dimensions. To this end, we first discuss how to connect
COVID-19 online and offline data in CoVaxNet. Furthermore, we conduct COVID-19
vaccine stance detection to (i) verify the quality and compatibility of our
datasets for various downstream tasks, including stance detection; (ii)
demonstrate the necessity for employing online-offline data; and (iii)
encourage researchers to develop state-of-the-art models and incorporate
various online and offline data provided by CoVaxNet to facilitate the
COVID-19 vaccine hesitancy study.

###  5.1 Data Connection

An online-offline dataset contains (i) an online dataset and (ii) an offline
dataset. Thus, we denote the online-offline dataset $\mathcal{D}_{online-
offline}$ as

$\mathcal{D}_{online-
offline}=(\mathcal{D}_{online},\,\mathcal{D}_{offline}),\vspace{-0.2cm}$ (1)

where
$\mathcal{D}_{online}=\\{{d}_{online}^{1},\,{d}_{online}^{2},\,...,\,{d}_{online}^{n},\,...,\,{d}_{online}^{N}\\}$
denotes the set of $N$ online data and
$\mathcal{D}_{offline}=\\{{d}_{offline}^{1},\,{d}_{offline}^{2},\,...,\,{d}_{offline}^{m},\,...,\,{d}_{offline}^{M}\\}$
is the set of $M$ offline data. Then an online-offline data connection
$\mathcal{C}$ which establishes a link between the online and offline dataset,
presented by implicit or explicit features, can be extracted. In this work, we
connect online social media dataset $\mathcal{D}_{online}$ with offline
COVID-19 statistics dataset $\mathcal{D}_{offline}$ through a geolocation-
based online-offline connection $\mathcal{C}_{geo}$. For example, given two
online social media posts $\mathcal{P}_{1}$ and $\mathcal{P}_{2}$ from
$\mathcal{D}_{online}$ and two offline COVID-19 statistics $\mathcal{S}_{1}$
and $\mathcal{S}_{2}$ from $\mathcal{D}_{offline}$, and a connection
$\mathcal{C}_{geo}$ such that

$\displaystyle\mathcal{C}_{geo}(\mathcal{P}_{1})=\mathcal{C}_{geo}(\mathcal{S}_{1})=\\{\textit{Tempe,
AZ, USA}\\}$
$\displaystyle\mathcal{C}_{geo}(\mathcal{P}_{2})=\mathcal{C}_{geo}(\mathcal{S}_{2})=\\{\textit{Pittsburgh,
PA, USA}\\}.\vspace{-0.2cm}$

Therefore, the online-offline dataset $\mathcal{D}_{online-offline}$, can be
formulated as

$\begin{split}\mathcal{D}_{online-
offline}&=(\mathcal{D}_{online},\,\mathcal{D}_{offline})\\\
&=\\{<\mathcal{P}_{1},\mathcal{S}_{1}>,\,<\mathcal{P}_{2},\mathcal{S}_{2}>,\,...\\},\end{split}\vspace{-0.2cm}$
(2)

where the tuple $<\mathcal{P}_{1},\mathcal{S}_{1}>$ indicates the post
$\mathcal{P}_{1}$ and the statistics $\mathcal{S}_{1}$ are linked to each
other through the geolocation indicator $\mathcal{C}_{geo}$.

###  5.2 COVID-19 Vaccine Stance Detection

We adopt $\mathcal{D}_{online}$, $\mathcal{D}_{offline}$, and
$\mathcal{D}_{online-offline}$ in the last subsection to perform vaccine
stance detection. For data preprocessing, we encode the tabular COVID-19
statistics into text representations so as to make them consistent with the
modality of tweets. First, we compute the weekly COVID-19 severity of each
city

$Severity=\sum_{i=1}^{7}(\mathcal{G}_{i}^{c}+\mathcal{G}_{i}^{d}+\mathcal{G}_{i}^{h}-\mathcal{G}_{i}^{r}-\mathcal{G}_{i}^{v}),\vspace{-0.2cm}$
(3)

where $i$ is the number of days before a tweet is posted.
$\mathcal{G}_{i}^{c}$, $\mathcal{G}_{i}^{d}$, $\mathcal{G}_{i}^{h}$,
$\mathcal{G}_{i}^{r}$, and $\mathcal{G}_{i}^{v}$ denote the daily growth rate
of cases, deaths, hospitalizations, recoveries, and vaccinations. Second, for
a given list of cities, we categorize the $<$severity level$>$ of each city as
low, medium, or high by sorting the Severity values in an ascending order and
then splitting the list based on the proportion of $0.33:0.33:0.33$. Third, we
encode the $<$severity level$>$ at the end of each tweet. For example, a
encoded tweet will be “mRNA is made from Human DNA. It’s designed to attach to
your DNA. The COVID-19 severity of my city is <severity level>.” In addition,
we encode the user demographics, including age, gender, and race, from the
user profile into the tweets. For instance, "I am a $<$age$>$ years old
$<$race$>$ $<$gender$>$." Now we have three types of data: (i) online data
(tweets); (ii) offline data (encoded COVID-19 statistics); and (iii) online-
offline data (tweets $+$ encoded COVID-19 statistics and user demographics).

We randomly selected 70% data as the training set and 30% as the testing set,
as well as a balanced distribution ($1:1$) between pro- and anti-vaccine
tweets. We deploy four well-established baselines including Support Vector
Machines (SVM), Naive Bayes (NB), Logistic Regression (LR), and BERT in our
experiments. The experimental results are listed in Table 3. We observe that
four baselines have reasonably good performance results in terms of F1-score
and accuracy. This validates the quality of labels and collected features in
CoVaxNet. Moreover, we can find that models trained on online-offline data
outperform models trained on either online or offline data. This demonstrates
that online-offline data can provide more critical information than online or
offline data, which can be learned by the baseline models. However, the
experimental results of these baselines are still limited and shows great
potential for improvements.

Table 3: COVID-19 vaccine stance detection performance on CoVaxNet Method | SVM | NB | LR | BERT
---|---|---|---|---
Metric | Accuracy | F1 | Accuracy | F1 | Accuracy | F1 | Accuracy | F1
Online | 0.736 | 0.763 | 0.758 | 0.751 | 0.776 | 0.777 | 0.784 | 0.788
Offline | 0.630 | 0.602 | 0.658 | 0.667 | 0.674 | 0.676 | 0.675 | 0.680
Online-Offline | 0.759 | 0.769 | 0.782 | 0.777 | 0.789 | 0.788 | 0.796 | 0.793

##  6 Conclusion

To combat the COVID-19 pandemic, we present a multi-source, multi-modal, and
multi-feature online-offline repository for COVID-19 vaccine hesitancy
research. We introduce our strategy for constructing a high-quality online-
offline data repository. We share sufficient insights, which detail textual,
visual, spatio-temporal, and network features in the data. Furthermore, we
propose a geolocation-based online-offline data generation approach. We
believe CoVaxNet would facilitate the research of analyzing the impact of
various online and offline data on COVID-19 vaccine hesitancy and help
policymakers prioritize resource allocation in the ongoing and future
pandemics. There are several promising future directions. First, we will
improve our labeling strategy to reduce noise and provide ground truth of
“vaccine hesitancy” in our Twitter dataset. Second, we intend to extend the
CoVaxNet to include more relevant data from other reliable sources. For
example, we will collect data from countries with similar and dissimilar
COVID-19 vaccine roll-out plans. This will help researchers study the
structural inequalities and the correlation among all information regarding
vaccine uptake during the pandemic.

## References

  * [1] Geopy, https://geopy.readthedocs.io/en/stable/, accessed June 22
  * [2] Googlenews, https://github.com/Iceloof/GoogleNews, accessed June 22
  * [3] Nameparser, https://nameparser.readthedocs.io/en/latest/, accessed June 22
  * [4] Bradski, G., et al.: Learning OpenCV: Computer vision with the OpenCV library
  * [5] Das, A.S., et al.: Google news personalization: scalable online collaborative filtering. In: WWW (2007)
  * [6] DeVerna, M.R., et al.: Covaxxy: A collection of english-language twitter posts about covid-19 vaccines. In: ICWSM (2021)
  * [7] Di Giovanni, M., et al.: Vaccineu: Covid-19 vaccine conversations on twitter in french, german and italian. In: ICWSM (2022)
  * [8] Dong, E., et al.: An interactive web-based dashboard to track covid-19 in real time. Lancet Infect. Dis. (2020)
  * [9] Guo, J., et al.: Vaccinations, mobility and covid-19 transmission. IJERPH (2021)
  * [10] Hale, T., et al.: A global panel database of pandemic policies (oxford covid-19 government response tracker). Nature human behaviour (2021)
  * [11] Hussain, A., et al.: Artificial intelligence–enabled analysis of public attitudes on facebook and twitter toward covid-19 vaccines in the united kingdom and the united states: Observational study. JMIR (2021)
  * [12] Kim, J.H., et al.: Looking beyond covid-19 vaccine phase 3 trials. Nat. Med. (2021)
  * [13] Lyu, H., et al.: Social media study of public opinions on potential covid-19 vaccines: informing dissent, disparities, and dissemination. Intelligent medicine (2021)
  * [14] Lyu, H., et al.: Misinformation versus facts: Understanding the influence of news regarding covid-19 vaccines on vaccine uptake. Health Data Science (2022)
  * [15] MacDonald, N.E., et al.: Vaccine hesitancy: Definition, scope and determinants. Vaccine (2015)
  * [16] McRee, A.L., et al.: Hpv vaccine hesitancy: findings from a statewide survey of health care providers. JPHC (2014)
  * [17] Mesch, G.S., et al.: Social and political determinants of vaccine hesitancy: Lessons learned from the h1n1 pandemic of 2009-2010. AJIC (2015)
  * [18] Mollalo, A., et al.: Spatial modeling of covid-19 vaccine hesitancy in the united states. IJERPH (2021)
  * [19] Muric, G., et al.: Covid-19 vaccine hesitancy on social media: Building a public twitter data set of antivaccine content, vaccine misinformation, and conspiracies. JPHS (2021)
  * [20] Oluyomi, A.O., et al.: Covid-19 community incidence and associated neighborhood-level characteristics in houston, texas, usa. IJERPH (2021)
  * [21] Sood, G., et al.: Predicting race and ethnicity from the sequence of characters in a name. arXiv preprint arXiv:1805.02109 (2018)
  * [22] Taigman, Y., et al.: Deepface: Closing the gap to human-level performance in face verification. In: CVPR (2014)