∎

11institutetext: Ø.S. Klemetsdal 22institutetext: Tel.: +47 98 43 86 39
Norwegian University of Science and Technology, Norway
22email<EMAIL_ADDRESS>33institutetext: A.F. Rasmussen
44institutetext: SINTEF Digital, Norway
44email<EMAIL_ADDRESS>55institutetext: O. Møyner
66institutetext: Norwegian University of Science and Technology/SINTEF
Digital, Norway
66email<EMAIL_ADDRESS>77institutetext: K.-A. Lie 88institutetext:
Norwegian University of Science and Technology/SINTEF Digital, Norway
88email<EMAIL_ADDRESS>
# Efficient Reordered Nonlinear Gauss–Seidel Solvers With Higher Order For
Black-Oil Models

Øystein S. Klemetsdal  Atgeirr F. Rasmussen  Olav Møyner  Knut-Andreas Lie

(Received: date / Accepted: date)

###### Abstract

The fully implicit method is the most commonly used approach to solve black-
oil problems in reservoir simulation. The method requires repeated
linearization of large nonlinear systems and produces ill-conditioned linear
systems. We present a strategy to reduce computational time that relies on two
key ideas: (i) a sequential formulation that decouples flow and transport into
separate subproblems, and (ii) a highly efficient Gauss–Seidel solver for the
transport problems. This solver uses intercell fluxes to reorder the grid
cells according to their upstream neighbors, and groups cells that are
mutually dependent because of counter-current flow into local clusters. The
cells and local clusters can then be solved in sequence, starting from the
inflow and moving gradually downstream, since each new cell or local cluster
will only depend on upstream neighbors that have already been computed.
Altogether, this gives optimal localization and control of the nonlinear
solution process.

This method has been successfully applied to real-field problems using the
standard first-order finite volume discretization. Here, we extend the idea to
first-order dG methods on fully unstructured grids. We also demonstrate proof
of concept for the reordering idea by applying it to the full simulation model
of the Norne oil field, using a prototype variant of the open-source OPM Flow
simulator.

##  1 Introduction

Reservoir simulation requires the solution of large systems of nonlinear
partial differential equations to compute fluid pressure, phase saturations,
and component concentrations/mass fractions. In industry-grade simulations, it
is most common to use implicit temporal discretization to overcome severe CFL
restrictions arising because of grids with high aspect ratios, orders of
magnitude variations in petrophysical properties, and large variations in
fluid velocities from stagnant and slow-moving flow regions to high-flow
regions around wells. For cases with weak coupling between the fluid pressure
and the transport of conserved quantities, one can observe significant
reduction in computational costs by splitting the flow model into a pressure
equation and a set of transport equations and solve them sequentially
[Trangenstein and Bell, 1989, Watts, 1986]. To this end, it is common to use
specialized solvers since the two subproblems typically have very different
mathematical character: pressure equations are often close to elliptic,
whereas transport equations have a strong hyperbolic character; see e.g., Bell
et al. [1986] and Lie [2016].

The pressure equation is usually discretized by a first-order, two-point flux-
approximation scheme, giving a discrete problem that can be efficiently solved
by a multigrid [Trottenberg et al., 2000, Gries et al., 2014] or multiscale
method [Lie et al., 2017]. For the transport equations, a number of methods
aim to accelerate computations by utilizing inherent locality and co-current
flow properties of hyperbolic equations. Examples include streamline
simulation [Datta-Gupta and King, 2007, Bratvedt et al., 1996], a priori
estimation of nonzero update regions [Sheth and Younis, 2017], and use of
interface-localized trust regions to determine safe saturation updates
[Møyner, 2017, Klemetsdal et al., 2019a].

Reordering methods rely on the important observation that the transport of
fluid phases is always unidirectional along streamlines in the absence of
gravity and capillary forces. Consequently, it is possible to order the grid
cells so that the discretization matrix for the transport equations becomes
lower triangular and can be solved very efficiently in a cell-by-cell manner.
A branch of these methods builds upon the Cascade method [Appleyard and
Cheshire, 1982] and reorders grid cells based on the fluid potential [Kwok and
Tchelepi, 2007, Shahvali and Tchelepi, 2013]. Another branch uses topological
traversal of the intercell flux graph [Natvig and Lie, 2008, Lie et al.,
2012b, 2014]. The latter has the advantage that it is easy to implement using
methods from standard graph theory. When gravity and capillary forces are
present, the flow field will have regions of counter-current flow, which show
up as cycles of mutually connected cells in the discretized flow equations. If
the cells in these cycles are grouped into supernodes and solved for
simultaneously, we can still use the same reordering idea.

Numerical smearing and grid-orientation effects are well-known problems in
reservoir simulation. These effects are particularly evident for displacement
fronts with little or no self-sharpening effects, as are often found in
multicomponent models describing various mechanisms for enhanced oil recovery,
e.g., polymer injection and miscible flooding. We can mitigate these effects
by increasing the grid resolution and/or the formal order of the spatial
discretization. Both approaches increase the computational cost by increasing
the number of unknowns and/or the nonlinearity of the discretized equations,
and it is therefore even more important to have a nonlinear solver with high
efficacy. Most high-resolution methods—i.e., methods that deliver high formal
order on smooth parts of the solution and stable propagation of
discontinuities—rely on some kind of spatial reconstruction using cell-based
averaged quantities from a ring of immediate cell neighbors. The resulting
stencils involve cells both in the upstream and downstream direction. This
means that the resulting discretization graph does not reflect the same
unidirectional properties as the underlying flow equations for co-current
flow.

Alternatively, one can use discontinuous Galerkin methods (dG), which were
first extended to general systems of hyperbolic conservation laws by Cockburn
and Shu [1989, 1991]. When combined with a single-point upstream mobility
scheme for flux evaluation, these methods preserve the causality of the
continuous flow equations. That is, the corresponding stencil is restricted to
cells in the upstream direction in regions of co-current flow. After
application of flux-based reordering, the discrete nonlinear flow equations
are permuted to block-triangular form, with small blocks representing
individual cells from regions of co-current flow and larger blocks
representing regions of counter-current flow. This gives a natural
localization: You start at the inflow locations and solve the nonlinear
equations block-by-block toward the location of outflow. For blocks consisting
of multiple cells, you can either solve for all mutually dependent unknowns
simultaneously, or you can use an effective Gauss–Seidel solver that
decomposes the problem to an iteration over a sequence of single-cell problems
[Lie et al., 2014]. The advantage of the Gauss–Seidel approach is that it
relies entirely on single-cell nonlinear solvers, which are simpler to
optimize.

In this work, we combine intercell flux reordering methods with higher-order
dG methods. This idea has previously been studied by Natvig et al. [2007] and
Eikemo et al. [2009] for passive advection problems and by Natvig and Lie
[2008] for incompressible problems in two-phase, three-phase, and two-
phase–three-component flow without gravity and capillary forces. Later, the
first-order variant of the method was extended to include gravity [Lie et al.,
2012b] and to polymer flooding [Lie et al., 2014] with compressibility and
gravity effects. Herein, we study the widely used family of black-oil
equations and present a method that can handle all relevant flow effects
including hysteresis and capillary forces as seen in the simulation model of
the Norne oil field The Open Porous Media Initative(2015) [OPM]. We also
discuss how to formulate high-order dG discretization for general polyhedral
grids, introduce a simple order-reduction method to prevent creating spurious
oscillations, and present a new method based on blocks of cells, which is more
parallel and cache efficient than solving for a single cell at a time.

##  2 Governing equations

The black-oil model describes conservation of three pseudo-components (water,
oil, and gas), which at reservoir conditions can distribute in three phases
(aqueous, oleic, and gaseous):

$\displaystyle\partial_{t}\left(\phi
b_{w}S_{w}\right)+\nabla\cdot\left(b_{w}\vec{v}_{w}\right)-b_{w}q_{w}=0,$ (1)
$\displaystyle\partial_{t}\left(\phi\left[b_{o}S_{o}+b_{g}r_{v}S_{g}\right]\right)$
$\displaystyle\,\,+\nabla\cdot\left(b_{o}\vec{v}_{o}+b_{g}r_{v}\vec{v}_{g}\right)-(b_{o}q_{o}+b_{g}r_{v}q_{g})=0,$
$\displaystyle\partial_{t}\left(\phi\left[b_{g}S_{g}+b_{o}r_{s}S_{o}\right]\right)$
$\displaystyle\,\,+\nabla\cdot\left(b_{g}\vec{v}_{g}+b_{o}r_{s}\vec{v}_{o}\right)-(b_{g}q_{g}+b_{o}r_{s}q_{o})=0.$

Here, $S_{\alpha}$ denotes saturation, $\vec{v}_{\alpha}$ the macroscopic
Darcy velocity, and $q_{\alpha}$ sources and sinks of phase $\alpha$.
Shrinkage factors $b_{\alpha}$ model pressure-dependent densities
$\rho_{\alpha}$, and the gas-oil and oil-gas ratios $r_{s}$ and $r_{v}$ model
the volume of gas dissolved in oil and oil vaporized in gas, respectively,
both at standard conditions. The phase velocity $\vec{v}_{\alpha}$ is given by
the multiphase extension of Darcy’s law:

$\vec{v}_{\alpha}=-\lambda_{\alpha}\boldsymbol{\mathsf{K}}\left(\nabla
p_{\alpha}-\rho_{\alpha}g\nabla z\right),$

where $\lambda_{\alpha}=k_{r\alpha}/\mu_{\alpha}$ is the mobility of phase
$\alpha$; relative permeability $k_{r\alpha}$ models the reduced mobility of
one phase in the presence of another, whereas $\mu_{\alpha}$ as usual denotes
the viscosity of phase $\alpha$. The model is closed by assuming that the
fluid phases fill up the entire pore space, $S_{w}+S_{o}+S_{g}=1$, and that
the phase pressures are related through capillary pressures:

$p_{o}=p_{w}+P_{cow}(S_{w},S_{o}),\quad p_{g}=p_{o}+P_{cgo}(S_{o},S_{g}).$

##  3 Sequential splitting

We use backward Euler to discretize in time. On semi-discrete form, the
equation describing conservation of mass for water (1) at time $n+1$ then
reads

$\displaystyle\mathcal{R}_{w}=$ $\displaystyle\frac{1}{\Delta
t}\left(\left[\phi b_{w}S_{w}\right]^{n+1}-\left[\phi
b_{w}S_{w}\right]^{n}\right)$ (2)
$\displaystyle\qquad+\nabla\cdot\left(b_{w}\vec{v}_{w}\right)^{n+1}-(b_{w}q_{w})^{n+1}=0,$

where we have introduced the time-step length $\Delta t=t^{n+1}-t^{n}$, and
superscripts $n$ and $n+1$ refer to the time step. Semi-discrete forms for the
oil ($\mathcal{R}_{o}$) and gas ($\mathcal{R}_{g}$) equations are analogous.
We multiply each equation by the following factors

$\displaystyle\omega_{w}$ $\displaystyle=\frac{1}{b_{w}^{n+1}},$
$\displaystyle\omega_{o}$
$\displaystyle=\frac{1}{1-r_{s}^{n+1}r_{v}^{n+1}}\left(\frac{1}{b_{o}^{n+1}}-\frac{r_{s}^{n+1}}{b_{g}^{n+1}}\right),$
$\displaystyle\omega_{g}$
$\displaystyle=\frac{1}{1-r_{s}^{n+1}r_{v}^{n+1}}\left(\frac{1}{b_{g}^{n+1}}-\frac{r_{v}^{n+1}}{b_{o}^{n+1}}\right),$

and sum to eliminate all terms involving saturations at time step $n+1$. This
gives us a pressure equation

$\mathcal{R}_{p}=\omega_{w}\mathcal{R}_{w}+\omega_{o}\mathcal{R}_{o}+\omega_{g}\mathcal{R}_{g}.$
(3)

To obtain a fully discrete formulation, we introduce a grid covering our
computational domain with $N$ non-overlapping polyhedral cells
$\\{\Omega_{i}\\}_{i=1}^{N}$. We denote the common interface of cell $i$ and
$j$ by $\Gamma_{ij}$ and the volumetric flux from cell $i$ to cell $j$ by

$v_{\alpha,ij}=\int_{\Gamma_{ij}}\vec{v}_{\alpha}\cdot\vec{n}_{ij}\,\mathrm{d}\sigma.$

Here, $\vec{n}_{ij}$ is the unit normal from cell $i$ to cell $j$. Notice that
conservation of mass requires that $v_{\alpha,ij}=-v_{\alpha,ji}$. In the
following, we assume that we have a pressure solver that can solve the fully
discrete version of (3) to yield a total velocity field
$\vec{v}=\vec{v}_{w}+\vec{v}_{o}+\vec{v}_{g}$, given as a set of constant
fluxes $v_{ij}=v_{w,ij}+v_{o,ij}+v_{g,ij}$ over each cell interface. In the
numerical experiments, unless otherwise noted, we use the standard black-oil
pressure solver from the open-source MATLAB Reservoir Simulation Toolbox [Lie,
2016].

In a sequential solution procedure, we solve (3) with fixed saturations to
obtain pressures and the total Darcy velocity $\vec{v}$. Given the total flux,
we compute new phase fluxes using a fractional flow formulation. These are
defined from the following formula in the semi-continuous case (neglecting
capillary effects)

$\displaystyle\vec{v}_{\alpha}=f_{\alpha}\left(\vec{v}+\boldsymbol{\mathsf{K}}\sum_{\beta\neq\alpha}\lambda_{\beta}\left(\rho_{\alpha}-\rho_{\beta}\right)g\nabla
z\right),$ (4)

where we have introduced the fractional flow function

$f_{\alpha}(S_{w},S_{o},S_{g})=\frac{\lambda_{\alpha}(S_{\alpha})}{\lambda_{w}(S_{w})+\lambda_{o}(S_{o})+\lambda_{g}(S_{g})}.$

Capillary pressure differences are accounted for with a term similar to the
gravity term in (4). Using this, we solve the transport equations (1) with
fixed pressure to obtain new saturations (or solution ratios). This introduces
a splitting error proportional to the time step $\Delta t$ that can be
significant for cases with strong coupling between pressure and saturation.
One can also make the solution converge towards the fully implicit solution by
adding outer iterations, as described by e.g, Jenny et al. [2006].

##  4 Discontinuous Galerkin discretization

For simplicity, we only describe the weak formulation for the water equation
without capillary and gravity effects, even though gravity effects are
included in our solver. Capillary effects are also included, but currently
only supported in our first-order transport solver. Moreover, we drop the
phase subscript $w$ and the time superscript $n+1$. Equation (2) now reads

$\mathcal{R}_{w}=\tfrac{1}{\Delta t}\left(\left[\phi bS\right]-\left[\phi
bS\right]^{n}\right)+\nabla\cdot\left(bf\vec{v}\right)-(bq)=0,$ (5)

where we recall that the fractional flow function $f$ depends on all phase
saturations. We multiply by a test function $\psi$ in a function space $V$ of
arbitrarily smooth functions and integrate (by parts) over cell $\Omega_{i}$:

$\displaystyle\frac{1}{\Delta t}\int_{\Omega_{i}}\Big{(}(\phi bS)-(\phi
bS)^{n}\Big{)}\psi\,\mathrm{d}V-$
$\displaystyle\int_{\Omega_{i}}bf\vec{v}\cdot\nabla\psi\,\mathrm{d}V$
$\displaystyle+\int_{\partial\Omega_{i}}\left[bf\psi\right]\vec{v}\cdot\vec{n}\,\mathrm{d}\sigma-$
$\displaystyle\int_{\Omega_{i}}(bq)\psi\,\mathrm{d}V=0.$

Here, the square brackets in the surface integral signify that the integrand
is possibly discontinuous across the cell boundary. Note that we can express
the surface integral as

$\int_{\partial\Omega_{i}}\left[bf\psi\right]\vec{v}\cdot\vec{n}\,\mathrm{d}\sigma=\sum_{j\in\mathcal{N}(i)}\int_{\Gamma_{ij}}\left[bf\psi\right]\vec{v}\cdot\vec{n}_{ij}\,\mathrm{d}\sigma,$

where $\mathcal{N}(i)$ denotes the set of all cells sharing a common interface
with cell $i$. We define the upwind operator

$\displaystyle u[y,\vec{v}_{\alpha}\cdot\vec{n}_{ij};x]$
$\displaystyle=\begin{cases}y_{ij}\vec{v}_{\alpha}\cdot\vec{n}_{ij},&\text{if}\quad\vec{v}_{\alpha}\cdot\vec{n}_{ij}>0,\\\
y_{ji}\vec{v}_{\alpha}\cdot\vec{n}_{ij},&\text{if}\quad\vec{v}_{\alpha}\cdot\vec{n}_{ij}\leq
0,\end{cases}$ $\displaystyle\text{where }y_{ij}$
$\displaystyle=\lim_{x^{\prime}\to x}y(x^{\prime})\text{ for
}x^{\prime}\in\Omega_{i}.$

Subscripts $ij$ and $ji$ denote the limiting interface value of
$y(x^{\prime})$ for $x\in\Gamma_{ij}$ as $x^{\prime}$ approaches $x$ from
inside cell $i$ or cell $j$, respectively. If the coordinate is not important,
we simply write $u[y,\vec{v}_{\alpha}\cdot\vec{n}_{ij}]$. Using this, we write

$\int_{\partial\Omega_{i}}\left[bf\psi\right]\vec{v}\cdot\vec{n}\,\mathrm{d}\sigma=\sum_{j\in
N(i)}\int_{\Gamma_{ij}}u\left[b,f\vec{v}\cdot\vec{n}_{ij}\right]\psi\,\mathrm{d}\sigma.$

Since the fractional flow function is always non-negative, $\vec{v}$ and
$\vec{v}_{\alpha}=f_{\alpha}\vec{v}$ will point in the same direction in the
absence of gravity and capillary effects, and the upwind definition is
explicit. However, referring back to (4), we see that the upwind definition is
generally implicit in the presence of gravity/capillary effects, since
quantities in the parenthesis depend on the mobilities $\lambda_{\alpha}$. An
equivalent explicit upwind definition has been derived by Brenier and Jaffré
[1991] for finite-volume methods. Herein, we employ this definition at each
cubature point involved in evaluating the surface integrals.

With this notation, we define the following weak form of the residual equation
(5) in cell $i$

$\displaystyle\mathcal{R}_{w,i}(S_{w},S_{o},S_{g},\psi)$ (6)
$\displaystyle=\tfrac{1}{\Delta
t}\mathcal{A}_{w,i}(S_{w},\psi)+\mathcal{F}_{w,i}(S_{w},S_{o},S_{g},\psi)-\mathcal{Q}_{w,i}(\psi),$

where the accumulation, flux and source/sink functionals are defined as

$\displaystyle\mathcal{A}_{w,i}(S_{w},\psi)=$
$\displaystyle\int_{\Omega_{i}}\Big{(}(\phi b_{w}S_{w})-(\phi
b_{w}S_{w})^{n}\Big{)}\psi\,\mathrm{d}V,$
$\displaystyle\mathcal{F}_{w,i}(S_{w},S_{g},S$ ${}_{o},\psi)$
$\displaystyle=\sum_{j\in N(i)}$
$\displaystyle\int_{\Gamma_{ij}}u\left[b_{w},f_{w}(S_{w},S_{o},S_{g})\vec{v}\cdot\vec{n}_{ij}\right]\psi\,\mathrm{d}\sigma$
$\displaystyle-$
$\displaystyle\int_{\Omega_{i}}b_{w}f_{w}(S_{w},S_{o},S_{g})\vec{v}\cdot\nabla\psi\,\mathrm{d}V,$
$\displaystyle\mathcal{Q}_{w,i}(\psi)=$
$\displaystyle\int_{\Omega_{i}}b_{w}q_{w}\psi\,\mathrm{d}V.$

Weak form residuals for the oil ($\mathcal{R}_{o,i}$) and gas
($\mathcal{R}_{g,i}$) can be derived in a similar fashion. To obtain a
discrete weak formulation, we replace the function space $V$ with a finite-
dimensional subspace $V_{h}$ consisting of functions that are smooth on each
cell $\Omega_{i}$, but possibly discontinuous across cell interfaces. We
replace the saturations and the test function $\psi$ with approximations
$S_{\alpha,h}\in V_{h}$ and $\psi_{h}\in V_{h}$, and arrive at the following
discrete weak formulation:

$\displaystyle\text{Find }S_{w,h},$ $\displaystyle S_{o,h},S_{g,h}\in
V_{h}\text{ such that }$ $\displaystyle\mathcal{R}_{\alpha,i}($ $\displaystyle
S_{w,h},S_{o,h},S_{g,h},\psi_{h})=0\text{ for all }\psi_{h}\in V_{h}.$

If we introduce the basis $\\{\psi_{k}\\}_{k=1}^{n_{\mathrm{dof}}}$ for
$V_{h}$, we may express the saturations as
$S_{\alpha,h}=\sum_{k=1}^{n_{\mathrm{dof}}}S_{\alpha,k}\psi_{k}$, where
$S_{\alpha,k}\in$ is referred to as the $k$th degree of freedom of
$S_{\alpha}$. The discrete weak formulation then takes the form

$\displaystyle\text{Find
}(S_{\alpha,1},\dots,S_{\alpha,n_{\mathrm{dof}}})\in^{n_{\mathrm{dof}}},\,\alpha=w,o,g,\text{
such that }$
$\displaystyle\,\mathcal{R}_{\alpha,i}\left(\sum_{k=1}^{n_{\mathrm{dof}}}S_{w,k}\psi_{k},\sum_{k=1}^{n_{\mathrm{dof}}}S_{o,k}\psi_{k},\sum_{k=1}^{n_{\mathrm{dof}}}S_{g,k}\psi_{k},\psi_{\ell}\right)=0$
$\displaystyle\,\text{for all }\ell=1,\dots,n_{\mathrm{dof}}.$

Note that we typically use the closure relation $S_{w}+S_{o}+S_{g}=1$ to
eliminate one variable and one equation from the above.

There are several important factors that need careful consideration to
implement the above discretization. In the following, we briefly discuss the
most important of these. A more detailed analysis is subject of further
research.

###  4.1 Basis functions

In finite-element methods, it is common to use orthogonal basis functions
$\\{\psi_{k}\\}_{k=1}^{n_{\mathrm{dof}}}$, i.e.,
$\int_{\Omega_{i}}\psi_{k}\psi_{\ell}=\delta_{k,\ell}$. The main reason for
this choice is that it typically results in a significantly less dense
discretization matrix. Defining orthogonal basis functions for general
polyhedral grids is not a trivial problem. For the types of problems we are
interested in, however, the weak form residual (6) will generally depend
nonlinearly on the saturations (and thus the basis functions) due to the
fractional flow function $f_{\alpha}$ in the flux terms
$\mathcal{F}_{\alpha,i}$. Effectively, orthogonal basis functions will
generally not lead to less dense discretization matrices. Herein, we will
therefore simply use tensor products of Legendre polynomials as our basis
functions. These form a convenient basis for the space of polynomials and are
orthogonal for cuboid grid cells. The first two Legendre basis functions are
$\ell_{0}(x)=1$ and $\ell_{1}(x)=x$, and higher-order basis functions are
defined in our dG implementation by exploiting Bonnet’s recursion formula,

$(k+1)\ell_{k+1}(x)-(2k+1)x\ell_{k}(x)+k\ell_{k-1}(x)=0.$

Basis function $j$ of polynomial order $r+s+t$ for cell $i$ is then defined as

$\psi_{j}(\vec{x})=\begin{cases}\ell_{r}\left(\frac{x-x_{i}}{\Delta
x_{i}/2}\right)\ell_{s}\left(\frac{y-y_{i}}{\Delta
y_{i}/2}\right)\ell_{t}\left(\frac{z-z_{i}}{\Delta
z_{i}/2}\right),&\vec{x}\in\Omega_{i},\\\
0,&\vec{x}\not\in\Omega_{i},\end{cases}$

where

$\vec{x}=(x,y,z),\quad\vec{x}_{i}=(x_{i},y_{i},z_{i}),\quad\Delta\vec{x}=(\Delta
x_{i},\Delta y_{i},\Delta z_{i})$

denote the spatial coordinate, cell centroid and the dimensions of the
smallest cuboid aligned with the coordinate axes that completely contains cell
$i$, as shown in Figure 1. For a dG method of degree $k$, the basis consists
of all tensor products on this form such that $0\leq r+s+t\leq k$, with
$r,s,t\geq 0$. In $d$ space dimensions, this gives a total number of basis
functions

$n_{\mathrm{dof}}=\binom{k+d}{d}=\frac{(k+d)!}{d!k!}.$

We use the notation dG($k$) to refer to a dG method of degree $k$. The error
of a method of degree $k$ will for smooth solutions decay with order $k+1$, so
that a dG($k$) method is of formal order $k+1$.

Figure 1: Dimensions of bounding boxes used to define the dG basis functions
for a polygon in 2D and a polyhedron in 3D.

###  4.2 Cubature rules

The integrals in (6) must be numerically evaluated. Efficient cubature rules
for general polyhedral grid cells are crucial to obtain an efficient
implementation. A straightforward approach to numerically evaluate integrals
over polyhedral cells is to partition the cell into non-overlapping simplices
and apply a known cubature rule for these. This uses significantly more
cubature points than strictly needed for the cubature to be correct and is
thus highly inefficient from a computational point of view. Instead, we apply
an approach based on moment-fitting, see e.g., Müller et al. [2013]. This
approach defines the quadrature rule for a cell $\Omega$ (omitting subscript
$i$) as the solution to a small linear system

$\displaystyle\begin{bmatrix}\psi_{1}(\vec{x}_{1})&\dots&\psi_{1}(\vec{x}_{n_{\mathrm{dof}}})\\\
\vdots&\ddots&\vdots\\\
\psi_{n_{\mathrm{dof}}}(\vec{x}_{1})&\dots&\psi_{n_{\mathrm{dof}}}(\vec{x}_{n_{\mathrm{dof}}})\\\
\end{bmatrix}\begin{bmatrix}w_{1}\\\ \vdots\\\
w_{n_{\mathrm{dof}}}\end{bmatrix}$ $\displaystyle\hskip
60.00009pt=\frac{1}{|\Omega|}\begin{bmatrix}\int_{\Omega}\psi_{1}\,\mathrm{d}V\\\
\vdots\\\
\int_{\Omega}\psi_{n_{\mathrm{dof}}}\,\mathrm{d}V\end{bmatrix},\quad\text{or}\quad\Psi\vec{w}=\vec{b}.$

The integrals (or moments) on the right-hand side can be evaluated using a
sub-optimal cubature rule. If we pick the cubature points
$\vec{x}_{0},\dots\vec{x}_{n_{\mathrm{dof}}}$ so that the matrix $\Psi$ is
invertible, we can solve for the quadrature weights $w_{1},\dots
w_{n_{\mathrm{dof}}}$ to obtain a cubature rule with $n_{\mathrm{dof}}$ points
for each cell. Note that it is possible to construct a quadrature rule of
precision $k$ with less than $n_{\mathrm{dof}}$ points by eliminating points
with marginal significance. This computation can be done once for each grid
cell in a preprocessing step.

In 3D, the surface integrals in the second term of $\mathcal{F}$ in (6) are
evaluated in the same way as areal integrals in 2D, using the moment-fitting
approach described above. In 2D, the surface integrals are simple line
integrals, for which a standard Gauss cubature rule is employed.

###  4.3 Velocity interpolation

As stated above, we assume that the solution to the pressure equation yields a
total velocity as a set of fluxes that are constant on each interface. We see
from (6) that for linear and higher-order basis functions, we need information
about the velocity $\vec{v}$ also inside the cell, meaning that we must
interpolate from the interface fluxes. Herein, we apply a simple scheme
inspired by the mimetic finite difference method [Lie et al., 2012a]. This
gives a constant velocity inside the cell, consistent with the volumetric
interface fluxes $v_{ij}$:

$\vec{v}_{i}=\frac{1}{|\Omega_{i}|}\sum_{j\in
N(i)}v_{ij}(\vec{x}_{ij}-\vec{x}_{i}).$

Here, $\vec{x}_{i}$ and $\vec{x}_{ij}$ refer to the centroid of cell $i$ and
interface $ij$, respectively. The velocity on an interface is assumed to be
constant. More sophisticated interpolation schemes like the extended corner-
velocity interpolation scheme (ECVI) [Klausen et al., 2012] have earlier been
investigated for dG methods for flow diagnostics [Rasmussen and Lie, 2014],
but will for simplicity not be applied herein.

###  4.4 Order reduction

Higher-order methods usually require a slope limiter, or some other means to
prevent the creation of over- and under-shoots or other types of nonphysical
oscillations near spatial discontinuities in the solution. Natvig et al.
[2007] previously showed that order reduction followed by local and dynamic
grid refinement could be used as an alternative strategy to avoid oscillations
and give high-resolution in a robust and cost-efficient manner. Herein, we use
a simplified version: Whenever the jump across an intercell interface exceeds
a prescribed tolerance, we reduce the local approximation space to dG(0). The
same is done when the solution exceeds its domain of definition by a small
factor $\varepsilon$.

##  5 Localized nonlinear solvers

This section discusses how we can utilize certain unidirectional flow
properties to formulate efficient nonlinear solution strategies that localize
the Newton procedure to avoid spending unnecessary iterations in regions of
the reservoir where many iterations are not needed.

###  5.1 Reordering based on intercell fluxes

An important physical property of the transport (1) is that the flow is
unidirectional along streamlines if no gravity or capillary effects are
present. Each dG basis function is restricted to a single cell only, and hence
the only intercell interactions in the weak form $\mathcal{R}_{w,i}$ in (6)
are due to the upwind operator. We split the neighboring cells
$\mathcal{N}(i)$ of $i$ into upstream cells $\mathcal{U}(i)$ and downstream
cells $\mathcal{D}(i)$ based on the total flux $\vec{v}\cdot\vec{n}_{ij}$:

$\displaystyle\sum_{j\in\mathcal{N}(i)}$
$\displaystyle\int_{\Gamma_{ij}}u\left[b,f\vec{v}\cdot\vec{n}_{ij}\right]\psi\,\mathrm{d}\sigma$
$\displaystyle=\sum_{j\in\mathcal{U}(i)}\int_{\Gamma_{ij}}(bf)_{j}\,\vec{v}\cdot\vec{n}_{ij}\psi\,\mathrm{d}\sigma$
$\displaystyle(\text{Upstream})$
$\displaystyle+\sum_{j\in\mathcal{D}(i)}\int_{\Gamma_{ij}}(bf)_{i}\,\vec{v}\cdot\vec{n}_{ij}\psi\,\mathrm{d}\sigma,$
$\displaystyle(\text{Downstream})$

where $(\cdot)_{i}$ indicates that quantities inside the parentheses should be
evaluated from cell $i$. Recall that all other unknowns in the weak-form
residual $\mathcal{R}_{w,i}$ have support limited to cell $i$, and hence the
inherent unidirectional property of the continuous equations (1) is preserved
by the discretization. In practice, this means that if we know the solution in
all upstream cells $\mathcal{U}(i)$, the only unknowns in $\mathcal{R}_{w,i}$
are the ones associated with cell $i$. That is, if we now perform a
topological sort of the directed, acyclic graph (DAG) induced by the intercell
fluxes and use this to reorder the cells, we can solve the transport equations
cell-by-cell by traversing the sorted graph. The algebraic interpretation of
this is the following: If we linearize the nonlinear transport equations for
all cells simultaneously, and permute the system according to the topological
order, we obtain a lower-triangular matrix. Note, however, that we never
assemble the discretization matrix for the full system in the reordering
solution procedure. Instead we solve the nonlinear transport equations
$\mathcal{R}_{\alpha,i}=0$ cell-by-cell. This way, we avoid expensive
linearizations of large systems of nonlinear equations.

Figure 2 illustrates the reordering principle applied to a dG(0) scheme for a
quarter-five spot test case formulated on a 2D Voronoi grid. In the original
grid, the cells are ordered almost, but not entirely, by the $x$-coordinate of
their centroids, giving an irregular sparsity pattern with entries both above
and below the diagonal. After reordering, the cells are ordered such that if
$j\in\mathcal{N}(i)$, then $j\in\mathcal{U}(i)$ for $j<i$ and
$j\in\mathcal{D}(i)$ for $j>i$, where we recall that $\mathcal{N}(i)$ denotes
neighbors and $\mathcal{U}(i)$ and $\mathcal{D}(i)$ are cells that lie
upstream/downstream of cell $i$.

Figure 2: Sparsity pattern for a dG(0) method formulated on a Voronoi grid
with original (left) and reordered (right) numbering of cells. The inset shows
grid cell 35 and its upstream $\mathcal{U}(35)$ and downstream
$\mathcal{D}(35)$ neighbors.

Figure 3 similarly illustrates the difference between a dG(0) and a dG(1)
discretization on a Cartesian grid. The two matrices have exactly the same
flux graph, and thus obtain the same cell ordering. However, whereas dG(0)
gives a nonlinear scalar problem in each cell, dG(1) gives a $3\times 3$
nonlinear system.

Figure 3: Sparsity pattern after reordering for dG(0)/dG(1) for a five-spot
problem on a rectangular grid. Figure 4: Examples of different degrees of
coupling, where the cycles can be seen as blocks on the diagonal of the dG(0)
discretization matrix. In the first sparsity pattern, the system is almost
completely reorderable with only two small cycles; in the second, a small
fraction of the cells belongs to four small cycles; the third patterns shows
an example where almost all cells belongs to eight larger cycles; finally the
last pattern shows a case where all cells belong to one huge cycle. All cells
in the same cycle must be solved simultaneously. Figure from Lie et al.
[2014].

When gravity and capillary effects are included, the flow will no longer be
unidirectional along streamlines. As a result, the intercell flux graph may
contain cycles, so that the reordering method cannot be applied directly.
However, we may construct a DAG from the intercell flux graph by grouping
cells that are part of the same cycle into a single supernode. By
topologically sorting this DAG, we end up with a block-triangular system in
which each block consists of degrees of freedom that are mutually coupled and
must be solved for simultaneously. Note that as explained above, we evaluate
the upstream direction for each quadrature point at the interface. For a
higher-order method, the upstream direction may thus vary along a single
interface when gravity/capillary effects are included, particularly when we
use high polynomial order and/or cells with large vertical extent. This means
that each quadrature point on interfaces with changing upstream direction must
be treated as an edge in the original intercell flux graph, and such
interfaces will automatically introduce a local cycle and require the
neighboring cells to be grouped into a supernode. In general, the reordering
procedure will result in a DAG with $n_{c}$ nodes, where $n_{c}$ equals the
number of connected components of the original intercell flux graph.
Topologically sorting this graph gives us an ordering

$\mathcal{P}=\left(\mathcal{P}_{1},\dots,\mathcal{P}_{n_{c}}\right),\quad\mathcal{P}_{k}=(i_{1},\dots,i_{n_{k}}),$

where each $\mathcal{P}_{k}$ is a set of cell indices, and $n_{k}$ is the
number of cells in connected component $k$. The number of connected components
$n_{c}$ may vary from one in the worst case, where all cells are connected in
a single cycle, to the total number of cells $N$ in the best case, where the
problem can be solved cell by cell without the need of any supernodes. Given
two connected components $\mathcal{P}_{k}$ and $\mathcal{P}_{\ell}$, we will
then have that if $i>j$ for a cell $i\in\mathcal{P}_{k}$ and a cell
$j\in\mathcal{P}_{\ell}$ then $i>j$ for all cells $i\in\mathcal{P}_{k}$,
$j\in\mathcal{P}_{\ell}$. Figure 4 illustrates different scenarios of
coupling, from an almost completely reordered system, to a system where all
cells belong to a single cycle.

###  5.2 Block-wise processing

For test cases with only viscous forces, the reordered system can be solved
cell-by-cell in a single pass. This is optimal in the sense that it minimizes
the number of cell-wise nonlinear solves over the simulation. However, it is
not necessarily optimal in terms of runtime; most modern workstations contain
several CPU cores, as well as any number of SIMD-type registers. In this
setting, feeding sufficient computations to the CPU is often a large
bottleneck and a barrier to rapid execution time. For this reason, our solver
has an optional block-wise algorithm that solves a prescribed number of cells
simultaneously. If the largest ordering index of the previous solve was $i$,
this algorithm then solves for indices $i+1$ to $i+n_{b}$, where $n_{b}$ is
the block size. This gives a larger system, and individual cells can follow a
sub-optimal solution path when solved by a multi-cell Newton–Raphson solver.
On the other hand, the startup cost is incurred only once per block of cells
and not for every cell, and the block solve utilizes cache better and leaves
less computational resources idle. Hence, simultaneous solve of $n_{b}$ cells
will typically be significantly faster than $n_{b}$ single-cell solvers.

Figure 5: Schematic overview of the solution procedure. The pressure equation
is solved with a suitable solver, resulting in intercell fluxes $v_{ij}$,
which we may interpret as a directed graph with grid cells as nodes, and
intercell fluxes as edge weights. By grouping all cells that are part of the
same cycle into a (super)node, we obtain a directed, acyclic graph (DAG). A
topological sort of this graph gives us a node (or component) ordering
$\mathcal{P}$. By traversing the nodes in this order, we can solve the
transport problem $\mathcal{R}_{t}(S)=0$ node by node, since the cells in each
node in the DAG only depend on each other and cells from upstream (and already
resolved) nodes.

###  5.3 Nonlinear Gauss–Seidel solvers

As explained above, gravity will result in a block-triangular system in which
each block consists of degrees of freedom that are mutually coupled and must
be solved for simultaneously. The size of each block depends on whether the
block represents degrees of freedom that are part of a co-current or counter-
current flow region. Blocks in co-current regions represent the discrete
system posed inside a single cell and are treated as described above. Blocks
in counter-current regions consist of multiple cells that are coupled in
cycles. In principle, each of these cycles of cells should be solved by a
localized Newton iteration. Another possibility is to first solve the water
equations cell-by-cell in the order of descending water phase potential, and
then solve the oil equations in the order of descending oil phase potential
Kwok and Tchelepi [2007]. This method requires that the oil component equation
only depends on water and oil saturation, which is the case for black-oil
models. Herein, we use a nonlinear Gauss–Seidel strategy instead that simply
loops through all the cells of a cycle and solves them one at a time, using
the values in all neighboring cells as if they were all known upstream values.
This will typically not give the correct solution after a single iteration
through the loop, but a large number of numerical experiments indicate that
the process converges after only a few repetitions. Capillary forces will in
the worst case couple all cells in one big cycle, as seen in the rightmost
sparsity pattern in Figure 4. However, the coupling is oftentimes so weak that
a Gauss–Seidel type iteration is still very effective.

##  6 Solution procedure

The entire procedure to evolve the solution one time step can now be
summarized as follows: (i) Solve the pressure equation (3) to obtain intercell
fluxes $v_{ij}$ and oil phase pressure $p_{o}$; (ii) construct a DAG from the
intercell flux graph by grouping cells that are part of the same cycle into a
single supernode; (iii) starting from the first node in the sorted graph,
solve the transport equations (6) cell-by-cell or cycle-by-cycle in
topological order. The procedure is summarized in Figure 5. Notice the
important property that if the residual of a cell or cycle is converged at the
beginning of the time step, and none of the upstream neighbors were updated in
this step, we do not need to call the nonlinear Gauss–Seidel solver and can
proceed to the next component. Effectively, we thus avoid solving for a large
number of zeros, and can instead focus our computational resources on parts of
the domain where updates are nonzero.

Figure 6: Simulation responses for two of the many wells in the Norne field
model. We plot the bottom-hole-pressure and production rates for water, oil
and gas at surface conditions.

##  7 Numerical examples

We have implemented our reordering approach in two different open-source
simulator platforms. The MATLAB Reservoir Simulation Toolbox (MRST) [Lie,
2016] is a general toolbox aimed at rapid development of proof-of-concept
simulators and workflow tools. We have used MRST to develop a combined
dG–reordering method that allows for cell-based refinement. The solver uses
automatic differentiation to linearize the localized systems and is as such
applicable to a wide variety of multiphase flow models. At the time of
writing, the implementation in MRST has only been tested for two-phase models
with varying degrees of compressibility.

OPM Flow [OPM, 2019] is developed with the goal of providing industry-grade
simulations and generally offers single-core computational performance that is
comparable with contemporary commercial simulators. We have used this as a
platform to implement dG(0) with reordering for a general three-phase black-
oil model with compressibility, dissolution and vaporization, gravity and
capillary forces, and hysteresis. Our implementation uses a simple, global
nonlinear Gauss–Seidel method that traverses all cells repeatedly until the
residual of the transport equations are below a prescribed tolerance in all
cells.

###  7.1 The Norne oil field

Figure 7: Initial fluid saturations and wells (left); aggregate statistics for
the number of cycles (right) for the Norne field model in Example 7.1.

The Norne oil field in the Norwegian Sea is one of the few real assets for
which simulation models and other data have been made openly available The
Open Porous Media Initative(2015) [OPM]. Over the 3,312 days of production
history, 36 wells are active. We study a water-alternating-gas (WAG) injection
scenario, in which the wells are controlled by reservoir-volume rates to match
the historical rates to obtain a reasonable pressure development. The purpose
of the example is to demonstrate that the sequential splitting method is
capable of capturing the correct behavior of an industry-grade simulation run
in history-matching mode.

We ran the case with two simulators from OPM: the fully implicit Flow
simulator and the experimental Flow-reorder simulator, which uses sequential
splitting and reordering of the transport equations. We note that these
simulators both support a slightly more complex fluid model than that
described in the introduction, incorporating hysteresis effects for relative
permeability. Both simulators are available from OPM in source and binary form
[OPM, 2019].

We highlight the behavior of two different wells to compare the fully implicit
to the sequential/reordering solution. Well “D-3AH” is a typical producer and
shows good match between the fully implicit and reordering solutions for both
bottom-hole pressure and rates; see the upper row in Figure 6. Well “B-2H” is
also a producer, which we have chosen to highlight because it shows the worst
match between the two solutions. As can be seen in the lower row of Figure 6,
gas and oil production rates match well, but water-production rate and bottom-
hole pressure show significant deviation some time after water breakthrough.
There are two main factors that contribute to the difference: First of all,
the sequential splitting in Flow-reorder does not use outer iterations. The
overall solution algorithm is thus not guaranteed to reduce the overall
residual at the end of each time step to the same tolerance as the fully
implicit method in Flow. Moreover, convergence is measured somewhat
differently in the global Gauss–Seidel and the Newton–Raphson methods.

In our opinion, this simulation is proof of concept for the reordering idea in
a production-grade code running a full simulation of a real field in history
mode. This is a first step, and there are several ways in which this
implementation can be improved. The important observation is that we made no
simplifications in the description of geology, fluid behavior, and well and
production facilities.

Figure 7 shows the fluids in place and reports the occurrence of cycles
throughout the simulation. Altogether, the model has approximately 44,000
active cells. The number of cycles and the average number of cells in each
cycle vary largely throughout the simulation because of significant changes in
well controls and fluid composition. In any time step, the number of cells
belonging to cycles only represents a fraction of the total cells. This shows
the inherent locality and the speedup potential, since most of the cells can
be solved independently.

###  7.2 Subset of the SPE10 Model 2

Figure 8: Petrophysical properties for Layer 50 from the SPE10 Benchmark. The
right plot shows the cell index after reordering.

Figure 9: Evolution of the water saturation for three different times in the
channelized reservoir in Example 7.2. The upper row is for dG(0) and the lower
row for dG(1).

Figure 10: Number of nonlinear iterations per cell at three instances in time
corresponding to the solution profiles shown in Figure 11. Zero iterations
were performed for cells in white. The upper row is for dG(0) and the lower
row for dG(1).

Figure 11: Comparison of average number of nonlinear iterations per cell as
function of time for Example 7.2. Note that the $y$-axis is logarithmic.
Figure 12: Petrophysical properties and well setup for the Voronoi-example.

Next, we consider a 2D simulation model with petrophysical properties sampled
from Layer 50 from the SPE10 Model 2 Benchmark [Christie and Blunt, 2001].
This layer is part of a fluvial formation and consists of an intertwined
pattern of high-permeable sand channels on a background of low-permeable
mudstone. We inject 0.2 pore volumes of water at a constant rate in one end of
a channel, and produce at a fixed bottom-hole pressure of 275 bars at the
opposite end of the same channel. The reservoir is initially filled with a
mixture of water and oil, and the rock and fluids are weakly compressible.
Water and oil have quadratic relative permeabilities with residual saturations
of 0.2, and viscosities of 2.85 and 3.0 cP, respectively. Permeability and
porosity are depicted in Figure 8. The figure also shows the cell index after
reordering. As expected, the cell number increases from the injector and along
the high-flow channels to the producer. Interestingly, the model contains
cycles of cells that are not connected to the injector–producer pair by
nonzero fluxes. This is an effect of fluid compressibility, and the reordering
algorithm assigns lower indices to cells which essentially should remain
inactive throughout the transport time step.

We simulate using dG(0) and dG(1) with reordering implemented in MRST. Figure
11 shows the resulting water saturation at selected time steps, whereas Figure
11 shows the iterations used per active cell in the corresponding time steps.
By active cells, we mean cells in which the water saturation was updated
during the time step either because the residual at the beginning of the time
step was above the nonlinear tolerance, or because one or more of the cell’s
upstream neighbors were updated during the time step. Notice that iterations
are performed only in a small fraction of the cells in each time step and that
the updates are chiefly located to the fluvial channel. We also see that the
higher-order solution profile is much less diffuse. Effectively, the water
channel fills faster, so that dG(1) predicts earlier water breakthrough than
dG(0).

For comparison, Figure 11 also reports the (average) number of nonlinear
Newton–Raphson iterations used to solve each time step. For the reordered
solver, we report the average number of nonlinear iterations used per cell,
defined over all cells, and the maximum number of iterations observed in the
cell with slowest nonlinear convergence. Because large fractions of the cells
remain inactive in most of the time steps, the average number of cellwise
nonlinear iterations is far below one in all time steps. This confirms
previous observations by Natvig and Lie [2008]. Interestingly, we see that the
maximum number of nonlinear iterations for the reordered solver in most time
steps is less than the number of Newton–Raphson solves used when solving for
all cells simultaneously for the same time step. This indicates that the
iteration path followed by individual cells in a Newton–Raphson solve of all
cells simultaneously is not necessarily optimal.

Figure 13: Thresholded water saturation profiles for the two dG solutions at
three selected time steps for Example 7.3.

Figure 14: Average number of cellwise nonlinear iterations per time step used
for each of the solvers in Example 7.3, with dG(0) in the left plot and dG(1)
in the right plot.

###  7.3 3D Voronoi grid

In this example, we study the effect of using higher-order dG methods and
reordering to simulate a case posed on a 3D, fully unstructured PEBI/Voronoi
grid generated using the upr module in MRST [Berge et al., 2019, Klemetsdal et
al., 2017]. The domain is comprised of a rectangular box of dimensions
$300\times 100\times 100$ m3 with approximately 20, 10, and 10 cells in each
of the axial directions, respectively. Porosity and permeability are sampled
from the top layers of Model 2 from SPE10 [Christie and Blunt, 2001]. The
injector is placed in a cell at the center of the left boundary. We inject
water at constant bottom-hole pressure of 600 bars. A producer operated at a
fixed bottom-hole pressure of 50 bars is placed at the center of the right
boundary. Figure 12 shows petrophysical properties and well positions.

We simulate the case using the dG(0) and dG(1) solvers from MRST with and
without reordering. For the reordering algorithm, we use the block version
with a block size $n_{b}$ of 100 cells. Figure 14 reports saturation profiles
at three selected time steps. Compared with dG(1), the dG(0) saturation
profiles are more diffusive and fill more of the reservoir cross-section. As a
result, dG(1) predicts a higher water-cut in the producer, as reported in
Figure 15.

The solutions computed with the ordered and the original Newton–Raphson
solvers are identical to plotting resolution, which verifies the correctness
of the localization and the new block algorithm. Finally, we look at the
number of iterations used by the different solvers, shown in Figure 14. We
observe that even with a block size of $n_{b}=100$ cells, the average number
of nonlinear iterations per cell is significantly smaller than the number of
iterations used when solving for all cells simultaneously.

Figure 15: Water cut in the production well in Example 7.3.

##  8 Concluding remarks

The transport substep in a sequentially implicit solution procedure for black-
oil models can be significantly accelerated by ordering the grid cells in the
simulation model according to total velocity computed in a proceeding pressure
step, so that the residual transport equations can be solved cell-by-cell from
inflow to outflow. This not only localizes the nonlinear Newton solver, so
that the computational effort can be focused on cells where the transported
quantities change significantly during each time step, but also contributes to
improve the general robustness of the Newton solver. As a verification of this
concept in an industry-grade setting, we have used a sequential simulator from
the open-source OPM initiative to demonstrate that a reordered, global
nonlinear Gauss–Seidel strategy computes an acceptable solution for the full
simulation model of the Norne oil field. Additional verifications have also
been obtained using similar industry-grade simulators from MRST.

The ordering approach is also applicable to transport solvers based on
discontinuous Galerkin discretizations in combination with single-point
upstream mobility weighting. Such discretizations are detailed herein for
general 3D cases with (almost) no assumptions on the polyhedral cell
geometries and grid topology. By reordering cells optimally, one can ensure
that the added computational cost of solving for more degrees of freedom to
obtain higher spatial order remains local to each cell. Our results obtained
using MRST show that a second-order dG(1) method has significantly improved
accuracy compared with the standard first-order volume method used in industry
at the expense of only a moderate increase in the number of nonlinear
iterations.

Herein, we have only presented simulations of two-phase oil-water and three-
phase black-oil systems, but the same principles are applicable to a wide
variety of different black-oil type and compositional models used to described
secondary and tertiary recovery processes. As an example, Klemetsdal et al.
[2019b] discuss several cases with significantly more complex flow physics,
including compositional flow and gas injection with strong gravity effects, as
well as use of dynamically adaptive grid coarsening.

##  9 Acknowledgements

The work of Klemetsdal, Rasmussen, and Lie is funded by the Research Council
of Norway through grant no. 244361. Møyner is funded by VISTA, a basic
research program funded by Equinor and conducted in close collaboration with
the Norwegian Academy of Science and Letters.

## References

  * Appleyard and Cheshire [1982] J. R. Appleyard and I. M. Cheshire.  The cascade method for accelerated convergence in implicit simulators.  In _European Petroleum Conference, 25-28 October, London, United Kingdom_. Society of Petroleum Engineers, October 1982.  doi: 10.2118/12804-MS.
  * Bell et al. [1986] J. B. Bell, J. A. Trangenstein, and G. R. Shubin.  Conservation laws of mixed type describing three-phase flow in porous media.  _SIAM Journ. Appl. Math._ , 46(6):1000–1017, 1986\.  doi: 10.1137/0146059.
  * Berge et al. [2019] R. L. Berge, Ø. S. Klemetsdal, and K.-A. Lie.  Unstructured voronoi grids conforming to lower dimensional objects.  _Comput. Geosci._ , 23(1):169–188, 2019.  doi: 10.1007/s10596-018-9790-0.
  * Bratvedt et al. [1996] F. Bratvedt, T. Gimse, and C. Tegnander.  Streamline computations for porous media flow including gravity.  _Transp. Porous Media_ , 25(1):63–78, 1996.  doi: 10.1007/BF00141262.
  * Brenier and Jaffré [1991] Y. Brenier and J. Jaffré.  Upstream differencing for multiphase flow in reservoir simulation.  _SIAM J. Numer. Anal._ , 28(3):685–696, 1991\.  doi: 10.1137/0728036.
  * Christie and Blunt [2001] M. A. Christie and M. J. Blunt.  Tenth SPE comparative solution project: A comparison of upscaling techniques.  _SPE Reservoir Eval. Eng._ , 4:308–317, 2001.  doi: 10.2118/72469-PA.
  * Cockburn and Shu [1989] B. Cockburn and C.-W. Shu.  TVB runge-kutta local projection discontinuous galerkin finite element method for conservation laws ii: General.  _Math. Comp._ , 52(186):411–435, 1989.  doi: 10.1090/S0025-5718-1989-0983311-4.
  * Cockburn and Shu [1991] B. Cockburn and C.-W. Shu.  The runge-kutta local projection $p^{1}$\- discontinuous-galerkin finite element method for scalar conservation laws.  _ESAIM–Math. Model. Num._ , 25(3):337–361, 1991\.
  * Datta-Gupta and King [2007] A. Datta-Gupta and M. J. King.  _Streamline simulation: theory and practice_ , volume 11 of _SPE Textbook Series_.  Society of Petroleum Engineers, 2007.
  * Eikemo et al. [2009] B. Eikemo, K.-A. Lie, H. K. Dahle, and G. T. Eigestad.  Discontinuous Galerkin methods for transport in advective transport in single-continuum models of fractured media.  _Adv. Water Resour._ , 32(4):493–506, 2009.  doi: 10.1016/j.advwatres.2008.12.010.
  * Gries et al. [2014] S. Gries, K. Stüben, G. L. Brown, D. Chen, and D. A. Collins.  Preconditioning for efficiently applying algebraic multigrid in fully implicit reservoir simulations.  _SPE J._ , 19(04):726–736, 2014.  doi: 10.2118/163608-PA.
  * Jenny et al. [2006] P. Jenny, S. H. Lee, and H. A. Tchelepi.  Adaptive fully implicit multi-scale finite-volume method for multi-phase flow and transport in heterogeneous porous media.  _J. Comput. Phys._ , 217(2):627–641, 2006.  doi: 10.1016/j.jcp.2006.01.028.
  * Klausen et al. [2012] R. A. Klausen, A. F. Rasmussen, and A. Stephansen.  Velocity interpolation and streamline tracing on irregular geometries.  _Comput. Geosci._ , 16:261–276, 2012.  doi: 10.1007/s10596-011-9256-0.
  * Klemetsdal et al. [2017] Ø. S. Klemetsdal, R. L. Berge, K.-A. Lie, H. M. Nilsen, and O. Møyner.  Unstructured gridding and consistent discretizations for reservoirs with faults and complex wells.  In _SPE Reservoir Simulation Conference_. Society of Petroleum Engineers, 2017.
  * Klemetsdal et al. [2019a] Ø. S. Klemetsdal, O. Møyner, and K.-A. Lie.  Robust nonlinear newton solver with adaptive interface-localized trust regions.  _SPE J._ , 2019a.  doi: 10.2118/195682-PA.
  * Klemetsdal et al. [2019b] Ø. S. Klemetsdal, O. Møyner, and K.-A. Lie.  Implicit high-resolution compositional simulation with optimal ordering of unknowns and adaptive spatial refinement.  In _SPE Reservoir Simulation Conference, 10-11 April, Galveston, Texas, USA_. Society of Petroleum Engineers, 2019b.  doi: 10.2118/193934-MS.
  * Kwok and Tchelepi [2007] F. Kwok and H. Tchelepi.  Potential-based reduced newton algorithm for nonlinear multiphase flow in porous media.  _J. Comput. Phys._ , 227(1):706 – 727, 2007.  doi: 10.1016/j.jcp.2007.08.012.
  * Lie [2016] K.-A. Lie.  _An Introduction to Reservoir Simulation Using MATLAB/GNU Octave: User guide for the Matlab Reservoir Simulation Toolbox (MRST)_.  Cambridge University Press, 2016.
  * Lie et al. [2012a] K.-A. Lie, S. Krogstad, I. S. Ligaarden, J. R. Natvig, H. M. Nilsen, and B. Skaflestad.  Open source MATLAB implementation of consistent discretisations on complex grids.  _Comput. Geosci._ , 16:297–322, 2012a.  ISSN 1420-0597.  doi: 10.1007/s10596-011-9244-4.
  * Lie et al. [2012b] K. A. Lie, J. R. Natvig, and H. M. Nilsen.  Discussion of dynamics and operator splitting techniques for two-phase flow with gravity.  _Int. J. Numer. Anal. Mod._ , 9(3):684–700, 2012b.
  * Lie et al. [2014] K.-A. Lie, H. M. Nilsen, A. F. Rasmussen, and X. Raynaud.  Fast simulation of polymer injection in heavy-oil reservoirs on the basis of topological sorting and sequential splitting.  _SPE J._ , 19(06):0991–1004, 2014.  doi: 10.2118/163599-PA.
  * Lie et al. [2017] K.-A. Lie, O. Møyner, J. R. Natvig, A. Kozlova, K. Bratvedt, S. Watanabe, and Z. Li.  Successful application of multiscale methods in a real reservoir simulator environment.  _Comput. Geosci._ , 21(5):981–998, Dec 2017.  doi: 10.1007/s10596-017-9627-2.
  * Møyner [2017] O. Møyner.  Nonlinear solver for three-phase transport problems based on approximate trust regions.  _Comput. Geosci._ , 21(5-6):999–1021, 2017.  doi: 10.1007/s10596-017-9660-1.
  * Müller et al. [2013] B. Müller, F. Kummer, and M. Oberlack.  Highly accurate surface and volume integration on implicit domains by means of moment-fitting.  _Int. J. Numer. Meth. Eng._ , 96(8):512–528, 2013\.  doi: 10.1002/nme.4569.
  * Natvig and Lie [2008] J. R. Natvig and K.-A. Lie.  Fast computation of multiphase flow in porous media by implicit discontinuous Galerkin schemes with optimal ordering of elements.  _J. Comput. Phys._ , 227(24):10108–10124, 2008\.  doi: 10.1016/j.jcp.2008.08.024.
  * Natvig et al. [2007] J. R. Natvig, K.-A. Lie, B. Eikemo, and I. Berre.  An efficient discontinuous Galerkin method for advective transport in porous media.  _Adv. Water Resour._ , 30(12):2424–2438, 2007\.  doi: 10.1016/j.advwatres.2007.05.015.
  * OPM [2019] OPM.  The open porous media (OPM) initative.  https://opm-project.org/, 2019.
  * Rasmussen and Lie [2014] A. F. Rasmussen and K.-A. Lie.  Discretization of flow diagnostics on stratigraphic and unstructured grids.  In _16th European Conference on the Mathematics of Oil Recovery, Catalania, Sicily, Italy_. European Association of Geoscientits and Engineers, 2014.
  * Shahvali and Tchelepi [2013] M. Shahvali and H. A. Tchelepi.  Efficient coupling for non-linear multiphase flow with strong gravity.  In _SPE Reservoir Simulation Symposium, 18-20 February, The Woodlands, Texas, USA_. Society of Petroleum Engineers, February 2013.  doi: 10.2118/163659-MS.
  * Sheth and Younis [2017] S. M. Sheth and R. M. Younis.  Localized solvers for general full-resolution implicit reservoir simulation.  In _SPE Reservoir Simulation Conference_. Society of Petroleum Engineers, 2017.  doi: 10.2118/182691-MS.
  * The Open Porous Media Initative(2015) [OPM] The Open Porous Media (OPM) Initative.  The Norne dataset.  https://github.com/OPM/opm-data, 2015.
  * Trangenstein and Bell [1989] J. A. Trangenstein and J. B. Bell.  Mathematical structure of the black-oil model for petroleum reservoir simulation.  _SIAM J. Appl. Math._ , 49(3):749–783, 1989.  doi: 10.1137/0149044.
  * Trottenberg et al. [2000] U. Trottenberg, C. W. Oosterlee, and A. Schuller.  _Multigrid_.  Academic press, 2000.
  * Watts [1986] J. Watts.  A compositional formulation of the pressure and saturation equations.  _SPE Reservoir Eng._ , 1(03):243–252, 1986.  doi: 10.2118/12244-PA.


# On the Ricci curvature of K$\ddot{a}$hler-Ricci Flow

Fung Cheuk Yan

(Nov 24 2021)

###### Abstract

In this paper, we consider $n$-dimensional compact K$\ddot{a}$hler manifold
with semi-ample canonical line bundle under the long time solution of
K$\ddot{a}$hler Ricci Flow. In particular, if the Kodaira dimension is one,
Ricci curvature converge to negative of generalized K$\ddot{a}$hler Einstein
metric $\omega_{B}$ locally away from singular set in $C^{0}_{loc}(\omega(t))$
topology.

###### Contents

  1. 1 Introduction
    1. 1.1 Background
    2. 1.2 Major Result
  2. 2 Gradient estimate
    1. 2.1 Notation
    2. 2.2 Preliminary
    3. 2.3 Preparation for Theorem 2.1
    4. 2.4 Proof of Theorem 2.1
  3. 3 Ricci curvature estimate
    1. 3.1 Preliminary
    2. 3.2 Preparation for Proof of Theorem 1.2
    3. 3.3 Proof of Theorem 1.2

##  1 Introduction

###  1.1 Background

The setup is taken from [11] and [5].

Let $(X,\omega_{0})$ be a compact $n$-dimensional K$\ddot{a}$hler manifold
with semi-ample canonical line bundle $K_{X}$ and
$0<m:=\operatorname{Kod}(X)<n$.
Let $\omega(t)$ be the long time solution of the normalized K$\ddot{a}$hler-
Ricci flow

(1.1) $\dfrac{\partial}{\partial t}\omega=-\operatorname{Ric}(\omega)-\omega\
,\ \omega|_{t=0}=\omega_{0}$

There exists a holomorphic map:

$f:X\rightarrow
B=X_{can}\subset\mathbb{CP}^{n}:=\mathbb{C}H^{0}(X,K_{X}^{\otimes l})$

for sufficiently large $l\in{\mathbb{Z}^{+}}$ where $B=X_{can}$ is the
canonical model of $X$. We have $dimB=m$.
Let $S^{\prime}$ be the singular set of B together with the set of critical
values of $f$, and we define $S=f^{-1}(S^{\prime})\subset X$.
Note $f^{*}(\mathcal{O}(1))=K_{X}^{\otimes l}$. Let
$\chi=\frac{1}{l}\omega_{FS}$. $f^{*}\chi$ (also denoted by $\chi$) is a
smooth semi-positve representative of $-c_{1}(X)$ where ${\omega_{FS}}$ is the
Fubini-Study metric. Let $\Omega$ be the smooth volume form on $X$ with

$\sqrt{-1}\partial\bar{\partial}\log(\Omega)=\chi,\int_{X}\Omega=\binom{n}{m}\int_{X}\omega_{0}^{n-m}\wedge\chi^{m}$

Since $X_{y}:=f^{-1}(y)$ are Calabi-Yau for $y\in B\setminus S^{\prime}$,
there exists a unique smooth function $\rho_{y}$ on $X_{y}$ with
$\int_{X_{y}}\rho_{y}\omega_{0}^{n-m}=0$ and such that

$\omega_{0}|_{X_{y}}+\sqrt{-1}\partial\bar{\partial}\rho_{y}=:\omega_{y}$

is the unique Ricci-flat K$\ddot{a}$hler metric on ${X_{y}}$. Moreover,
$\rho_{y}$ depends smoothly on $y$, and so define a global smooth function on
$X\setminus S$. We define a a semi Ricci-flat form
$\omega_{SRF}=\omega_{0}+\sqrt{-1}\partial\bar{\partial}\rho$ in which its
restriction on fibre $X_{y}$ is Ricci-flat.
Define F on $X\setminus S$ by

(1.2) $F:=\dfrac{\Omega}{\omega^{n-m}_{SRF}\wedge\binom{n}{m}\chi^{m}}$

There exists a unqiue solution $v\in PSH(\chi)\cap C^{0}(B)\cap
C^{\infty}(B\setminus S^{\prime})$ of

(1.3) $(\chi+\sqrt{-1}\partial\bar{\partial}v)^{m}=Fe^{v}\chi^{m}$

Define a K$\ddot{a}$hler metric on $B$

$\omega_{B}:=\chi+\sqrt{-1}\partial\bar{\partial}v$

Note that $\omega_{B}$ is smooth and satisfies the generalized Einstein
equation

$\operatorname{Ric}(\omega_{B})=-\omega_{B}+\omega_{WP}\ on\ B\setminus
S^{\prime}$

where $\omega_{WP}$ is the Weil-Petersson metric induced by the Calabi-Yau
fibration $f$.

We also denote $f^{*}\omega_{B}$ by $\omega_{B}$.

In Theorem 1.1 and Theorem 1.2 of [11], it is shown that the metric
$\omega(t)$ converges to $\omega_{B}$ in the $C^{0}_{loc}(\omega_{0})$
topology on $X\setminus S$.

We consider the reference metric

(1.4) $\hat{\omega}=e^{-t}\omega_{0}+(1-e^{-t})\chi$

Some would consider another reference metric alternatively

(1.5) $\tilde{\omega}=e^{-t}\omega_{B}+(1-e^{-t})\omega_{SRF}$

Writing the metric as

$\omega(t)=\hat{\omega}(t)+\sqrt{-1}\partial\bar{\partial}\varphi$

The K$\ddot{a}$hler-Ricci Flow (1.1) is equivalent to the Monge-
Amp$\grave{e}$re equation

(1.6) $\dfrac{\partial}{\partial
t}\varphi=\log\dfrac{e^{(n-m)t}(\hat{\omega}(t)+\sqrt{-1}\partial\bar{\partial}\varphi(t))^{n}}{\Omega}-\varphi,\varphi(0)=0$

And we consider the function

(1.7) $u:=\frac{\partial\varphi}{\partial t}+\varphi-v$

The norm $|\cdot|$ in this article, unless specified, mean the norm
$\|\cdot\|_{\omega(t)}$.

In this article, we consider compact sets $K\subset\subset X\setminus S$.
Choose some open subset $U^{\prime}$ such that

(1.8) $f(K)\subset\subset U^{\prime}\subset\subset B\setminus S^{\prime}$

Define

(1.9) $U=f^{-1}(U^{\prime})$

Then

(1.10) $K\subset\subset U\subset\subset X\setminus S$

On any non-singular fibre $X_{y}$, define

$g_{F}(t):=g(t)|_{X_{y}}$

Locally, consider holomorphic coordinates as $(z_{1},z_{2},...,z_{n})$. Since
$f$ is a submersion in a neighborhood of point away from singularity, we can
write the coordinate of fiber as $\\{z_{\alpha}\\}_{\alpha=1}^{n-m}$. At a
fixed point $x$, if we can assume $\\{z_{\alpha}\\}_{\alpha=1}^{n-m}$
restricted on $X_{y}$ are normal coordinates and assume $g(t)$ to be identity.
And we can write:

$g_{F}(t)^{\alpha\bar{\beta}}=(g_{F}(t)_{\alpha\bar{\beta}})^{-1}$

###  1.2 Major Result

The uniform Ricci curvature bound on compact set away from singularities has
been studied for a long time. When $\operatorname{Kod}(X)=0$ and
$\operatorname{Kod}(X)=n$, the bound is shown by Cao [1] and Tsuji [12]
respectively. When $\operatorname{Kod}(X)=n-1$, since general fibre of $X$
over $X_{can}$ are complex tori, the estimate is given by the result of
Tosatti-Weinkove-Yang in [11]. The bound when $\operatorname{Kod}(X)=1$ is
shown by a recent paper of Jian-Song [7]. In general Kodaira dimension, when
regular fibers are biholomorphic to each other, Fong-Lee [2] showed the bound.
Jian-Shi [6] further extended the result and showed the local convergence of
Ricci curvature in such case:

$\norm{\operatorname{Ric}+\omega_{B}}_{C^{k}(K,\tilde{\omega}(t))}\leq
C(K)h_{K,k+2}(t)$

where $h_{K,k+2}(t)$ are positive functions which tends to zero as
$t\rightarrow\infty$, depending only on $k$ and the domain $K$ away from
singularities. And $\tilde{\omega}(t))$ is defined in (1.5). Jian-Shi
conjectured that will hold in the general case.

In this paper, inspired by the idea of Jian-Song [7] and the conjecture of
Jian-Shi [6] , we will show in Corollary 1.4 that, when
$\operatorname{Kod}(X)=1$, the Ricci curvature will converge to negative of
the generalized K$\ddot{a}$hler-Einstein metric $-\omega_{B}$ in
$C^{0}_{loc}(\omega(t))$. We will assume $\operatorname{Kod}(X)=1$ in the
whole paper.

There are two technical difficulties restricting us to consider the case
$\operatorname{Kod}(X)>1$. The first reason is that, when
$\operatorname{Kod}(X)=1$, it can be shown that Ricci curvature is negative
definite in the base $\times$ base direction. For detailed explanation, see
the paragraph at the end of Corollary 1.4. The second reason comes from a
computational obstacle of Laplacian of fibrewise average of the ricci
potential $u$. See the proof of Propositions 2.4 for detailed explanation.

The convergence of $\tr_{\omega}\omega_{B}-m$ is one of the most important
ingredients of the paper.

######  Propositions 1.1.

When $\operatorname{Kod}(X)=1$, for any compact set $K\subset\subset
X\setminus S$. There exists a positive decreasing function:

$h(t):\left[0,\infty\right)\rightarrow\left[0,\infty\right)$

such that the following holds:

$2|\tr_{\omega}{\omega_{B}}-1|\leq h(t)\ on\ K\times\left[0,\infty\right)$
$h(t)\rightarrow 0\ as\ t\rightarrow\infty\>$
$0\leq-h^{\prime}(t)\leq\frac{1}{2}h(t)$

###### Proof.

It is shown by Jian in Proposition 3.2 of [5]. ∎

The requirement $0\leq-h^{\prime}(t)\leq\frac{1}{2}h(t)$ is equivalent to
$(e^{\frac{1}{2}t}h(t))^{\prime}>0$. In particular, $h(t)^{-1}\leq e^{t/2}$
and $h(t)^{-1}\leq e^{t}h(t)$ .

The following theorem is an improvement of Proposition 7.1 of [7], where Jian-
Song showed that
$P:=g^{i\bar{j}}g^{\alpha\bar{\beta}}_{F}R_{i\bar{\beta}}R_{\alpha\bar{j}}$
where $i,j=1,...,n$ include both the base and fibre direction and
$\alpha,\beta=1,...,n-1$ represent the fibre direction, is bounded. We can
regard $2P$ as an upper bound of Ricci curvature in the fibre $\times$ fibre
direction and the fibre $\times$ base direction if we consider $g(p,t)=I$ at a
point.

The following theorem show $P$ converge to 0 locally away from singularities.

######  Theorem 1.2.

Let $(X,\omega_{0})$ be a compact $n$-dimensional K$\ddot{a}$hler manifold
with semi-ample canonical line bundle $K_{X}$ such that
$\operatorname{Kod}(X)=1$. Let $\omega(t)$ satisfy the K$\ddot{a}$hler-Ricci
Flow equation (1.1) on $X$ and $h(t)$ be defined as in Propositions 1.1. Then
for any compact set $K\subset\subset X\setminus S$, there exists a constant
$C=C(K)$ such that

$P:=g^{i\bar{j}}g^{\alpha\bar{\beta}}_{F}R_{i\bar{\beta}}R_{\alpha\bar{j}}\leq
C(K)h(t)^{\frac{1}{2}}$

The breakthrough comes from a sharper estimate of $|\nabla^{F}u|^{2}$ where it
is shown in Proposition 6.1 of [7] that $|\nabla^{F}u|^{2}\leq C(K)e^{-t}$
locally away from singularities. In Theorem 2.1, we will show a sharper
estimate. The main breakthrough comes from computational techniques which
consider as well the rate of convergence of $|\tr_{\omega}\omega_{B}-1|$.

For the reason why we have to consider $\operatorname{Kod}(X)=1$. See the
proof of Propositions 2.4 for detailed explanation.

The proof of Theorem 1.2 will be shown in Section 3.

The idea of the following corollaries comes from the proof of Theorem 7.1 of
[7].

######  Corollary 1.3.

When $\operatorname{Kod}(X)=1$, for any $K\subset\subset X\setminus S$, there
exists a positive decreasing function $h_{1}(t)$ such that on
$K\times\left[0,\infty\right)$, we have

$|\|\operatorname{Ric}\|_{\omega(t)}^{2}-1|\leq h_{1}(t)$ $h_{1}(t)\rightarrow
0\ as\ t\rightarrow\infty\>$

###### Proof.

By Theorem 1.1 of [5] and Theorem 1.2, there exists a positive decreasing
function $A(t)$ that tends to 0 as $t\rightarrow\infty$ such that on
$K\times[0,\infty)$, we have :

$P+|R+1|\leq A(t)$

Fix a point $(p,t)$. Consider a local coordinate where
$\\{z_{\alpha}\\}_{\alpha=1}^{n-1}$ is the coordinate of fibre and
$\\{z_{n}\\}$ is the coordinate of base. And at $(p,t)$, $g(p,t)$ is identity.

In such coordinate, we have :

$P=\sum_{i=1}^{n}\sum_{\alpha=1}^{n-1}|R_{i\bar{\alpha}}|^{2}$
$\displaystyle|\|\operatorname{Ric}\|_{\omega(t)}^{2}-1|$
$\displaystyle=||R_{n\bar{n}}|^{2}+\sum_{(i,j)\neq(n,n)}|R_{i\bar{j}}|^{2}-1|$
$\displaystyle\leq||R-\sum_{i\neq n}R_{i\bar{i}}|^{2}-1|+2P$
$\displaystyle\leq|R^{2}-1|+2|R\sum_{i\neq n}R_{i\bar{i}}|+|\sum_{i\neq
n}R_{i\bar{i}}|^{2}+2P$ $\displaystyle\leq C|(R+1)|+CP^{\frac{1}{2}}+CP$
$\displaystyle\leq CA(t)^{\frac{1}{2}}$

∎

######  Corollary 1.4.

When $\operatorname{Kod}(X)=1$, for any $K\subset\subset X\setminus S$, there
exists a positive decreasing function $h_{2}(t)$ such that on
$K\times\left[0,\infty\right)$, we have

$\|\operatorname{Ric}+\omega_{B}\|_{\omega(t)}^{2}\leq h_{2}(t)$
$h_{2}(t)\rightarrow 0\ as\ t\rightarrow\infty\>$

###### Proof.

By Theorem 1.1 of [5], Proposition 3.2 in [5] and Theorem 1.2, there exists a
positive decreasing function $A(t)$ that tends to 0 as $t\rightarrow\infty$
such that on $K\times[0,\infty)$, we have :

$P+|R+1|+|\|\omega_{B}\|_{\omega(t)}^{2}-1|\leq A(t)$

Choose coordinate as in Corollary 1.3 and let $n$ to be the base direction, we
have :

$\displaystyle\|\operatorname{Ric}+\omega_{B}\|_{\omega(t)}^{2}$
$\displaystyle=(R_{n\bar{n}}+\omega_{Bn\bar{n}})^{2}+\sum_{(i,j)\neq(n,n)}|{R_{i\bar{j}}}|^{2}$
$\displaystyle\leq(R_{n\bar{n}}+1+\omega_{Bn\bar{n}}-1)^{2}+2P$
$\displaystyle\leq 2(R_{n\bar{n}}+1)^{2}+2(\omega_{Bn\bar{n}}-1)^{2}+2P$
$\displaystyle\leq C|(R+1-\sum_{i\neq
n}R_{i\bar{i}})|^{2}+2(\omega_{Bn\bar{n}}-1)^{2}(\omega_{Bn\bar{n}}+1)^{2}+2P$
$\displaystyle\leq C|R+1|^{2}+2C|(R+1)\sum_{i\neq
n}R_{i\bar{i}}|+C|\sum_{i\neq n}R_{i\bar{i}}|^{2}$
$\displaystyle+2(\|\omega_{B}\|_{\omega(t)}^{2}-1)^{2}+2P$ $\displaystyle\leq
A(t)$

∎

Notice that even if Theorem 1.2 is generalized to other
$\operatorname{Kod}(X)$ i.e. Ricci curvature along Fibre $\times$ Fibre and
Fibre $\times$ Base direction tends to 0 in general $\operatorname{Kod}(X)$,
it is not sufficient to generalize Corollary 1.4. The reason is that, when
$\operatorname{Kod}(X)=1$, we make use of the convergence of scalar curvature
to -1 to conclude Ricci curvature is negative definite in the base $\times$
base direction. And this argument does not hold for $\operatorname{Kod}(X)>1$.

The outline of paper is as follows. In Section 2, we will show Theorem 2.1,
which is required to show Theorem 1.2. In Section 3, the proof of Theorem 1.2
will be shown.

Acknowledgement: The author would like to thank his advisor Frederick Fong
Tsz-Ho for his guide and support. The author would also like to thank Wangjian
Jian for his useful discussion.

##  2 Gradient estimate

In this section, we will improve the gradient estimate shown in Propositions
6.1 of [7], where Jian-Song showed that $|\nabla^{F}u|^{2}\leq C(K)e^{-t}$
locally.

######  Theorem 2.1.

Assume $\operatorname{Kod}(X)=1$, for any $K\subset\subset X\setminus S$, let
$h(t)$ be defined in Propositions 1.1. There exists a constant $C(K)$
depending compact set $K$ away from singularities such that on
$K\times\left[0,\infty\right)$, we have:

$|\nabla^{F}u|^{2}\leq C(K)e^{-t}h(t)$

Before proving the theorem, We have to introduce certain notation and compute
certain quantity first.

###  2.1 Notation

Recall in (1.7), $u$ is defined as

$u:=\frac{\partial\varphi}{\partial t}+\varphi-v$

Please note that there are slight difference in definition of ”$u$” in this
paper and [7]. The ”$u$” defined in [7], denoted by $f$ here is defined as
below:

(2.1) $f:=\frac{\partial\varphi}{\partial t}+\varphi$

Define $\bar{u}$ to be fibrewise average of u :

(2.2) $\bar{u}:=\dfrac{\int_{X_{y}}u\ \omega^{n-1}}{\int_{X_{y}}\omega^{n-1}}$

Define

$a(t):=\dfrac{1}{\int_{X_{y}}\omega^{n-1}}$

Define the fibrewise gradient :

$\nabla^{F}u:=(\nabla|_{X_{y}})(u|_{X_{y}})$

The norm of fibrewise gradient of u is given by the following :

$|\nabla^{F}u|^{2}=(g_{F})^{\alpha\bar{\beta}}u_{\alpha}u_{\bar{\beta}}$

We also use the following notation:

$|\nabla\nabla^{F}u|^{2}:=g^{i\bar{j}}g^{\alpha\bar{\beta}}_{F}(\nabla_{i}\nabla_{\alpha}u)(\nabla_{\bar{j}}\nabla_{\bar{\beta}}u)$
$|\nabla\bar{\nabla}^{F}u|^{2}:=g^{i\bar{j}}g^{\alpha\bar{\beta}}_{F}(\nabla_{\bar{j}}\nabla_{\alpha}u)(\nabla_{i}\nabla_{\bar{\beta}}u)$
$\langle\nabla^{F}\tr_{\omega}\omega_{B}\cdot\bar{\nabla}^{F}u\rangle=g_{F}^{\alpha\bar{\beta}}\nabla_{\alpha}\tr_{\omega}\omega_{B}\cdot\nabla_{\bar{\beta}}u$

###  2.2 Preliminary

######  Propositions 2.2.

For any $K\subset\subset X\setminus S$, on $K\times\left[0,\infty\right)$, we
have:

$\left(\dfrac{\partial}{\partial t}-\Delta\right)u=\tr_{\omega}\omega_{B}-1$
$\left(\dfrac{\partial}{\partial t}-\Delta\right)\tr_{\omega}\omega_{B}\leq
C(K)$

###### Proof.

They are shown in Lemma 3.2 and Lemma 3.3 in [11] by Tosatti-Weinkove-Yang.
The original proof comes from Schwarz Lemma computation from Theorem 4.3 of
[9] by Song-Tian. ∎

######  Propositions 2.3.

For any $K\subset\subset X\setminus S$, on $K\times\left[0,\infty\right)$, we
have:

$|u-\bar{u}|+|\nabla^{F}u|^{2}\leq C(K)e^{-t}$

###### Proof.

Let $\bar{f}$ denote fiberwise average of $f$ as in the definition of
$\bar{u}$ in (2.2). By Proposition 6.1 of [7] by Jian-Song, we have :

$|f-\bar{f}|+|\nabla^{F}f|^{2}\leq C(K)e^{-t}$

Note that $f-u=v$ and $v$ defined in (1.3) is a function depends on bases
only. Note $|\nabla^{F}v|=0$ and the fiberwise average $\bar{v}$ of $v$ is
equal to $v$. So we have

$|u-\bar{u}|+|\nabla^{F}u|^{2}\leq C(K)e^{-t}$

∎

######  Propositions 2.4.

For any $K\subset\subset X\setminus S$, on $K\times\left[0,\infty\right)$, we
have:

$\left(\dfrac{\partial}{\partial
t}-\Delta\right)\bar{u}=a(t)\int_{X_{y}}(\tr_{\omega}\omega_{B}-1)\omega^{n-1}+a(t)\int|\nabla^{F}u|^{2}\omega^{n-1}$

###### Proof.

It is an analogue of Lemma 6.2 in [7] by Jian-Song. Note that we require
$\operatorname{Kod}(X)=1$ here. We have to make use of Lemma 6.1 in [7] which
assumed $\operatorname{Kod}(X)=1$ to show

$\Delta\bar{u}=\dfrac{\int_{X_{y}}\Delta u\
\omega^{n-1}}{\int_{X_{y}}\omega^{n-1}}$

If this obstacle can be overcome, Theorem 1.2 can be generalized to other
$\operatorname{Kod}(X)$. ∎

######  Propositions 2.5.

For any $K\subset\subset X\setminus S$, on $K\times\left[0,\infty\right)$, we
have:

$\left(\dfrac{\partial}{\partial
t}-\Delta\right)|\nabla^{F}u|^{2}=|\nabla^{F}u|^{2}-|\nabla\nabla^{F}u|^{2}-|\nabla\bar{\nabla}^{F}u|^{2}+2\real\langle\nabla^{F}\tr_{\omega}\omega_{B}\cdot\bar{\nabla}^{F}u\rangle$

###### Proof.

It is an analogue of Lemma 6.4 from [7] by Jian-Song. ∎

The following allows us to consider cutoff function.

######  Lemma 2.6.

Let $K\subset\subset X\setminus S$ be a compact set. Choose some open subset
$U^{\prime}$ such that $f(K)\subset\subset U^{\prime}\subset\subset B\setminus
S^{\prime}$. Let $U=f^{-1}(U^{\prime})$. Then there exists a smooth cutoff
function $\rho$ such that $\rho$ is compactly supported on $U$, $0\leq\rho\leq
1$ $\rho\equiv 1$ on $K$, and for each K$\ddot{a}$hler-Ricci flow $\omega(t)$
there exists $1\leq C$ depending on the initial metric $\omega_{0}$ and $U$
with

$(|\partial\rho|^{2}_{\omega(t)}+|\Delta_{\omega(t)}\rho|)\leq C\ on\
U\times\left[0,\infty\right)$

###### Proof.

It is taken from Lemma 4.2 in [5] by Jian. ∎

###  2.3 Preparation for Theorem 2.1

######  Propositions 2.7.

For any $K\subset\subset X\setminus S$, let $h(t)$ be defined as in
Proposition 1.1, on $K\times\left[0,\infty\right)$, we have:

$|\left(\dfrac{\partial}{\partial t}-\Delta\right)(u-\bar{u})|\leq C(K)h(t)$

###### Proof.

Using formulas for evolution of $u$ and $\bar{u}$ in Propositions 2.2 and
Propositions 2.3, we have:

$\displaystyle|\left(\dfrac{\partial}{\partial t}-\Delta\right)(u-\bar{u})|$
$\displaystyle=|\tr_{\omega}\omega_{B}-1-a(t)\int_{X_{y}}(\tr_{\omega}\omega_{B}-1)\omega^{n-1}$
$\displaystyle-a(t)\int|\nabla^{F}u|^{2}\omega^{n-1}|$ $\displaystyle\leq
h(t)+h(t)a(t)\int_{X_{y}}\omega^{n-1}+Ce^{-t}a(t)\int_{X_{y}}\omega^{n-1}$
$\displaystyle\leq Ch(t)$

In the second last line, we used estimate of $|\nabla^{F}u|^{2}$ in
Propositions 2.4. ∎

######  Propositions 2.8.

For any $K\subset\subset X\setminus S$, let $h(t)$ be defined as in
Proposition 1.1, on $K\times\left[0,\infty\right)$, we have:

$\left(\dfrac{\partial}{\partial
t}-\Delta\right)(e^{t}h(t)^{-1}|\nabla^{F}u|^{2})\leq
Ce^{2t}|\nabla^{F}u|^{2}+h(t)^{-2}|\nabla^{F}\tr_{\omega}\omega_{B}|^{2}$

###### Proof.

Using formula of evolution of $|\nabla^{F}u|^{2}$ from Propositions 2.5, we
have:

$\displaystyle\left(\dfrac{\partial}{\partial
t}-\Delta\right)(e^{t}h(t)^{-1}|\nabla^{F}u|^{2})$
$\displaystyle=(e^{t}h(t)^{-1}+e^{t}\left(-\dfrac{h^{\prime}(t)}{h(t)}\right)h(t)^{-1})|\nabla^{F}u|^{2}$
$\displaystyle+e^{t}h(t)^{-1}\left(|\nabla^{F}u|^{2}-|\nabla\nabla^{F}u|^{2}-|\nabla\bar{\nabla}^{F}u|^{2}+2\real\langle\nabla^{F}\tr_{\omega}\omega_{B}\cdot\bar{\nabla}^{F}u\rangle\right)$
$\displaystyle\leq Ce^{t}h(t)^{-1}|\nabla^{F}u|^{2}$
$\displaystyle+e^{t}h(t)^{-1}\left(|\nabla^{F}u|^{2}-0-0\right)+2e^{t}h(t)^{-1}\real\langle\nabla^{F}\tr_{\omega}\omega_{B}\cdot\bar{\nabla}^{F}u\rangle$
$\displaystyle\leq
Ce^{t}h(t)^{-1}|\nabla^{F}u|^{2}+2h(t)^{-1}|\nabla^{F}\tr_{\omega}\omega_{B}|e^{t}|\nabla^{F}u|$
$\displaystyle\leq
Ce^{2t}|\nabla^{F}u|^{2}+h(t)^{-2}|\nabla^{F}\tr_{\omega}\omega_{B}|^{2}$

In the last step, we use Cauchy inequality and $e^{2t}|\nabla^{F}u|^{2}$
absorb the term $Ce^{t}h(t)^{-1}|\nabla^{F}u|^{2}$. ∎

######  Propositions 2.9.

For any $K\subset\subset X\setminus S$, let $h(t)$ be defined as in
Proposition 1.1. Choosing large enough $C(K)$ such that
$|u-\bar{u}|\leq\frac{1}{2}C(K)e^{-t}$ by Propositions 2.3 and applying
Propositions 2.7 for estimate of $u-\bar{u}$, on
$K\times\left[0,\infty\right)$, we have:

$\left(\dfrac{\partial}{\partial
t}-\Delta\right)\log(\dfrac{Ce^{-t}}{Ce^{-t}-(u-\bar{u})})\leq
e^{t}h(t)-Ce^{2t}|\nabla^{F}u|^{2}+C$

###### Proof.

$\displaystyle\left(\dfrac{\partial}{\partial
t}-\Delta\right)\log(\dfrac{Ce^{-t}}{Ce^{-t}-(u-\bar{u})})$
$\displaystyle=-\dfrac{\left(\frac{\partial}{\partial
t}-\Delta\right)(Ce^{-t}-(u-\bar{u}))}{Ce^{-t}-(u-\bar{u})}-\dfrac{|\nabla(u-\bar{u})|^{2}}{(Ce^{-t}-(u-\bar{u}))^{2}}-1$
$\displaystyle=-\dfrac{(-Ce^{-t}-\left(\frac{\partial}{\partial
t}-\Delta\right)(u-\bar{u}))}{Ce^{-t}-(u-\bar{u})}-\dfrac{|\nabla(u-\bar{u})|^{2}}{(Ce^{-t}-(u-\bar{u}))^{2}}-1$
$\displaystyle\leq\dfrac{e^{-t}+h(t)}{\frac{1}{2}e^{-t}}-\dfrac{|\nabla^{F}u|^{2}}{4Ce^{-2t}}+0$
$\displaystyle\leq e^{t}h(t)-Ce^{2t}|\nabla^{F}u|^{2}+C$

∎

Note that, by our choice of $C$, we have

$\log\frac{1}{2}\leq\log(\dfrac{Ce^{-t}}{Ce^{-t}-(u-\bar{u})})\leq\log 2$

######  Propositions 2.10.

For any $K\subset\subset X\setminus S$. let $h(t)$ be defined as in
Proposition 1.1. Using upper bound of evolution of $\tr_{\omega}\omega_{B}$ in
Propositions 2.2, on $K\times\left[0,\infty\right)$, we have:

$\left(\dfrac{\partial}{\partial
t}-\Delta\right)(h(t)^{-2}(\tr_{\omega}\omega_{B}-1+h(t))^{2})\leq
C+Ch(t)^{-1}-2h(t)^{-2}|\nabla^{F}\tr_{\omega}\omega_{B}|^{2}$

###### Proof.

$\left(\dfrac{\partial}{\partial
t}-\Delta\right)(h(t)^{-2}(\tr_{\omega}\omega_{B}-1+h(t))^{2})$

$\displaystyle=2h(t)^{-2}(\frac{-h^{\prime}(t)}{h(t)})(\tr_{\omega}\omega_{B}-1+h(t))^{2}$
$\displaystyle+h(t)^{-2}\left(2(\tr_{\omega}\omega_{B}-1+h(t))(\frac{\partial}{\partial
t}-\Delta)(\tr_{\omega}\omega_{B}+h(t))-2|\nabla\tr_{\omega}\omega_{B}|^{2}\right)$
$\displaystyle\leq Ch(t)^{-2}(Ch(t))^{2}$
$\displaystyle+h(t)^{-2}\left(2(Ch(t))(C)-2|\nabla\tr_{\omega}\omega_{B}|^{2}\right)$
$\displaystyle=C+Ch(t)^{-1}-2h(t)^{-2}|\nabla\tr_{\omega}\omega_{B}|^{2}$
$\displaystyle\leq
C+Ch(t)^{-1}-2h(t)^{-2}|\nabla^{F}\tr_{\omega}\omega_{B}|^{2}$

Note that we only have an upper bound of evolution of $\tr_{\omega}\omega_{B}$
and we do not have a lower bound for it. Note that in the first inequality,
$\tr_{\omega}\omega_{B}-1+h(t)$ is a non-negative function and we can apply
upper bound of evolution of $\tr_{\omega}\omega_{B}$. ∎

Note also, by Propositions 1.1, $h(t)^{-2}(\tr_{\omega}\omega_{B}-1+h(t))^{2}$
is bounded.

###  2.4 Proof of Theorem 2.1

Now we are ready to prove Theorem 2.1.

###### Proof.

Let $K\subset\subset X\setminus S$. Define

$Q:=e^{t}h(t)^{-1}|\nabla^{F}u|^{2}+C_{2}\log(\dfrac{Ce^{-t}}{Ce^{-t}-(u-\bar{u})})+C_{1}h(t)^{-2}(\tr_{\omega}\omega_{B}-1+h(t))^{2}$

By Proposition 2.8, Proposition 2.9 and Proposition 2.10, and choosing large
$C_{1}$ and $C_{2}$,

$\displaystyle\left(\dfrac{\partial}{\partial t}-\Delta\right)Q$
$\displaystyle\leq-e^{2t}|\nabla^{F}u|^{2}+C+h(t)^{-1}+e^{t}h(t)$
$\displaystyle\leq-e^{2t}|\nabla^{F}u|^{2}+e^{t}h(t)$

In the last step, $e^{t}h(t)$ absorb $h(t)^{-1}$ as explained at the end of
Propositions 1.1.

By Lemma 2.6, let $\rho$ be cutoff function s.t. $|\nabla\rho|\leq
C,|\Delta\rho|\leq C$ and $\rho\equiv 1$ on $K$, $supp\ \rho\subset\subset U.$

$\displaystyle(\dfrac{\partial}{\partial t}-\Delta)\left(\rho^{4}Q\right)$
$\displaystyle=\rho^{4}(\dfrac{\partial}{\partial
t}-\Delta)Q-Q\Delta\rho^{4}-2\real\langle\nabla
Q\cdot\bar{\nabla}\rho^{4}\rangle$
$\displaystyle\leq\rho^{4}(-e^{2t}|\nabla^{F}u|^{2}+e^{t}h(t))$
$\displaystyle-Q\Delta\rho^{4}-2\real\langle\nabla
Q\cdot\bar{\nabla}\rho^{4}\rangle$

$-Q\Delta\rho^{4}\leq C\rho^{2}Q\leq Ch(t)^{-1}\leq Ce^{t}h(t)$

$\displaystyle-2\real(\nabla Q\cdot\bar{\nabla}\rho^{4})$
$\displaystyle=-8\rho^{-1}\real\langle\nabla(\rho^{4}Q)\cdot\bar{\nabla}\rho\rangle+32\rho^{2}Q|\nabla\rho|^{2}$
$\displaystyle\leq-8\rho^{-1}\real\langle\nabla(\rho^{4}Q)\cdot\bar{\nabla}\rho\rangle+Ce^{t}h(t)$

Combining above, we have:

$(\dfrac{\partial}{\partial
t}-\Delta)\left(\rho^{4}Q\right)\leq-e^{2t}\rho^{4}|\nabla^{F}u|^{2}+e^{t}h(t)-8\rho^{-1}\real\langle\nabla(\rho^{4}Q)\cdot\bar{\nabla}\rho\rangle$

At max point, we have :

$\rho^{4}|\nabla^{F}u|^{2}\leq Ce^{-t}h(t)$

As explained at the end of Proposition 2.9 and Proposition 2.10, the two terms
below is bounded :

$C_{2}\log(\dfrac{Ce^{-t}}{Ce^{-t}-(u-\bar{u})})+C_{1}h(t)^{-2}(\tr_{\omega}\omega_{B}-1+h(t))^{2}$

we have finished the proof. ∎

##  3 Ricci curvature estimate

###  3.1 Preliminary

######  Propositions 3.1.

For any $K\subset\subset X\setminus S$, let $P$ be defined as in Theorem 1.2.
on $K\times\left[0,\infty\right)$, we have:

$P\leq C(K)$

###### Proof.

It is taken from Proposition 7.1 of [7] by Jian-Song. ∎

######  Propositions 3.2.

For any $K\subset\subset X\setminus S$, let $P$ be defined as in Theorem 1.2.
on $K\times\left[0,\infty\right)$, we have:

$\left(\dfrac{\partial}{\partial t}-\Delta\right)P\leq C(K)e^{t}\\\ $

###### Proof.

By Lemma 7.2 of [7] by Jian-Song,

$\left(\dfrac{\partial}{\partial t}-\Delta\right)P\leq C(K)e^{t}+C(K)e^{t}P$

By Propositions 3.1, we have show the Propositions. ∎

###  3.2 Preparation for Proof of Theorem 1.2

######  Propositions 3.3.

For any $K\subset\subset X\setminus S$, let $h(t)$ be defined as in
Proposition 1.1 and $P$ be defined as in Theorem 1.2. on
$K\times\left[0,\infty\right)$, we have:

$\left(\dfrac{\partial}{\partial t}-\Delta\right)(h(t)^{-1}P^{2})\leq
Ch(t)^{-1}e^{t}P+Ch(t)^{-1}$

###### Proof.

$\displaystyle\left(\dfrac{\partial}{\partial
t}-\Delta\right)(h(t)^{-1}P^{2})$
$\displaystyle=\left(-\dfrac{h^{\prime}(t)}{h(t)}\right)h(t)^{-1}P^{2}+h(t)^{-1}\left(\dfrac{\partial}{\partial
t}-\Delta\right)P^{2}$ $\displaystyle\leq
Ch(t)^{-1}+h(t)^{-1}(2P\left(\dfrac{\partial}{\partial
t}-\Delta\right)P-2|\nabla P|^{2})$ $\displaystyle\leq
Ch(t)^{-1}+Ch(t)^{-1}e^{t}P$

In the first inequality,we have used property of $h(t)$ from its definition in
Propositions 1.1. In the last step we have used Propositions 3.2 for the
estimate of evolution of $P$ and Propositions 3.1 for an upper bound for $P$.
∎

The idea of considering $P^{2}$ comes from Proof of Theorem 1.6 in [6] by
Jian-Shi.

######  Propositions 3.4.

For any $K\subset\subset X\setminus S$,let $h(t)$ be defined as in Proposition
1.1 and $P$ be defined as in Theorem 1.2. on $K\times\left[0,\infty\right)$,
we have:

$\left(\dfrac{\partial}{\partial
t}-\Delta\right)(h(t)^{-1}e^{t}|\nabla^{F}u|^{2})\leq-
Ch(t)^{-1}e^{t}P+Ce^{t}h(t)+h(t)^{-2}|\nabla^{F}\tr_{\omega}\omega_{B}|^{2}$

###### Proof.

Using formula of evolution of $|\nabla^{F}u|^{2}$ from Propositions 2.5, we
have:

$\displaystyle\left(\dfrac{\partial}{\partial
t}-\Delta\right)(h(t)^{-1}e^{t}|\nabla^{F}u|^{2})$
$\displaystyle=(h(t)^{-1}e^{t}+e^{t}\left(-\dfrac{h^{\prime}(t)}{h(t)}\right)h(t)^{-1})|\nabla^{F}u|^{2}$
$\displaystyle+h(t)^{-1}e^{t}\left(|\nabla^{F}u|^{2}-|\nabla\nabla^{F}u|^{2}-|\nabla\bar{\nabla}^{F}u|^{2}+2\real\langle\nabla^{F}\tr_{\omega}\omega_{B}\cdot\bar{\nabla}^{F}u\rangle\right)$
$\displaystyle\leq Ch(t)^{-1}e^{t}|\nabla^{F}u|^{2}$
$\displaystyle+h(t)^{-1}e^{t}\left(|\nabla^{F}u|^{2}-0-|\nabla\bar{\nabla}^{F}u|^{2}\right)+2(h(t)^{-1}|\nabla^{F}\tr_{\omega}\omega_{B}|)(e^{t}|\nabla^{F}u|)$
$\displaystyle\leq-
Ch(t)^{-1}e^{t}|\nabla\bar{\nabla}^{F}u|^{2}+Ce^{2t}|\nabla^{F}u|^{2}+h(t)^{-2}|\nabla^{F}\tr_{\omega}\omega_{B}|^{2}$
$\displaystyle\leq-
Ch(t)^{-1}e^{t}|\nabla\bar{\nabla}^{F}u|^{2}+Ce^{t}h(t)+h(t)^{-2}|\nabla^{F}\tr_{\omega}\omega_{B}|^{2}$
$\displaystyle=-Ch(t)^{-1}e^{t}P+Ce^{t}h(t)+h(t)^{-2}|\nabla^{F}\tr_{\omega}\omega_{B}|^{2}$

In the second last line, we applied Theorem 2.1 for the estimate of
$|\nabla^{F}u|^{2}$.

In the last step, by direct computation, we have

$\partial\bar{\partial}u=-\operatorname{Ric}-\omega_{B}$

Also, note $\omega_{B}$ vanish on fibre $\times$ fibre and fibre $\times$ base
direction. We have:

$|\nabla\bar{\nabla}^{F}u|^{2}=P$

∎

Note by Theorem 2.1, $h(t)^{-1}e^{t}|\nabla^{F}u|^{2}$ is bounded.

###  3.3 Proof of Theorem 1.2

###### Proof.

Let $K\subset\subset X\setminus S$. Define:

$Q:=h(t)^{-1}P^{2}+C_{1}h(t)^{-1}e^{t}|\nabla^{F}u|^{2}+C_{2}h(t)^{-2}(\tr_{\omega}\omega_{B}-1+h(t))^{2}$

By Proposition 2.10, Proposition 3.3, Proposition 3.4 and choosing large
enough $C_{1}$ and $C_{2}$, we have:

$\displaystyle\left(\dfrac{\partial}{\partial t}-\Delta\right)Q$
$\displaystyle\leq-h(t)^{-1}e^{t}P+Ce^{t}h(t)+Ch(t)^{-1}$
$\displaystyle\leq-h(t)^{-1}e^{t}P+Ce^{t}h(t)$

By Lemma 2.6, let $\rho$ be cutoff function s.t. $|\nabla\rho|\leq
C,|\Delta\rho|\leq C$ and $\rho\equiv 1$ on $K$, $supp\ \rho\subset\subset U.$

$\displaystyle(\dfrac{\partial}{\partial t}-\Delta)\left(\rho^{4}Q\right)$
$\displaystyle=\rho^{4}(\dfrac{\partial}{\partial
t}-\Delta)Q-Q\Delta\rho^{4}-2\real\langle\nabla
Q\cdot\bar{\nabla}\rho^{4}\rangle$
$\displaystyle\leq\rho^{4}(-h(t)^{-1}e^{t}P+e^{t}h(t))$
$\displaystyle-Q\Delta\rho^{4}-2\real\langle\nabla
Q\cdot\bar{\nabla}\rho^{4}\rangle$

$-Q\Delta\rho^{4}\leq C\rho^{2}Q\leq Ch(t)^{-1}\leq Ce^{t}h(t)$

$\displaystyle-2\real\langle\nabla Q\cdot\bar{\nabla}\rho^{4}\rangle$
$\displaystyle=-8\rho^{-1}\real\langle\nabla(\rho^{4}Q)\cdot\bar{\nabla}\rho\rangle+32\rho^{2}Q|\nabla\rho|^{2}$
$\displaystyle\leq-8\rho^{-1}\real\langle\nabla(\rho^{4}Q)\cdot\bar{\nabla}\rho\rangle+Ce^{t}h(t)$

Combining above, we have:

$(\dfrac{\partial}{\partial
t}-\Delta)\left(\rho^{4}Q\right)\leq-\rho^{4}h(t)^{-1}e^{t}P+e^{t}h(t)-8\rho^{-1}\real\langle\nabla(\rho^{4}Q)\cdot\bar{\nabla}\rho\rangle$

At max point, we have

$\rho^{4}P\leq h(t)^{2}$

In particular,

$\rho^{4}P^{2}\leq C\rho^{4}P\leq Ch(t)^{2}\leq h(t)$

Since other terms in Q are bounded (explained in Proposition 2.10 and
Proposition 3.4), we have proved the theorem. ∎

## References

  * [1] Cao, H.-D. Deformation of K$\ddot{a}$hler metrics to K$\ddot{a}$hler-Einstein metrics on compact K$\ddot{a}$hler manifolds, Invent. Math. 81 (1985), no. 2, 359–372
  * [2] Fong, F.T.-H., Lee, M.C. Higher-order estimates of long-time solutions to the K$\ddot{a}$hler-Ricci flow, J. Funct. Anal. 281 (2021), no. 11, Paper No. 109235
  * [3] Fong, F.T.-H., Zhang, Y. Local curvature estimates of long-time solutions to the K$\ddot{a}$hler-Ricci flow, Adv. Math. 375. (2020) 107416. MR 4170232
  * [4] Gill, M. Collapsing of products along the K$\ddot{a}$hler-Ricci flow, Trans. Amer. Math. Soc. 366 (2014), no. 7, 3907–3924.
  * [5] Jian, W. Convergence of scalar curvature of K$\ddot{a}$hler-Ricci flow on manifolds of positive Kodaira dimension, Adv. Math. 371. (2020) 107253. MR 4108223
  * [6] Jian, W., Shi, Y. A ”boundedness implies convergence” principle and its applications to collapsing estimates in K$\ddot{a}$hler geometry arXiv.1909.05521, accepted by Nonlinear Analysis
  * [7] Jian, W., Song J. Diameter and ricci curvature estimates for long-time solutions of the K$\ddot{a}$hler-ricci flow. arxiv.2101.04277
  * [8] Song, J., Tian, G. Canonical measures and K$\ddot{a}$hler-Ricci flow, J. Amer. Math. Soc. 25 (2012), no. 2, 303–353.
  * [9] Song, J., Tian, G. The K$\ddot{a}$hler-Ricci flow on surfaces of positive Kodaira dimension, Invent. Math. 170 (2007), no. 3, 609–653.
  * [10] Tosatti, V. Adiabatic limits of Ricci-flat K$\ddot{a}$hler metrics J. Differential Geom. 84 (2010), no.2, 427–453.
  * [11] Tosatti, V., Weinkove, B. and Yang, X. The K$\ddot{a}$hler-Ricci flow, Ricci-flat metrics and collapsing limit, Amer. J. Math. 140 (2018), no. 3, 653–698
  * [12] Tsuji, H. Existence and degeneration of K$\ddot{a}$hler-Einstein metrics on minimal algebraic varieties of general type, Math. Ann. 281 (1988), 123–133