# Modeling of Ex-Situ Dissolution for Geologic Sequestration of Carbon Dioxide
in Aquifers

Federico Cao  Department of Applied Mathematics, University of Waterloo,
Waterloo, Ontario, Canada  Dmitry Eskin  Skolkovo Institute of Science &
Technology, Moscow, Russia  Yuri Leonenko  Department of Earth and
Environmental Sciences, University of Waterloo, Waterloo, Ontario, Canada
<EMAIL_ADDRESS>
###### Abstract

Underground carbon dioxide ($CO_{2}$) sequestration is considered to be one of
the main methods to mitigate greenhouse gas (GHG) emissions. In this
technology, pure $CO_{2}$ is injected into an underground geological formation
and since it is less dense than residual fluids, there is always a risk of
leakage to the surface. To increase security of underground $CO_{2}$ disposal,
ex-situ dissolution can be implemented. When $CO_{2}$ is dissolved in brine
before injection, it significantly reduces the risks of leakage. In this
approach, pure $CO_{2}$ is dissolved on the surface before injection. Surface
dissolution could be achieved in a pipeline operating under the pressure of a
target aquifer into which the $CO_{2}$ is injected. In a pipeline, $CO_{2}$
droplets are dissolved being dispersed in a brine turbulent flow. In this
paper, a comprehensive model of droplet dissolution along a pipeline is
presented. The model accounts for droplet breakup and coalescence processes
and is validated against available experimental data.

###### keywords:

ex-situ dissolution; geologic sequestration; modeling

Department of Geography and Environmental Management, University of Waterloo,
Waterloo, Ontario, Canada CO_2, CCS, GHG, IPCC

##  1 Introduction

As rapidly developing economies require higher energy consumption, it is clear
that major greenhouse gas (GHG) emitting sources cannot be avoided in future
decades. Currently, there are serious limitations to alternative/sustainable
energy sources as they are still cost-prohibitive for many industries and in
developing countries. As a result, fossil fuel consumption will continue to be
the main source in the near future. Annual global carbon emissions from fossil
fuels have increased to nearly 10 billion metric tons in 2014 (Boden et al.,
2017). With $CO_{2}$ as the most common GHG and responsible for 65% of
anthropogenic global warming, it is crucial to determine feasible mitigation
measures. In particular, the Intergovernmental Panel of Climate Change (IPCC)
reports that carbon capture and storage (CCS) methods can be an effective
solution in significantly lowering the amount of $CO_{2}$ in the atmosphere
(Metz et al., 2005).

Various CCS technologies exist; however, $CO_{2}$ sequestration in sedimentary
basins is of particular interest among many. This form of sequestration relies
upon depleted oil and gas reservoirs, (Herzog, 2001 and Jenkins et al., 2012)
unmineable coal bed reservoirs (Shi and Durucan, 2005) and deep aquifers
(Celia et al., 2015) where the saline water (brine) is not suitable for
agricultural or consumption purposes. Another sequestration option is ocean
storage, where $CO_{2}$ would be injected into the ocean at depths of over one
thousand meters (Haugan and Joos, 2004). Among the above options, deep saline
aquifers represent the largest long term potential for CCS (International
Energy Agency (2008), IEA). An IPCC special report (Metz et al., 2005) has
suggested that deep saline formations have a storage capacity of around 2000
Gigatons (Gt) of $CO_{2}$. It is approximately two orders of magnitude higher
than the total annual worldwide emissions amount, making saline aquifers the
most viable disposal option. Although it was recognized that deep saline
aquifers offer very large potential storage capacity, significant
uncertainties remain regarding storage security. $CO_{2}$ injected into a
saline aquifer is less dense than the resident brine and, driven by buoyancy,
will flow horizontally, spreading under the cap-rock which should confine
$CO_{2}$ for thousands of years until it is fully dissolved. Cap-rocks of
aquitards have not been proven to hold buoyant $CO_{2}$ for geologic time
scales as in the case for cap-rocks that have confined buoyant oil and gas
(van der Meer, 1993; Lindeberg, 1997). It also may flow upward, leaking
through any high permeability zones such as natural fractures or artificial
penetrations such as abandoned wells. Therefore, approaches which allow an
increase of storage security are of great importance for developing and
implementing CCS technologies. One of the ways to increase storage security is
to accelerate dissolution of $CO_{2}$ in a formation brine and several methods
have been reported in the literature. After $CO_{2}$ is injected into a
reservoir it starts to dissolve naturally (natural convection). Both the onset
of convection and the dissolution rate are strongly dependent on reservoir
parameters, and an overview of these phenomena is studied in detail in the
work of Emami-Meybodi et al. (2015). The time scale for complete dissolution,
in this case, could reach $10^{3}-10^{5}$ years. To speed up the dissolution
process a variety of engineering options has been proposed. For example,
Hassanzadeh et al. (2009) and Keith et al. (2005) suggested a method for
enhancing $CO_{2}$ dissolution in saline aquifers by injecting brine on top of
the injected $CO_{2}$. This technique significantly increases the dissolution
rate, and half of the injected $CO_{2}$ is dissolved within 200 years. To
further improve the dissolution process, employment of mixing devices, such as
a static mixer (Zirrahi et al., 2013a) or a downhole mixer (Zirrahi et al.,
2013b) has been studied. Shafaei et al. (2012) investigated using an injection
well of a special design, where brine and $CO_{2}$ are co-injected. Those
devices have the potential to significantly increase the dissolution rate. In
the current study, we investigate a possibility of nearly $100\%$ dissolution
before it is injected underground.

The idea was first proposed in our previous study (Leonenko and Keith, 2008),
where we adopted the view that the only relevant risk of leakage arises from
mobile free-phase $CO_{2}$, which is not immobilized by residual or chemical
trapping or dissolution. Therefore storage security mainly depends on two
factors: (a) the likelihood that free-phase $CO_{2}$ will leak out of the
reservoir and (b) the rate at which free-phase $CO_{2}$ is immobilized by one
of the trapping mechanisms. Storage security then can be increased either by
reducing the probability of leakage or by increasing the rate at which
$CO_{2}$ is immobilized by residual gas trapping, dissolution in reservoir
fluids, or geochemical reactions. In the same study we proposed some options
to reduce the time scale of free phase of $CO_{2}$: in-situ and ex-situ
dissolution. The latter could be achieved within a surface pipeline where two
phase $CO_{2}$-brine mixture flow takes place. The generation of $CO_{2}$
droplets, which are sufficiently small to achieve rapid dissolution, occurs in
a turbulent pipe flow. We present a diagram of the ex-situ dissolution
procedure in Fig. 1.

Figure 1: A representation of ex-situ dissolution. Captured $CO_{2}$ is mixed
with brine along a turbulent pipeline flow. The $CO_{2}$ saturated brine is
then injected underground into saline aquifers.

In our former studies, (Zendehboudi et al., 2011; Cholewinski and Leonenko,
2013) mass transfer from $CO_{2}$ droplets into brine during co-current
pipeline flow was modeled to investigate the effectiveness of the proposed
method. The models, however, did not include droplet breakup and coalescence.
In a subsequent study (Zendehboudi et al., 2013), a very simplistic model of
droplet breakup was employed where coalescence was entirely ignored. Let us
briefly discuss the flaws of the model of Zendehboudi et al. (2013). Other
mentioned prior models are based on an even larger number of assumptions.
First, the model of Zendehboudi et al. (2013) describes mean droplet diameter
evolution along a pipe, whereas a real dispersed system of droplets is
strongly polydispersed. Second, this model assumes a constant turbulence
energy dissipation rate over a pipe cross-section, whereas it is maximum at
the pipe wall and very strongly decreases toward the pipe center. The energy
dissipation rate is the predominant driver of all the dispersion processes
considered. Third, the assumption made by Zendehboudi et al. (2013) that mass
transfer between a droplet and an ambient fluid can be calculated through
turbulent diffusion flux is incorrect. Levich (1962) showed that the diffusion
flux is calculated through the mean relative droplet-fluid velocity, caused by
turbulent fluctuations. The dispersion-dissolution model proposed in the
present work is free of the aforementioned drawbacks because it takes into
account droplet polydispersity, non-uniformity of flow parameters across a
pipe cross-section and employs a validated empirical correlation to calculate
the mass transfer coefficient for a droplet moving in a turbulent flow.

Thus, the model of Zendehboudi et al. (2013) could be considered only as an
engineering estimation of the dissolution process, whereas the model we
developed represents a mathematical description of the complex dissolution
process based on first principles.

In the current paper we perform comprehensive modeling of ex-situ dissolution
by incorporating all three phenomena which take place in a pipeline: droplet
coalescence, breakup, and dissolution. Dispersion of droplets is modeled by an
advection-diffusion population balance equation. The numerical results
obtained are validated versus the available experimental data.

##  2 Modeling

To simulate dispersion of $CO_{2}$ droplets in the pipeline, let us formulate
the major assumptions:

  1. 1.

Droplet dispersion over a pipe cross-section is intense enough to neglect
gravity induced droplet stratification.

  2. 2.

Flow is steady-state.

Assumption 1 appears to be reasonable taking into account that the density of
liquid carbon dioxide for typical aquifer conditions ($700-800$ kg/m3) is not
significantly different from that of water. Furthermore, $CO_{2}$ droplets are
rapidly fragmented in a turbulent pipe flow down to small sizes. Due to
turbulent diffusion, such droplets are rather uniformly distributed across a
pipe even if it is horizontal. Also, dissolution significantly contributes to
a droplet size reduction, additionally supporting the validity of the latter
assumption. We would like to emphasize that we calculated a steady-state
hydrotransport of particles in a water pipe flow by the model of Eskin (2012)
assuming that the particle density is equal to that of liquid $CO_{2}$. For
droplets of size $150$ $\mu$m transported at a holdup $0.1$ in a pipe of the
diameter $0.15$ m with the mean velocity $3$ m/s, a droplet volume
concentration variation across a pipe turned out to be only about $20$%. This
variation is relatively small, whereas even initially large droplets in a
turbulent pipe flow are rapidly dispersed (Eskin et al., 2017a) reaching the
aforementioned size of $150$ $\mu$m; therefore, assumption 1 is acceptable.
Assumption 2 is straightforward because we consider steady-state operational
conditions with a constant flow rate.
Let us formulate the advection-diffusion-population balance equation for a
steady-state pipe flow in cylindrical coordinates for a dispersed phase
discretized by $M$ size fractions as follows (Eskin et al., 2017a):

$\displaystyle u(r)\frac{\partial{N_{i}}}{\partial{x}}$
$\displaystyle=\frac{1}{r}\frac{\partial}{\partial{r}}\left(rD_{i}(r)\frac{\partial{N_{i}}}{\partial{r}}\right)+u(r)\left(\frac{\partial{N_{i}}}{\partial{x}}\right)_{PB},\hskip
14.22636pti=1,...,M.$ (1)

Here, $i$ is the size fraction number, $r$ is the radial coordinate, $u(r)$ is
the flow velocity, $N_{i}$ is the number concentration of droplets of the
$i$-th size fraction in a computational cell, $D_{i}(r)$ is the turbulent
diffusivity of a droplet of the $i$-th size fraction, and
$\left(\dfrac{\partial{N_{i}}}{\partial{x}}\right)_{PB}$ is the number
concentration derivative accounting for coalescence, breakup and dissolution
for the $i$-th size fraction. The boundary conditions for Eq. (1) are
formulated as

  1. 1.

The volume flux through the pipe wall is zero:

$q_{i}(R)=-D_{i}\dfrac{\partial{N_{i}}}{\partial{r}}\Bigr{|}_{r=R}=0.$ (2)

  2. 2.

The dispersed phase concentration gradient at the pipe axis is zero:

$\dfrac{\partial{N_{i}}}{\partial{r}}\Bigr{|}_{r=0}=0.$ (3)

  3. 3.

The droplet size distribution at the initial pipe cross-section is:

$N_{i}(0,r)=\Psi_{i}(r)$ (4)

where $\Psi_{i}(r)$ is a some known function.

Eq. (1) can be rewritten in a one-dimensional time-dependent form as follows:

$\dfrac{u(r)}{U}\dfrac{\partial{N_{i}}}{\partial{t}}=\dfrac{1}{r}\dfrac{\partial}{\partial{r}}r\left(D_{i}(r)\dfrac{\partial{N_{i}}}{\partial{r}}\right)+\dfrac{u(r)}{U}\left(\dfrac{\partial
N_{i}}{\partial t}\right)_{PB}$ (5)

where $U$ is the mean pipe flow velocity and $dt=dx/U$.

Further, we need to describe physical processes accompanying droplet size
evolution in a pipe. Fortunately, dispersed pipeline flows have been rather
intensely studied in the past; therefore, for our modelling we will use an
engineering approach based on mainly validated ideas. The following sections
describe models employed to calculate the terms of Eq. (5).

###  2.1 Flow velocity distribution across a pipe

The steady-state velocity distribution across a pipe $u(r)$ is assumed to be
consisting of the two regions: 1) the viscous layer in the wall vicinity,
characterized by a linear velocity distribution; 2) the turbulent boundary
layer, extended from the viscous layer to the pipe center and characterized by
a logarithmic velocity distribution. In dimensionless coordinates, this
velocity distribution is written as (Schlichting and Gersten, 2000):

$\displaystyle u^{+}$ $\displaystyle=y^{+},\hskip 14.22636pty^{+}\leq 11.6$
(6) $\displaystyle u^{+}$ $\displaystyle=2.5\ln y^{+}+5.5,\hskip
14.22636pty^{+}>11.6$ (7)

where $u^{+}=u/u_{*}$ is the dimensionless flow velocity,
$u_{*}=(\tau_{w}/\rho_{f})^{0.5}$ is the friction velocity,
$y^{+}=u_{*}y/\nu_{f}$ is the dimensionless coordinate, $y=R-r$ is the
distance from the wall, $\nu_{f}$ is the fluid kinematic viscosity, and
$\tau_{w}$ is the wall shear stress.

###  2.2 Turbulent diffusivity of droplets

The droplet turbulent diffusivity can be determined as:

$\displaystyle D_{p}=\dfrac{D_{t}}{Sc_{pt}}$ (8)

where $D_{t}=v_{t}/Sc_{t}$ is the turbulent diffusivity, $v_{t}$ is the eddy
diffusivity, $Sc_{t}$ is the Schmidt number for a fluid particle in a
turbulent flow, and $Sc_{pt}$ is the turbulent Schmidt number for a particle
(droplet).
In the present work, to determine the droplet turbulent diffusivity
distribution across a pipe, we will assume that the droplet turbulent
diffusivity is equal to the eddy diffusivity:

$D_{i}=v_{t}.$ (9)

Since for the dispersion system considered in the present work droplet sizes
as well as a difference between densities of a fluid and a dispersed phase are
relatively small, Eq. (9) is justified.
The dimensionless eddy diffusivity distribution across a pipe can be
calculated using the empirical equations suggested by Johansen (1991):

$\displaystyle v_{t}^{+}$
$\displaystyle\approx\dfrac{v_{t}}{\nu_{f}}=\left(\dfrac{y^{+}}{11.15}\right)^{3}\hskip
14.22636pt\text{for}\hskip 14.22636pty^{+}\leq 3$ (10) $\displaystyle
v_{t}^{+}$
$\displaystyle\approx\left(\dfrac{y^{+}}{11.4}\right)^{2}-0.049774\hskip
14.22636pt\text{for}\hskip 14.22636pt3<y^{+}\leq 52.108$ (11) $\displaystyle
v_{t}^{+}$ $\displaystyle\approx\kappa y^{+}\hskip 14.22636pt\text{for}\hskip
14.22636pty^{+}>52.108$ (12)

where $\kappa=0.406$ is the von Karman constant.

The derivative, expressing the population balance term in Eq. (5) can be
represented as a sum of the derivatives determining the contributions of
breakup, coalescence and dissolution respectively:

$\displaystyle\left(\dfrac{\partial N_{i}}{\partial
t}\right)_{PB}=\left(\dfrac{\partial N_{i}}{\partial
t}\right)_{break}+\left(\dfrac{\partial N_{i}}{\partial
t}\right)_{coal}+\left(\dfrac{\partial N_{i}}{\partial t}\right)_{diss}.$ (13)

For numerical calculation of both breakup and coalescence terms, we employed
the Fixed Pivot Approach of Kumar and Ramkrishna (1996).

###  2.3 Breakup term

For droplet breakup modeling, we employ a usual binary breakup assumption:
fragmentation of a mother droplet leads to formation of two daughter droplets.
The derivative associated with droplet breakup is calculated as (Kumar and
Ramkrishna, 1996):

$\displaystyle\left(\dfrac{\partial N_{i}}{\partial
t}\right)_{break}=\sum_{k=1}^{M}n_{i,k}G_{k}N_{k}-G_{i}N_{i}$ (14)

where $G_{i}$ is the breakup rate of a droplet belonging to the $i$-th size
fraction (a model will be given further). The function $n_{i,k}$ is calculated
as follows:

$\displaystyle
n_{i,k}=\int_{x_{i}}^{x_{i+1}}\dfrac{x_{i+1}-v}{x_{i+1}-x_{i}}\beta(v,x_{k})dv+\int_{x_{i-1}}^{x_{i}}\dfrac{v-x_{i-1}}{x_{i}-x_{i-1}}\beta(v,x_{k})dv$
(15)

where $\beta(v,x_{k})$ is the droplet breakup density function characterizing
the probability of formation of a droplet of the $k$-th size fraction of the
volume $x_{k}$ at the breakup of a droplet of volume $v$. Both first and
second integrals for $i=1$ and $i=k$ are zero respectively.

To use (14), the equations of both the droplet breakup rate and breakup
density function should be specified. In the present research, we will use the
functions employed by Eskin et al. (2017a and 2017b) in their work on modeling
droplet dispersion in a pipe. The equation for breakup rate at a low droplet
concentration is as follows:

$\displaystyle G(d)=K\dfrac{(\epsilon
d)^{1/3}}{d}\left[erfc(\Phi^{1/2})+\dfrac{2}{\pi^{1/2}}\Phi^{1/2}\exp(-\Phi)\right]$
(16)

where $erfc$ is the complimentary error function, $K$ is the model parameter,
$\epsilon$ is the energy dissipation rate per unit mass,
$\Phi=\dfrac{3}{2}\dfrac{We_{cr}}{We}$ is the dimensionless parameter,
$We=\dfrac{2\rho_{f}(\epsilon d)^{2/3}d}{\gamma}$ is the Weber number for a
droplet, $We_{cr}$ is the critical Weber number (experimental parameter).
Eskin et al. (2017b) identified $We_{cr}$ from the Couette device experiments
as $We_{cr}=0.5$. This parameter determines the steady-state droplet size
distribution. Due to a very rapid dispersion process in a Couette device, the
parameter $K$, determining the rate of a size distribution change, was not
identified accurately and was assumed to be $K=1$ (Eskin et al., 2017b). In
the present work, we will employ the critical Weber number recommended by
Eskin et al. (2017b) and allow the parameter $K$ to be varied to fit the
$CO_{2}$ droplet dissolution experimental data. We will also employ the same
breakup density function as Eskin et al. (2017a) used in their research on
droplet dispersion in a pipe:

$\displaystyle\beta(f_{bv})=12f_{bv}(1-f_{bv})$ (17)

where $f_{bv}=v/x$ is the breakup fraction, $x$ and $\beta$ are the mother and
smaller droplet volumes respectively. According to Eskin et al. (2017b), this
function weakly affects particle size distributions calculated by solving the
population balance equation (Eq. (14)).

The energy dissipation rate for the employed velocity profile (Eq. (6), (7))
is calculated as the specific power spent on friction between concentric fluid
layers (Eskin et al., 2017a):

$\displaystyle\epsilon(\tilde{r})=\xi\dfrac{\tilde{r}^{3/2}}{1-\tilde{r}}$
(18)

where $\tilde{r}=r/R$ is the dimensionless radius, $\xi=(-0.5\nabla
p/\rho_{f})^{1.5}\sqrt{R}/\kappa$ is the dimensional complex, $\nabla p$ is
the pressure gradient, $\kappa$ is the von Karman constant.
The pressure gradient in a pipe flow is calculated as:

$\displaystyle\nabla p=-2\rho_{f}fU^{2}/D.$ (19)

In the present work, we will calculate the Fanning friction factor $f$ by the
Blausius equation, valid for hydraulically smooth pipes, as (e.g., Bird et
al., 2002):

$\displaystyle f=\dfrac{0.079}{\text{Re}^{0.25}}$ (20)

where Re $=UD/\nu_{f}$ is the pipe Reynolds number.

###  2.4 Coalescence term

The derivative we used to account for coalescence in Eq. (13) is calculated as
(Kumar and Ramkrishna, 1996):

$\displaystyle\left(\dfrac{\partial N_{i}}{\partial
t}\right)_{coal}=\sum_{\begin{subarray}{c}j,k\\\ x_{i-1}\leq x_{j}+x_{k}\leq
x_{i+1}\end{subarray}}^{j\geq k}(1-0.5\delta_{jk})\eta
Q_{j,k}N_{j}N_{k}-N_{i}\sum\limits_{k=1}^{M}Q_{i,k}N_{k}$ (21)

where $\delta_{jk}$ is the Kronecker delta function and $Q_{j,k}$ is the
coalescence rate of droplets belonging to the $j$ and $k$ size fractions.
The variable $\eta$ is calculated by the following equations:

$\displaystyle\eta=$ $\displaystyle\dfrac{x_{i+1}-v}{x_{i+1}-x_{i}},$
$x_{i}\leq v\leq x_{i+1}$ (22) $\displaystyle\eta=$
$\displaystyle\dfrac{v-x_{i-1}}{x_{i}-x_{i-1}},$ $x_{i-1}\leq v\leq x_{i}$
(23)

There are many models for calculation of the coalescence rate, which are
available in the literature (Liao and Lucas, 2010). However, literature
analysis (Liao and Lucas, 2010) shows that different models predict
significantly different results. In the present work, we employed a
coalescence model (kernel) suggested by Coulaloglou and Tavlarides (1977):

$Q_{j,k}=\alpha(d_{j},d_{k})\omega(d_{j},d_{k})$ (24)

where $\alpha(d_{j},d_{k})$ is the coalescence efficiency, and
$\omega(d_{j},d_{k})$ is the collision frequency of randomly fluctuating
spheres that is calculated by a well-known equation as (e.g. Coulaloglou and
Tavlarides, 1977):

$\displaystyle\omega(d_{j},d_{k})=C_{1}(d_{j}+d_{k})^{2}\epsilon^{1/3}(d_{j}^{2/3}+d_{k}^{2/3})^{1/2}$
(25)

where $C_{1}$ is the model parameter.
Coulaloglou and Tavlarides (1977) suggested the following equation for the
coalescence efficiency:

$\displaystyle\alpha(d_{j},d_{k})=\exp\left(-C_{2}\dfrac{\mu_{f}\rho_{f}\epsilon}{\gamma^{2}}\left(\dfrac{d_{j}d_{k}}{d_{j}+d_{k}}\right)^{4}\right)$
(26)

where $C_{2}$ is the model parameter.
The coalescence model parameters, $C_{1}$ and $C_{2}$, identified under
specific conditions can be found in the literature (e.g. Laakkonen et al.,
2006). Note, the parameter $C_{2}$ in Eq. (26) is dimensional.
Since unique values of the model coalescence parameters do not exist, they
usually need to be tuned to fit experimental data. In our further analysis we
assumed $C_{1}=1$ that is close to the value 0.88 suggested by Laakkonen et
al. (2006), whereas the parameter $C_{2}$ was used as a tunable parameter to
fit the experimental data. The parameter $C_{2}$ was chosen to be tunable
because experiments, employed for model validation, were conducted for a
brine-liquid $CO_{2}$ system, where coalescence rate was expected to be
significantly lower than in the air-water system studied by Laakkonen et al.
(2006).

###  2.5 Dissolution term

The droplet-fluid mass transfer (dissolution) process is a key phenomenon
defining behavior of soluble droplets in a turbulent flow.
We calculated the derivative, expressing the dissolution term in the
population balance equation as follows:

$\left(\dfrac{\partial N_{i}}{\partial
t}\right)_{diss}=\dfrac{N_{i+1}}{x_{i+1}-x_{i}}\bigg{|}\dfrac{dx_{i+1}}{dt}\bigg{|}_{diss}-\dfrac{N_{i}}{x_{i}-x_{i-1}}\bigg{|}\dfrac{dx_{i}}{dt}\bigg{|}_{diss}$
(27)

where $(d{x_{i}}/dt)_{diss}$ is the rate of size change of a droplet of the
$i$-th size fraction only due to dissolution.
This equation is obtained from mass balance formulated for a droplet of the
$i$-th size fraction and its neighbouring size fractions during dissolution.
The derivation of Eq. (27) is given in the Appendix. The following limitation
is imposed on application of Eq. (27): a volume reduction of a droplet of the
size fraction $i$ during a single time step should be smaller than a
difference of volumes of droplets of the fractions $i$ and $i-1$ respectively.

The dissolution rate for the $i$-th size fraction droplet is calculated as
follows

$\displaystyle\left(\dfrac{dx_{i}}{dt}\right)_{diss}=-\dfrac{k\pi^{1/3}(6x_{i})^{2/3}}{\rho_{d}}(C_{s}-C_{\infty})$
(28)

where $C_{s}$ is the saturation concentration of carbon dioxide in a bulk
water fluid, $C_{\infty}$ is the concentration of dissolved carbon dioxide in
water, $k$ is the mass transfer coefficient, $\rho_{d}$ is the droplet
density.

The mass transfer coefficient is determined as:

$\displaystyle k=\dfrac{ShD_{CO_{2}}}{d}$ (29)

where $D_{CO_{2}}$ is the molecular diffusivity of carbon dioxide in water and
$Sh$ is the Sherwood number.
The Sherwood number for a droplet transported in a turbulent pipe flow is
calculated by an empirical correlation as follows (Kress and Keyes, 1973):

$\displaystyle Sh=0.34\left(\dfrac{d_{p}}{D}\right)^{2}Re^{0.94}Sc^{0.5}$ (30)

where $Sc=\nu_{f}/D$ is Schmidt number.
Note that both the concentration and saturation concentration of a dissolved
gas in bulk water change along a pipeline. At a moderate pressure, the
saturation concentration is calculated by Henry’s law (e.g. Bird et al.,
2002), according to which the saturation gas concentration is proportional to
the pressure. The dissolved gas concentration in bulk fluid is calculated from
the mass conservation for a gas phase. If at the initial time moment, no gas
is dissolved in a liquid, the dissolved gas concentration evolution with a
decrease in a droplet holdup is calculated as:

$\displaystyle C_{\infty}=\dfrac{\rho_{d}(\phi_{0}-\phi)}{1-\phi}$ (31)

where $\phi=\displaystyle\sum_{i=1}^{M}N_{i}x_{i}$ is the volume concentration
of a dispersed phase, $\phi_{0}$ is the volume concentration of a dispersed
phase at the initial time moment.

##  3 Results and discussion

The advection-diffusion population balance equation Eq. (5) has been solved
numerically. A MATLAB code was developed for this purpose. It is to be noted
that in the simulations we have seen that calculations with the time steps of
$10^{-4}$ and $10^{-5}$ provided virtually indistinguishable results, and thus
the calculations are independent of the grid size. We also found that a
further reduction in the radial mesh size and an increase in a number of
droplet size fractions do not noticeably affect the computational results
either. The code was validated against the experimental data given in
Zendehboudi et al. (2013). Zendehboudi et al. (2013) measured changes in
Sauter diameter along the pipe length for the volume fraction of $0.05$ under
the following conditions: Pressure is $70$ bar, temperature is $25^{\circ}C$
and brine flow rate in the range $0.008-0.064$ m3/s. At these conditions the
solubility of $CO_{2}$ is approximately $5\%$. The solubility decreases with
temperature and increases with pressure. The reservoir temperature is higher
than that used in the experiment, but the pressure is also higher; therefore,
the solubility value used for modeling is in agreement with possible reservoir
conditions. When designing real injection systems, the maximum dissolved
$CO_{2}$ amount should never exceed the solubility at reservoir conditions to
prevent escape of carbon dioxide from brine. In the experiments, the liquid
$CO_{2}$ was merged with a flow of brine phase. The droplets were recorded and
tracked using high-speed cameras. In the present work, for calculation of the
$CO_{2}$ concentration we employed the same assumption that Zendehboudi et al.
(2013) used in their calculations: The saturation concentration in brine was
evaluated as $Cs=0.85$ $Cs$ (in pure water). This assumption provides
reasonable data for $Cs$ at temperatures in the range of $20-100^{\circ}C$ and
pressures in the range of $0-80$ bar for salinities of $0.5-1.4$ mol/kg. The
computed Sauter diameter distributions along a pipe for different mean flow
velocities at the initial mean droplet concentration $\phi_{0}=0.05$ were
matched to the measured data by tuning the parameters $K$ in Eq. (16)
describing droplet breakup rate and $C_{2}$ in Eq. (26) defining the
coalescence probability. As it was discussed above, the parameter $K$ was not
reliably determined in the past, whereas the parameter $C_{2}$ depends on the
chemical composition of fluids composing a dispersion. The best fitting was
obtained at $K=0.1$ and $C_{2}=10^{13}$. One can see that the computational
results correlate well with the measured data in Fig. 2. The fitting is not
very accurate along the first half of the pipe length, whereas the second half
is characterized by a closer fit. These observations are primarily explained
by experimental inaccuracy. The droplet size distributions along the initial
pipe section are wide and rapidly changed; therefore, in the experiment, the
analysis of images obtained by using high speed cameras does not allow a
highly accurate evaluation of the droplet Sauter diameters. The second half of
the pipe is characterized with nearly steady-state droplet size distributions,
which are relatively narrow, and therefore, the Sauter diameters were
determined with a better fit than those in the first half of the pipe.

Figure 2: Distributions of computed droplet Sauter diameters along a pipe at
different flow velocities versus the experimental data of Zendehboudi et al.
(2013).

Both droplet breakup and coalescence strongly affect the dissolution process.
To illustrate the importance of these phenomena, in Fig. 3 we showed the rate
at which droplet sizes decrease along a pipe if breakup and coalescence are
absent and droplet size is reduced only due to dissolution.

Figure 3: Evolution of droplet sizes along a pipe at different flow velocities
accounting only for dissolution.

One can see that droplet sizes change slowly and do not approach steady-state.
The smaller the droplets, the higher the overall dissolution rate - which is,
to a large extent, determined by specific droplet surface area.

In Figs. 4 and 5 we illustrate an effect of the droplet volume fraction on the
dissolution process. The calculations were conducted at the three different
initial volume fractions $\phi_{0}=0.01,0.05,0.1$. The mean flow velocity was
assumed to be the same for all the computations, $v=3.63$ m/s. Fig. 4 shows
how the droplet volume fraction changes along the pipe. In Fig. 5, one can see
the Sauter diameter evolution. The higher the droplet concentration, the
closer the Sauter diameter approaches steady-state. This observation is
explained as follows: A higher droplet concentration leads to a larger mass
flux (from the dispersed to the continuous phase) which causes a faster
increase in dissolved $CO_{2}$ concentration. This leads to the dissolved gas
concentration rapidly approaching the saturation concentration, resulting in a
slower mass transfer rate.

Figure 4: Evolution of droplet volume fractions along a pipe with different
initial droplet volume concentrations at the fixed mean flow velocity $v=3.63$
m/s. Figure 5: Evolution of droplet sizes along a pipe for different initial
droplet volume concentrations at a fixed mean flow velocity of $v=3.63$ m/s.

In Fig. 6 one can see an effect of the flow velocity on the change in the
dissolved gas concentration along a pipe at the fixed initial droplet
concentration $\phi_{0}=0.05$. The higher the flow velocity, the higher the
mass transfer rate, resulting in a higher dissolved $CO_{2}$ concentration in
water. This observation is explained by the following factors: 1) the higher
the Reynolds number, the higher the mass transfer coefficient between a
droplet and a surrounding liquid; and 2) the higher the flow velocity, the
smaller the droplets - this leads to an increase in the droplet surface area.

Figure 6: Evolution of dissolved $CO_{2}$ concentration by mass ($kg$) with
different flow velocities at the fixed initial water content $\phi_{0}=0.05$.

Based on the numerical study conducted, it is possible to conclude that the
process of $CO_{2}$ dissolution in a pipe flow can be optimized by a proper
selection of parameters such as the droplet concentration and the flow
velocity. Note, the pipe diameter is also an important parameter allowing to
vary the flow velocity if the flow rate is constrained. Overall, the
dissolution process optimization is a relatively simple problem and the code
presented here makes the optimization a rather straightforwardly handled task.
Although the model developed is based on first principles and has a
significantly higher predictive capability, it is capable of running
relatively quickly on a regular laptop. Therefore, it is possible to formulate
a simple approach, which will allow using the code developed to design an
optimum dissolution process with a minimum time expense. Let us describe this
approach.

In engineering applications, a $CO_{2}$ flow rate is known. Pipe diameter is
not, usually, varied significantly. Flow rate, determining $CO_{2}$ droplet
concentration, is also restricted. Pipe length can usually be varied to a
certain extent. Thus, the optimal combination of a water flow rate and a pipe
length can be determined by a few computations. The following calculation
strategy could be employed:

  * •

Conduct a computation using the most reasonable pipe length and the maximum
possible flow rate. Results will show a degree of droplet dissolution.

  * •

If both the droplet volume fraction and droplet sizes are sufficiently small,
then further computations could be done at reduced flow rates. The smallest
possible flow rate which allows to get satisfactory mixture parameters at the
pipe end should be selected.

  * •

If a dissolution degree is not satisfactory, the only option available, in
this case, is an increase in the pipe length. Otherwise, if an acceptable
increase in the pipe length leads to acceptable droplet sizes, then the
problem is solved. If the acceptable pipe length is not sufficient for
dissolution, then a desired design is impossible to achieve. Maybe, using a
system composed of two pipes could be employed.

To demonstrate the design process, we conducted computation in a pipe of $700$
m at a constant $CO_{2}$ flow rate corresponding to $\phi_{0}=0.05$ of
$CO_{2}$ droplets at a mean flow velocity equal to $v=1.45$ m/s. The other
computations were done at the same $CO_{2}$ volume flow rate but for the
higher mixture flow velocities, $v=2.5$ m/s ($\phi_{0}=0.029$) and $v=3.63$
m/s ($\phi_{0}=0.01997$), respectively. The computational results are shown in
Fig. 7. We can see that very fine droplet sizes are achievable and conclude
that an increase in the brine flow rate along with an increase in the pipe
length can cause substantial reduction in droplet sizes due to three main
factors: 1) an increase in the flow velocity leads an increase in a turbulence
energy dissipation rate level that in turn causes an increase in the droplet
breakup rate; 2) the higher the flow velocity, the higher the droplet-fluid
mass transfer coefficient; 3) the lower the droplet concentration, the higher
the difference between the $CO_{2}$ saturation concentration and the
concentration of $CO_{2}$ dissolved in a fluid causing a higher droplet-fluid
mass transfer rate. Thus, in the majority of practical applications, a desired
droplet size can be provided.

Figure 7: Evolution of droplet sizes along a pipe for different flow
velocities with corresponding initial water contents.

##  4 Conclusion

Thus, as a critical part of optimization of carbon dioxide underground
disposal technology, a comprehensive model of ex-situ $CO_{2}$ droplet
dissolution in a turbulent pipe water flow, accompanied with droplet breakup
and coalescence, has been developed. Modeling has been reduced to a numerical
solution of the steady-state advection-diffusion-population balance equation.
The Prandtl mixing length model was employed for the modeling of the velocity
distribution across a pipe. The turbulent diffusivity distribution along a
pipe radius was calculated by the empirical correlation found in open
literature. The turbulence energy dissipation rate, needed for population
balance modelling, was calculated by the analytical equation. The population
balance model developed accounts for droplet breakup, coalescence and
turbulent mass transfer. Although both the breakup and coalescence terms of
the discretized population balance equation were computed using known
expressions, the dissolution term was derived in the present work. The semi-
empirical models of droplet breakup, coalescence and dissolution were
employed.
The model developed has been validated against the experimental data. The
computed distributions of mean droplet Sauter diameter of liquid $CO_{2}$
along a pipe were compared with the measured distributions. The two model
parameters were tuned for fitting the experimental results. The numerical
studies, conducted by the validated model, showed that the droplet dissolution
significantly speeds up with an increase in the flow velocity. An increase in
the initial droplet concentration slows the dissolution process.
The model developed can be used for an efficient dissolution process design
and optimization.

##  5 Appendix

###  5.1 Dissolution

The dissolution term of the population balance equation was derived by
considering the mass balance for a droplet of volume $x_{i}$. The mass
transfer rate is assumed to be described by the function $\dfrac{dx_{i}}{dt}$.
We can formulate the mass balance equation between two particles of volumes
$x_{i}$ and $x_{i+1}$ due to dissolution during the time $dt$ as follows:

$\displaystyle\dfrac{\partial\hat{N}_{i+1}}{\partial
t}x_{i}dt+N_{i+1}\dfrac{dx_{i+1}}{dt}dt=\dfrac{\partial\hat{N}_{i+1}}{\partial
t}x_{i+1}dt$ (32)

where $\dfrac{\partial\hat{N}_{i+1}}{\partial t}$ is the droplet number
concentration rate of change of the size fraction $i+1$ only due to
dissolution of droplets of this size fraction. The first term is equal to the
volume transferred from droplets of the volume $x_{i+1}$ to droplets of the
volume $x_{i}$ due to dissolution of droplets of the $i+1$ size fraction. The
second term is equal to the dissolved volume of droplets of the size fraction
$i+1$. The right-hand side term is equal to the volume of droplets of the size
fraction $i+1$ which disappeared due to dissolution. Note, the dissolution
rate is assumed to be positive in our analysis; therefore, absolute value
parentheses are used in Eq. (33) and further. From Eq. (32) we obtain:

$\displaystyle\dfrac{d\hat{N}_{i+1}}{dt}$
$\displaystyle=\dfrac{N_{i+1}}{x_{i+1}-x_{i}}\bigg{|}\dfrac{dx_{i+1}}{dt}\bigg{|}$
(33)

Then, the volume flux from droplets of the size fraction $i+1$ to droplets of
the size fraction $i$ due to dissolution is:

$\displaystyle
q_{i}^{+}=\dfrac{N_{i+1}x_{i}}{x_{i+1}-x_{i}}\bigg{|}\dfrac{dx_{i+1}}{dt}\bigg{|},$
(34)

and the volume flux from droplets of the size fraction $i$ (indicating
droplets disappearance due to dissolution) to droplets of the size fraction
$i-1$ is as follows:

$\displaystyle
q_{i}^{-}=\dfrac{N_{i}x_{i}}{x_{i}-x_{i-1}}\bigg{|}\dfrac{dx_{i}}{dt}\bigg{|}.$
(35)

Hence Eq. (34) and Eq.(35) gives the total rate of number concentration change
of $x_{i}$ due to dissolution,

$\left(\dfrac{\partial N_{i}}{\partial
t}\right)_{diss}=\dfrac{N_{i+1}}{x_{i+1}-x_{i}}\bigg{|}{\dfrac{dx_{i+1}}{dt}}\bigg{|}-\dfrac{N_{i}}{x_{i}-x_{i-1}}\bigg{|}{\dfrac{dx_{i}}{dt}}\bigg{|}.$
(36)

For both the smallest and the largest size fractions this equation is modified
as follows:

$\displaystyle\dfrac{\partial N_{1}}{\partial t}$
$\displaystyle=\dfrac{N_{2}}{x_{2}-x_{1}}\bigg{|}\dfrac{dx_{2}}{dt}\bigg{|}-\dfrac{N_{1}}{x_{1}}\bigg{|}\dfrac{dx_{1}}{dt}\bigg{|}\hskip
14.22636pt\text{for the first size fraction $x_{1}$}$ (37)
$\displaystyle\dfrac{\partial N_{m}}{\partial t}$
$\displaystyle=-\dfrac{N_{m}}{x_{m}-x_{m-1}}\bigg{|}\dfrac{dx_{m}}{dt}\bigg{|}\hskip
14.22636pt\text{for the last size fraction $x_{m}$}$ (38)

Financial support for this work provided by Natural Sciences and Engineering
Research Council of Canada (NSERC).

Declaration of interest: None

## References

  * Boden et al. 2017 Boden, T. A.; Marland, G.; Andres, R. J. _Global, Regional, and National Fossil-Fuel CO 2 Emissions_; Carbon Dioxide Information Analysis Center, Oak Ridge National Laboratory, U.S. Department of Energy: Oak Ridge, Tenn., U.S.A., 2017
  * Metz et al. 2005 Metz, B.; Davidson, O.; de Coninck, H.; Loos, M.; Meyer, L. _Carbon Dioxide Capture and Storage: Special Report of the Intergovernmental Panel on Climate Change_ ; Cambridge University Press: New York, 2005
  * Herzog 2001 Herzog, H. J. What future for carbon capture and sequestration? _Environ. Sci. Technol._ 2001, _35 (7)_ , 148A–153A
  * Jenkins et al. 2012 Jenkins, C. R. et al. Safe storage and effective monitoring of CO2 in depleted gas fields. _Proc. Natl. Acad. Sci. U.S.A._ 2012, _109 (2)_ , E35–E41
  * Shi and Durucan 2005 Shi, J.-Q.; Durucan, S. CO2 Storage in Deep Unminable Coal Seams. _Oil & Gas Science and Technology – Rev. IFP_ 2005, _60 (3)_ , 547–558
  * Celia et al. 2015 Celia, M. A.; Bachu, S.; Nordbotten, J. M.; Bandilla, K. W. Status of CO2 storage in deep saline aquifers with emphasis on modeling approaches and practical simulations. _Water Resour. Res._ 2015, _51_ , 6846–6892
  * Haugan and Joos 2004 Haugan, P. M.; Joos, F. Metrics to assess the mitigation of global warming by carbon capture and storage in the ocean and in geological reservoirs. _Geophys. Res. Lett._ 2004, _31_ , L18202
  * International Energy Agency (2008) IEA International Energy Agency (IEA), _CO 2 Capture and Storage_; IEA Publications, 2008
  * van der Meer 1993 van der Meer, L. G. H. The Conditions Limiting CO2 Storage in aquifers. _Energy Convers. Manage._ 1993, _34_ , 959–966
  * Lindeberg 1997 Lindeberg, E. Escape of CO2 from aquifers. _Energy Convers. Manage._ 1997, _38_ , S235–S240
  * Emami-Meybodi et al. 2015 Emami-Meybodi, H.; Hassanzadeh, H.; Green, C. P.; Ennis-King, J. Convective dissolution of $CO_{2}$ in saline aquifers: Progress in modeling and experiments. _Int. J. Greenh. Gas Con._ 2015, _40_ , 238–266
  * Hassanzadeh et al. 2009 Hassanzadeh, H.; Pooladi-Darvish, M.; Keith, D. W. Accelerating $CO_{2}$ dissolution in saline aquifers for geological storage - Mechanistic and sensitivity studies. _Energy & Fuels_ 2009, _23(6)_ , 3328–3336
  * Keith et al. 2005 Keith, D. W.; Hassanzadeh, H.; Pooladi-Darvish, M. _Greenhouse Gas Control Technologies_ , 7th ed.; Elsevier Science Ltd: Oxford, 2005; pp 2163 – 2167
  * Zirrahi et al. 2013a Zirrahi, M.; Hassanzadeh, H.; Abedi, J. Modeling of $CO_{2}$ dissolution by static mixers using back flow mixing approach with application to geological storage. _Chem. Eng. Sci._ 2013a, _104_ , 10–16
  * Zirrahi et al. 2013b Zirrahi, M.; Hassanzadeh, H.; Abedi, J. The laboratory testing and scale-up of a downhole device for $CO_{2}$ dissolution acceleration. _Int. J. Greenh. Gas Con._ 2013b, _16_ , 41–49
  * Shafaei et al. 2012 Shafaei, M. J.; Abedi, J.; Hassanzadeh, H.; Chen, Z. Reverse gas-lift technology for $CO_{2}$ storage into deep saline aquifers. _Energy_ 2012, _45(1)_ , 840–849
  * Leonenko and Keith 2008 Leonenko, Y.; Keith, D. W. Reservoir engineering to accelerate the dissolution of CO2 stored in aquifers. _Environ. Sci. Technol._ 2008, _42 (8)_ , 2742–2747
  * Zendehboudi et al. 2011 Zendehboudi, S.; Khan, A.; Carlisle, S.; Leonenko, Y. Ex Situ Dissolution of CO2: A New Engineering Methodology Based on Mass-Transfer Perspective for Enhancement of CO2 Sequestration. _Energy & Fuels_ 2011, _25_ , 3323–3333
  * Cholewinski and Leonenko 2013 Cholewinski, A.; Leonenko, Y. Ex-situ Dissolution of CO2 for Carbon Sequestration. _Energy Procedia_ 2013, _37_ , 5427–5434
  * Zendehboudi et al. 2013 Zendehboudi, S.; Shafiei, A.; Bahadori, A.; Leonenko, Y.; Chatzis, I. Droplets evolution during ex situ dissolution technique for geological CO2 sequestration: Experimental and mathematical modelling. _Int. J. Greenh. Gas Con._ 2013, _13_ , 201–214
  * Levich 1962 Levich, V. G. _Physicochemical Hydrodynamics_ ; Prentice-Hall, Englewood Cliffs, NJ, 1962
  * Eskin 2012 Eskin, D. A Simple Model of Particle Diffusivity in Horizontal Hydrotransport Pipelines. _Chem. Eng. Sci._ 2012, _82_ , 84–94
  * Eskin et al. 2017a Eskin, D.; Taylor, S.; Ma, S. M.; Abdallah, W. Modeling droplet dispersion in a vertical turbulent tubing flow. _Chem. Eng. Sci._ 2017a, _173_ , 12–20
  * Schlichting and Gersten 2000 Schlichting, H.; Gersten, K. _Boundary-Layer Theory_ ; Springer, Berlin, 2000
  * Johansen 1991 Johansen, S. T. The deposition of particles on vertical walls. _Int. J. Multiph. Flow_ 1991, _17_ , 355––376
  * Kumar and Ramkrishna 1996 Kumar, S.; Ramkrishna, D. On the solution of population balance equations by discretization-I. A fixed pivot technique. _Chem. Eng. Sci._ 1996, _51_ , 1311–1332
  * Eskin et al. 2017b Eskin, D.; Taylor, S.; Dingzheng, Y. Modeling of droplet dispersion in a turbulent Taylor-Couette flow. _Chem. Eng. Sci._ 2017b, _161_ , 36–47
  * Bird et al. 2002 Bird, B. R.; Stuart, W. E.; Lightfoot, E. N. _Transport Phenomena_ ; John Wiley & Sons Inc, New York/Toronto, 2002
  * Liao and Lucas 2010 Liao, Y.; Lucas, D. A literature review on mechanisms and models for the coalescence process of fluid particles. _Chem. Eng. Sci._ 2010, _65_ , 2851–2864
  * Coulaloglou and Tavlarides 1977 Coulaloglou, C.; Tavlarides, L. Description of interaction processes in agitated liquid–liquid dispersions. _Chem. Eng. Sci._ 1977, _32_ , 1289–1297
  * Laakkonen et al. 2006 Laakkonen, M.; Alopaeus, V.; Aittamaa, J. Validation of bubble breakage, coalescence and mass transfer models for gas–liquid dispersion in agitated vessel. _Chem. Eng. Sci._ 2006, _61_ , 218–228
  * Kress and Keyes 1973 Kress, T. S.; Keyes, J. J. Liquid phase controlled mass transfer to bubbles in cocurrent turbulent pipeline flow. _Chem. Eng. Sci._ 1973, _28_ , 1809–1823