∎

11institutetext: 1School of Earth Sciences and Engineering, Hohai University,
Nanjing, China.
2National Institute of Oceanography and Applied Geophysics OGS, Trieste,
Italy.
3KACST, PO Box 6086, Riyadh 11442, Saudi Arabia. 11email<EMAIL_ADDRESS>
($\ast$) corresponding author. 11email<EMAIL_ADDRESS>
# Rock acoustics of diagenesis and cementation

Short title: Seismic physics of rock burial history

José M. Carcione1,2 Davide Gei2 Stefano Picotti2 Ayman N. Qadrouh3 Mamdoh
Alajmi3 Jing Ba1(∗)

###### Abstract

We simulate the effects of diagenesis, cementation and compaction on the
elastic properties of shales and sandstones with four different petro-
elastical theories and a basin-evolution model, based on constant heating and
sedimentation rates. We consider shales composed of clay minerals, mainly
smectite and illite, depending on the burial depth, and the pore space is
assumed to be saturated with water at hydrostatic conditions. Diagenesis in
shale (smectite/illite transformation here) as a function of depth is
described by a 5th-order kinetic equation, based on an Arrhenius reaction
rate. On the other hand, quartz cementation in sandstones is based on a model
that estimates the volume of precipitated quartz cement and the resulting
porosity loss from the temperature history, using an equation relating the
precipitation rate to temperature. Effective pressure effects (additional
compaction) are accounted for by using Athy equation and the Hertz-Mindlin
model. The petro-elastic models yield similar seismic velocities, despite the
different level of complexity and physics approaches, with increasing density
and seismic velocities as a function of depth. The methodology provides a
simple procedure to obtain the velocity of shales and sandstones versus
temperature and pressure due to the diagenesis-cementation-compaction process.

###### Keywords:

shales, sandstone, diagenesis, cementation, compaction, seismic velocities,
granular media, Gassmann equation

††journal: Journal: Pure and Applied Geophysics

##  1 Introduction

Diagenesis is the chemical and physical process by which sediments which are
buried in the Earth crust are compacted by a process where minerals
precipitate from solution due to cementation, grain rearrangement, with a
reduction in porosity (Pytte and Reynolds, 1989; Walderhaug, 1996; Lander and
Walderhaug, 1999). Thus, sediments (clay and sand) become rocks by
lithification, i.e., the precipitated mineral create bonds between grains, and
sand becomes sandstone and clay-rich sediments become shale by this combined
diagenesis-cementation process. We do not consider the phenomenon by which
metamorphic rocks are formed, which typically occurs after that process,
beyond 200 ∘C, and before melting, roughly 800 ∘C, depending on the minerals
(Carcione et al., 2018). The diagenetic process requires the flow of water in
geological time for the minerals to precipitate and crystallize. Typical
cements include quartz, calcite and clay minerals.

The acoustics of diagenesis has scarcely been attacked. We can mention Draege
et al. (2006), who presented a methodology for the estimation of the shale
stiffnesses in the transition zone from mechanical compaction to chemical
compaction/cementing. The model showed consistent results when comparing
vertical P- and S-wave velocities with log data. Another type of diagenesis
process is the creation of kerogen and bitumen from which hydrocarbons (oil
and gas) are formed by the thermal alteration of these materials. A kinetic
model can describe this diagenetic process, to model the dissolution-
precipitation mechanism (e.g., Luo and Vasseur, 1996;). We do not consider
this mechanism in the present work, but it can easily be incorporated (Pinna
et al., 2011; Carcione et al., 2011; Carcione and Avseth, 2015).

We consider two dissimilar theories of the diagenesis-cementation process,
namely smectite/illite conversion based on a kinetic reaction given by the
Arrhenius equation (Pytte and Reynolds, 1989), and sandstone compaction and
cementation (Walderhaug, 1996; Avseth and Lehocki, 2016). The conversion from
smectite to illite (clay diagenesis) with increasing depth occurs in all
shales (Scotchman, 1987) and can be described by the widely accepted model
proposed by Pytte and Reynolds (1989) based on a 5th-order kinetic reaction of
the Arrhenius type. The result of the conversion is that the stiffnesses of
the minerals composing the shale increase with depth. To study quartz
cementation in sandstones, one has to estimate cement volumes mainly as a
function of temperature. According to Walderhaug (1996), the first stage is
precipitation, which occurs at temperatures above 60-80 ∘C, when the
cementation starts to be effective. The process depends on the geothermal
gradient and on the effective radii of the grains (and surface area), with
small quartz grains producing more cement (Walderhaug, 1996, Fig. 4).
Moreover, clay minerals or calcite cement covering the grains may reduce the
surface area, and inhibit dissolution. Roughly, the porosity loss is equal to
the volume of precipitated quartz, if the entire cement source is sourced from
outside the sediment volume and a volume of water equal to the cement volume
is removed from the sediment. An additional loss is accounted for compaction
through the intergranular volume index (IGV) that depends on effective stress
(Lander and Walderhaug, 1999). The model does not consider grain
interpenetration due to intergranular pressure solution, which is assumed to
be negligible, as supported by experimental evidence (Sippel, 1968; Paxton et
al., 2002). However, if required, such an effect can be implemented by using
the model of Weyl (1959) or that of Stephenson et al. (1992). Figures 1 and 2
show the microstructure of a shale and a sandstone, where the composition and
cementation can be appreciated. Other illustrative SEM (scanning electron
micrographs) images of several shales showing textures of montmorillonite
(smectite) and illitic shale can be seen in Keller et al. (1986, Figs. 4-22).
SEM images of several sandstones are shown in Walderhaug (2000, Fig 3), with
high and minor quartz-cement, pores surrounded by almost totally quartz-
cemented areas, and stylolite containing mica and detrital clay.

We assume a simple basin-evolution model with a constant sedimentation rate,
heating rate and geothermal gradient. The diagenesis process starts at a given
depth, pressure and temperature, with the pores filled with liquid water. Four
petro-elastical models yield the seismic velocity of the rocks. The first
model (Model 1) is based on an effective mineral moduli computed with the
Hashin-Shtrikmann average, Krief equation to obtain the dry-rock bulk and
shear moduli and Gassmann equations to estimate the wet-rock moduli. Model 2
is based on the self-consistent approximation to obtain the properties of
sand/clay mixtures (Gurevich and Carcione, 2000) and partially molten rocks
(Carcione et al., 2020). Model 3 is a generalization of Gassmann equation to
multi-mineral media (Carcione et al., 2005), a theory also used to describe
the properties gas-hydrate bearing sediments (Gei and Carcione, 2003). These
three models implement Athy equation (Athy, 1930) to account for the pressure
effects. The fourth approach (Model 4) is a modified version of the patchy
cement model introduced by Avseth et al. (2016), based on the contact-cement
theory (CCT) of Dvorkin and Nur (1996), assuming cement deposited at grain
contacts (scheme 1 in Mavko et al., 2009, p. 257) (see also Avseth et al.,
2010), and the Hertz-Mindlin theory and IGV concept to include pressure
effects. Then, the Voigt-Reuss Hill averages yield the wet-rock moduli, with
the self-consistent (SC) model used to obtain the properties of the mineral
mixture.

##  2 Temperature-pressure conditions

Let us assume a rock unit at depth $z$. The lithostatic (or confining)
pressure for an average sediment density $\bar{\rho}$ is $p_{c}=\bar{\rho}gz$,
where $g$ is the acceleration of gravity. On the other hand, the hydrostatic
pore pressure is approximately $p_{H}=\bar{\rho}_{w}gz$, where
$\bar{\rho}_{w}$ is the density of water (values of $\bar{\rho}$ = 2.5 g/cm3
and $\bar{\rho}_{w}$ = 1.04 g/cm3 are assumed here, and $g$ = 9.81 m/s2).

Effective pressure is defined here as

$p_{e}=p_{c}-p_{H},$ (1)

i.e., as differential pressure.

For a constant sediment burial rate, $S$, and a constant geothermal gradient,
$G$, the temperature variation of a particular sediment volume is

$T=T_{0}+Gz=T_{0}+Ht,\ \ \ z=St,\ \ \ H=GS$ (2)

with a surface temperature $T_{0}$ at time $t$ = 0, where $t$ is deposition
time and $H$ is the heating rate. Typical values of $G$ range from 20 to 40
oC/km, while $S$ may range between 0.02 and 0.5 km/m.y. (m.y. = million
years).

##  3 Diagenetic-compaction processes

###  3.1 Smectite/illite conversion in shales

To evaluate the amount of smectite/illite ratio forming the shale matrix is
important, since this ratio affects the density, stiffness moduli and wave
velocities of the rock. Shale mineralogy may include kaolinite,
montmorillonite-smectite, illite and chlorite, so the term smectite/illite as
used in this study may be representative for a mixture of clay minerals
(Mondol et al., 2008). The smectite/illite composite is subject to internal
hydration, so its mechanical properties such as the stiffnesses can vary
depending on the rock.

The conversion smectite/illite occurs in all shales with a general release of
bound water into the pore space (Scotchman, 1987). Smectite dehydration
implies a stiffer matrix due to the presence of more illite and therefore
higher velocities. The conversion depends on temperature and sedimentation
rate. A solution to this problem has been provided by Pytte and Reynolds
(1989). The process is pictorially explained in Figure 3 and mathematical
given in Appendix A.

Future works should consider another effect that could be important, i.e.,
quartz generation as a by-product of the diagenesis smectite/illite
conversion. Thyberg et al. (2009) show that this is another factor to explain
the velocity increase, due to micro-crystalline quartz-cementation of the rock
frame.

###  3.2 Cementation and compaction in sandstones.

The physico-chemical process of diagenesis-cementation-compaction in
sandstones is shown in Figure 4 and described mathematically in Appendix B by
using the theory developed by Walderhaug (1996) for a constant heating rate.
The model assumes that the quartz grains are spherical and have the same
radius. Surface area decreases as a function of porosity to account for
compaction and cementation. Compaction reduces quartz surface area by
increasing grain contact area and by injecting matrix material into pore
spaces (this last effect neglected here). Cementation can cause further
reduction of surface area when quartz grains are encased by pore-filling
cements.

##  4 Petro-elastical models

We consider four petro-elastical models to obtain the P- and S-wave velocities
of the rocks as a function of depth, pressure and temperature. In the
examples, Model 1, 2 and 3 are applied to shales and all four models to
sandstone. Some of the equations, such as Voigt, Reuss and Hashin-Shtrikmann
averages, are well-known in the rock-physics community. We refer to Mavko et
al. (2009) for details.

###  4.1 Model 1. Hashin-Shtrikmann-Gassmann model.

In the first model, smectite and illite (or cement and quartz) are “mixed” by
using the Hashin-Shtrikmann averages to obtain the bulk and shear moduli of
the mineral composing the frame, $K_{s}$ and $\mu_{s}$ (Appendix C),
respectively. Then, Krief equation (Krief et al., 1990) yields the dry-rock
moduli of the frame with porosity $\phi$,

$K_{m}=K_{s}(1-\phi)^{A/(1-\phi)},\ \ \ \mu_{m}=\mu_{s}(1-\phi)^{A/(1-\phi)},$
(3)

where $A$ is a dimensionless parameter (a value $A$ = 3 is assumed here).

The effect of the pore fluid can be accounted for by using Gassmann equations
(e.g., Mavko et al., 2009; Carcione, 2014). The Gassmann bulk and shear moduli
are given by

$K_{G}=K_{m}+\alpha^{2}M,\ \ \ \mu=\mu_{m},$ (4)

where

$M=\left(\displaystyle\frac{\alpha-\phi}{K_{s}}+\displaystyle\frac{\phi}{K_{f}}\right)^{-1},$
(5) $\alpha=1-\displaystyle\frac{K_{m}}{K_{s}}.$ (6)

and $K_{f}$ is the fluid modulus.

In shales, the reduction in volume due to release of bound water (see Figure
3) is a compaction effect that is modeled with Athy equation,

$\phi=\phi_{0}\exp(-\beta p_{e})$ (7)

where $\phi_{0}$ is the initial porosity and $\beta$ is an empirical constant
(Athy, 1930; Rieke III and Chilingarian, 1974 ) (see also Appendix B) (we
consider $\beta$ = 0.01/MPa in this work).

###  4.2 Model 2. Self-consistent (SC) scheme.

In the SC approximation, the elastic moduli of an unknown effective medium
have to be found implicitly. The model has been used by Gurevich and Carcione
(2000) to obtain the stiffnesses of sand-clay mixtures, where the inclusions
are spherical. Here, we consider spherical grains (aspect ratio $\gamma$ = 1)
and pores of aspect ratio $\gamma<1$. In this case, it is $n$ = 3, with
smectite, illite and water in the case of shales, and quartz, cement and water
in the case of sandstones. The proportion of phase 3 ($\phi_{3}$, the pore
space) is given by Athy equation (7).

The effective bulk and shear moduli of the composite medium ($K$ and $\mu$),
with $n$ phases and proportion $\phi_{i}$, are obtained as the roots of the
following system of equations

$\begin{array}[]{l}\sum_{i=1}^{n}\phi_{i}(K_{i}-K)P_{i}=0,\\\ \\\
\sum_{i=1}^{n}\phi_{i}(\mu_{i}-\mu)Q_{i}=0,\\\ \\\ \end{array}$ (8)

where

$\begin{array}[]{l}P_{i}=\displaystyle\frac{K+\frac{4}{3}\mu}{K_{i}+\frac{4}{3}\mu},\
\ \ i=1,\ldots,n\\\ \\\ Q_{i}=\displaystyle\frac{\mu+\zeta}{\mu_{i}+\zeta},\\\
\\\
\zeta=\displaystyle\frac{\mu}{6}\cdot\displaystyle\frac{9K+8\mu}{K+2\mu}\end{array}$
(9)

for the grains, i.e., $i$ =1, 2 (Mavko et al., 2009, p. 187; Carcione et al.,
2020) and $P_{3}=\frac{1}{3}T_{iijj}$ and $Q_{3}=\frac{1}{5}(T_{ijij}-P)$ (for
the pores), where $T_{iijj}$ and $T_{ijij}$ are given in Appendix A of
Berryman (1980) or in page 189 of Mavko et al. (2009). If $\gamma$ = 1,
$P_{3}$ and $Q_{3}$ are given by equation (9). A limitation of this theory is
that the inclusions are isolated, so that pore pressures are not equilibrated
and the model computes high-frequency velocities.

To solve equation (8), we use the algorithm developed by Goffe et al. (1994).
The Fortran code can be found in: https://econwpa.ub.uni-muenchen.de/econ-
wp/prog/papers/9406/9406001.txt.

###  4.3 Model 3. Generalized Gassmann (GG) theory

We also consider the composite model of Carcione et al. (2005), where the
porous medium is composed of $n$ = 2 two solids (smectite and illite in
shales; quartz and cement in sandstones) and a fluid. If $\phi_{i}$ is the
fraction of the $i$-th solid and $\phi$ is the porosity, it is
$\sum^{n}\phi_{i}+\phi$ = 1. The Gassmann modulus is

$K_{G}=\sum_{i=1}^{n}K_{mi}+\left(\sum_{i=1}^{n}\alpha_{i}\right)^{n}M,$ (10)

where

$M=\left(\sum_{i=1}^{n}\frac{\phi_{i}^{\prime}}{K_{i}}+\frac{\phi}{K_{f}}\right)^{-1},$
(11) $\phi_{i}^{\prime}=\alpha_{i}-\beta_{i}\phi,\ \ \
\alpha_{i}=\beta_{i}-\frac{K_{mi}}{K_{i}},\ \ \
\beta_{i}=\frac{\phi_{i}}{1-\phi}.$ (12)

$\beta_{i}$ is the fraction of solid $i$ per unit volume of total solid. Here
$K_{i}$, $i=1,\ldots n$ and $K_{f}$ are the solid and fluid bulk moduli,
respectively, and $K_{mi}$, $i=1,\ldots n$ are the frame moduli.

A generalization of Krief’s model for a multi-mineral porous medium is used to
obtain the frame moduli,

$K_{mi}=(K_{s}/V_{K})\beta_{i}K_{i}(1-\phi)^{A/(1-\phi)},\ \ \ \ i=1,\ldots,n$
(13)

(Carcione et al, 2005), $V_{K}=\sum_{i=1}^{n}\beta_{i}K_{i}$ is the bulk Voigt
average, and $K_{s}$ is given by equation (C.3). The expression (13) is such
that the composite modulus $K_{m}=\sum_{i}^{n}K_{mi}$ gives the Hashin-
Shtrikmann (HS) average when $\phi$ = 0.

Similarly,

$\mu_{mi}=(\mu_{s}/V_{\mu})\beta_{i}\mu_{i}(1-\phi)^{A/(1-\phi)},\ \ \ \
i=1,\ldots,n,$ (14)

where $V_{\mu}=\sum_{i=1}^{n}\beta_{i}\mu_{i}$ is the shear Voigt average, and
$\mu_{s}$ is given by equation (C.3). The dry-rock (and wet-rock) shear
modulus of the composite is $\mu_{m}=\sum_{i}^{n}\mu_{mi}$.

The pressure effects are described by Athy equation (7).

###  4.4 Model 4. Patchy cement (PC) model.

Dvorkin and Nur (1996) developed an elastic model [contact cement theory
(CCT)] based on spherical grains (see also Dvorkin et al., 1999). Let us
denote with subscripts 1 and 2 the properties of the grains and cement,
respectively. The model assumes that initially the sandstone is a random pack
of identical spherical grains with porosity near the critical one
($\phi_{0}\approx$ 0.36, a critical porosity) and an average number of grain
contacts equal to $C$ = 9, with bulk and shear moduli

$K_{a}=\frac{C}{6}(1-\phi_{0})\left(K_{2}+\frac{4}{3}\mu_{2}\right)S_{n}\ \ \
{\rm and}\ \ \
\mu_{a}=\frac{3}{5}\left[K_{a}+\frac{1}{4}{C}(1-\phi_{0})\mu_{2}S_{t}\right],$
(15)

respectively, where ${C}$ is the number of contact per grain and $S_{n}$ and
$S_{t}$ are related to the normal and shear stiffnesses, respectively, of a
cemented two-grain combination. The explicit expressions of these quantities
are given in Mavko et al. (2009, p. 256) and depend on the Poisson ratio of
the grains $\nu_{1}=(3K_{1}-2\mu_{1})/(6K_{1}+2\mu_{1})$, the shear modulus
$\mu_{1}$, and on the ratio of the radius of the cement layer to the grain
radius:

$\alpha=2\left[\frac{\bar{\phi}_{p}}{3{C}(1-\phi_{0})}\right]^{1/4},$ (16)

where we consider here that

$\bar{\phi}_{p}=\left\\{\begin{array}[]{ll}\phi_{p},&\ \ \
\phi_{p}<\phi_{p0},\\\ \phi_{p0},&\ \ \
\phi_{p}\geq\phi_{p0},\end{array}\right.$ (17)

where $\phi_{p}$ is given by equation (B.1) and $\phi_{p0}$ =0.05 is the
maximum amount of cement at the grain contacts (see below and Figure 5). Above
this value, the cement is free in the pore space with fraction
$\phi_{p}-\phi_{p0}$. Alternatively, for a uniform distribution of the cement
on the grain surface,

$\alpha=\sqrt{\frac{2\bar{\phi}_{p}}{3(1-\phi_{0})}}.$ (18)

A unified equation for $\alpha$ can be obtained by introducing a new parameter
that allows to interpolate between the two cases represented by equations (16)
and (18) (Allo, 2019; Table 1), but as we shall see in the examples, the
velocities obtained with these two cases do not differ significantly from a
practical viewpoint.

However, the CCT equation (15) holds for high porosity, small amounts of
cement (at the grain contacts) and does not consider the effect of pressure.
Here, we model pressure effects with a modified version of the patchy cement
model of Avseth et al. (2016) and include the compaction effect given by
equation (B.7). The high-porosity limit moduli, $K_{b}$ and $\mu_{b}$, are
obtained by “mixing” the CCT ($K_{a}$-$\mu_{a}$) and Hertz-Mindlin (HM)
($K_{u}$-$\mu_{u}$) uncemented moduli (see Appendix D and Figure 6) with the
HS upper bound (Appendix C), where the volume fraction of cemented rock to be
used in the HS bound is $f=\phi_{p}/\phi_{0}$, i.e., $f$ = 0 if there is no
cement and 1 if the whole pore-space is filled with cement.

Next, we apply the SC theory to obtain the moduli $K_{s}$ and $\mu_{s}$ of a
medium where all the pores are filled with cement. There are two phases, grain
and cement (as spheres). To this purpose, we use equation (8) with $n$ = 2,
$\phi_{1}=1-\phi_{0}$, $\phi_{2}=\phi_{0}$, $K_{1}=K_{1}$, $K_{2}=K_{2}$,
$\mu_{1}=\mu_{1}$, $\mu_{2}=\mu_{2}$, $K=K_{s}$ and $\mu=\mu_{s}$.

Finally, we interpolate between the effective high-porosity end member given
by the HS upper bound and the mineral point (i.e., zero porosity) using the
Voigt-Reuss-Hill average, i.e., an arithmetic average of the Voigt and Reuss
moduli, to obtain the dry-rock bulk and shear moduli:

$K_{m}=\frac{1}{2}(K_{V}+K_{W})\ \ \ {\rm and}\ \ \
\mu_{m}=\frac{1}{2}(\mu_{V}+\mu_{W}),$ (19)

where

$K_{V}=(1-\phi/\phi_{0})K_{s}+(\phi/\phi_{0})K_{b},\ \ \ \
\frac{1}{K_{W}}=\frac{1-\phi/\phi_{0}}{K_{s}}+\frac{\phi/\phi_{0}}{K_{b}},$
(20) $\mu_{V}=(1-\phi/\phi_{0})\mu_{s}+(\phi/\phi_{0})\mu_{b}\ \ \ {\rm and}\
\ \
\frac{1}{\mu_{W}}=\frac{1-\phi/\phi_{0}}{\mu_{s}}+\frac{\phi/\phi_{0}}{\mu_{b}}.$
(21)

The Hill average is close in accuracy to more sophisticated techniques such as
self-consistent schemes and are applicable to complex rheologies such as
general anisotropy and arbitrary grain topologies (e.g., Man and Huang, 2011)
and anelasticity (Picotti et al., 2018; Qadrouh et al., 2020).

The wet-rock moduli are obtained with Gassmann equations (4).

##  5 Seismic velocities

The P-wave modulus and velocity are

$E=K_{G}+\frac{4}{3}\mu_{m}.$ (22)

and

$V_{P}=\sqrt{\frac{E}{\rho}},$ (23)

respectively, where $\rho$ is the composite density, given by

$\rho=(1-\phi)\sum_{i=1}^{n}\beta_{i}\rho_{i}+\phi\rho_{f},$ (24)

where $\rho_{i}$ and $\rho_{f}$ are the densities of the $i$-th solid phase
and fluid, respectively.

Similarly, the S-wave velocity is

$V_{S}=\sqrt{\frac{\mu}{\rho}},$ (25)

##  6 Examples

The elastic properties of the minerals are given in Table 1. We assume $G$ =
30 oC/km, $S$= 0.04 km/m.y, which give a heating rate $H$ = 1.2 oC/m.y. and we
take $T_{0}$ = 15 oC at $z$ = 0. Figure 7 shows the relation between depth
temperature and pressure, corresponding to this simple (linear) basin
modeling, given by equations (1) and (2).

###  6.1 Shales

We consider an initial smectite/illite ratio $r_{0}$ = 0.99 and $n$ = 5. The
kinetic reaction corresponding to smectite/illite conversion assumes $E$ = 36
kcal/mol and $c$ = 1.217 $\times$ 1023/ m.y. (Pytte and Reynolds, 1989).
Figure 8 shows the effect of the geothermal gradient on the illite-smectite
fraction, indicating that a higher value accelerates the conversion. The
conversion ratio and density (a) and P- and S- wave velocities (b) of the
mineral mixture, as a function of depth, are shown in Figure 9, where, in this
case $G$ = 30 oC/km. As can be seen, the elastic properties increases with
depth due to the higher stiffness and density of illite compared to smectite.

Next, we consider that at $z$ = 0, the initial porosity is $\phi_{0}$ = 0.35,
and $\beta$ = 0.01/MPa in Athy equation (7). Then, we compute the porosity,
density and wave velocities for Models 1, 2 and 3 as a function of depth,
where $\gamma$ = 0.17 in Model 2 (aspect ratio of the pores). Figure 10 shows
the porosity and density, which have the same values for the three models,
whereas Figure 11 displays the P-wave (a) and S-wave (b) velocities, where the
three models yield similar results in practice. The results of Model 2 differ
at lower depths, whose theory is based on the assumption of idealized
geometries of the grains and pores. If we consider spheres, i.e., $\gamma$ =
1, the P and S-wave velocities predicted by Model 2 are higher by
approximately 0.5 km/s than those of Models 1 and 3.

To illustrate the effect of the geothermal gradient (and heating rate) on the
velocities, we compare the results from two values, $G$ = 20 and 30 oC/km ($H$
= 0.8 and 1.2 oC/m.y.), keeping the same sedimentation rate. Figure 12 shows
the velocities as a function of depth, where the black curves correspond to
the higher value of the heating rate (more illite implies higher velocities).
The lower velocities for 20 oC/km are due to the higher amount of smectite.

###  6.2 Sandstones

The dissolution and precipitation of dissolved silica are the source of
cement, which is quantified by equation (B.1) at geological time. The porosity
is given by equation (B.8), taken into account the pressure effects on
compaction.

We assume the same geothermal gradient, sedimentation rate and heating rate of
the first examples, so that Figure 7 describes the basin modeling in this
case. Let us consider the following properties: $\phi_{0}$ = 0.36, $D$ = 0.03
cm, $a$ = 1.98 $\times$ 10-22 mol/(cm2 s), $b$ = 0.022 ∘C-1, $\rho_{1}$ = 2.65
g/cm3, $M_{q}$ = 60.09 g/mol, $f$= 0.65, and $V$ = 1 cm3, $\beta$ = 0.01 MPa-1
and IGVi = 0.2. The time step to solve equation (B.1) is $dt$ = 0.1 m.y.
Figure 13 shows the cement fraction and porosity, with and without compaction.
In the latter case, the theory holds down to approximately 4.5 km depth, below
which the porosity becomes negative and the whole pore space is filled with
cement.

Let us now consider two values of $f$, the amount of detrital quartz in the
rock, i.e., $f$ = 0.4 (0.6 feldspar) and $f$ = 1. Figure 14 indicates that
nonquartz grains inhibit the cementation and the porosity loss. Using equation
(B.2) indicates that the exponent in (B.1) is proportional to $f/D$ so that
the effect of increasing the grain diameter is similar to that of decreasing
the detrital quartz content, i.e., the amount of precipitated quartz cement is
less in coarse-grain sandstones, because of the reduced surface area, compared
to fine-grain sandstones.

The cementation only depends on the heating rate $H=GS$, so that different
combinations of the geothermal gradient and sedimentation rate, but keeping
the same heating rate, yield the same fraction of cement. However, it is not
clear from equation (B.1) how $H$ affects the cementation, since it appears in
several terms of the exponent. We consider three extreme values of $H$ = 0.4,
1.2 and 2 oC/m.y. Figure 15 indicates that $H$ greatly affects the process,
where we assumed $f$ = 0.65 and $G$ = 30 oC/km. The plot shows that, for any
given depth, the greatest fraction of cement in the sandstone that had the
most time to form is that corresponding to a heating rate of 0.4 oC/m.y.,
which at 3 km depth would be 225 m.y. old and would have been above say, 80oC,
for 25 m.y.. On the contrary, at 3 km depth the sandstone that experienced a
heating rate of 2 oC/m.y. would be only 45 m.y. old and would have been above
80oC for only 5 m.y..

Next, we obtain the wave velocities. To compute the mineral HS averages in
Models 1 and 3, the proportions of quartz (grain) and cement are
$(1-\phi-\phi_{p})/(1-\phi)$ and $\phi_{p}/(1-\phi)$, respectively. The three
phases of Model 2 have proportions $1-\phi_{p}-\phi$ (grain), $\phi_{p}$
(cement) and $\phi$ (fluid or pore). The pore-aspect ratio in Model 2 is taken
$\gamma$ = 0.108 and the bonding parameter in Model 4 is assumed to be
$\delta$ = 100 MPa. Figure 16 shows the P- and S-wave velocities of the
sandstones as a function of depth, corresponding to the four petro-elastical
models. As can be observed, all the models are quite consistent with similar
trends and values of the velocities, with Models 1 and 3 giving almost
identical results. The research indicates that the use of simple models (e.g.,
Model 1) can be used in combination with the diagenesis-cementation process to
make predictions.

However, Models 2 and 4 have additional parameters, i.e., the pore aspect
ratio ($\gamma$) (Model 2), and the maximum contact cement $\phi_{p0}$ and
bonding parameters $\delta$ (Model 4). It can be shown that taking $\phi_{p0}$
= 0.1 (instead of 0.05) does not affect noticeably the results. On the other
hand, the other two parameters ($\gamma$ and $\delta$) have an effect, as can
be seen in Figure 17, where the velocities of Models 2 and 4 are also
represented for $\gamma$ = 1 (spherical pores) and $\delta$ = 40 MPa (Gangi
and Carlson, 1996). While the aspect ratio highly affects the velocity (e.g.,
Cheng et al., 2020), the initial bonding of the grains has a lesser influence.

Next, we consider the more simple and complex models (1 and 4, respectively)
and show the P-wave velocities for two values of the geothermal gradient $G$ =
20 and 30 oC/km ($H$ = 0.8 and 1.2 oC/m.y.) keeping the same sedimentation
rate. The results are shown in Figure 18, where, as expected, a smaller
temperature gradient (and heating rate) implies less diagenesis and
cementation and lower velocities. Moreover, the two models predict the same
trend and values of the velocities, indicating that the more simple Model 1 is
suitable for predictions as Model 4, which is more complex from a mathematical
point of view.

Model 4 allows to change the location of contact cement from the grain contact
to a uniform distribution on the grain surface. In this case, we use equation
(18) and assume $G$ = 30 oC/km, $\phi_{p0}$ = 0.05 and $\delta$ = 100 MPa.
Figure 19 shows the velocities when the bonding cement is deposited at grain
contacts and uniformly distributed. The differences are minimal in practice.

The fact that the properties of the cement are close to those of quartz in
Model 4 does not affect much the results. The CCT theory predicts that the
type of cement has not a relevant influence on the the dry-rock moduli.
Varying the cement stiffness one order of magnitude results in only an about
15 % increase in velocities, an effect much smaller than that of increasing
cement content [Fig. 11 in Dvorkin et al. (1994)].

##  7 Extensions of the models

The smectite/illite petro-elastical model can be extended to the anisotropic
case, as in Carcione and Avseth (2015), and to the presence of overpressure by
using the pore-volume balance equations developed by Carcione and Gangi (2000)
for disequilibrium compaction, theories also implemented by Qin et al. (2019).
Wangen (2000) developed an alternative theory of overpressure by cementation
of sediments. There is no overpressure if the released water caused by the
smectite/illite reaction mostly escapes from the pore space (due in part to a
relatively slow sedimentation rate), which is the case in the present work.

Another type of diagenesis process that can easily be incorporated into the
theory is the conversion of solid organic matter to hydrocarbons (Luo and
Vasseur, 1996; Carcione et al., 2011), as mentioned in the introduction.
Attenuation and velocity dispersion can be modeled with a theory based on
physical principles (not developed yet, to our knowledge) or
phenomenologically with the realistic assumption that cementation reduces
seismic losses (Gei and Carcione, 2003), since attenuation decreases with
depth and compaction. A completely different approach is the molecular dynamic
modeling (Garcia and Medina, 2007), where at each step of the algorithm,
cement is added at specific locations within the pores, in three different
ways that model distinct origins and microgeometric features.

##  8 Conclusions

We have developed a procedure to combine diagenesis and cementation processes
for the formation of shales and sandstones from non-consolidated clay and
sand, respectively, to different petro-elastical models to obtain the seismic
velocities as a function of depth, pressure and temperature, on the basis of a
linear basin modeling. The amount of illite in shales and quartz cementation
in sandstones and resulting porosity loss and compaction affect the
velocities, obtained with Gassmann-like equations, models based on oblate
spheroid inclusions, and theories that assume a random packing of identical
spherical grains (granular media).

The proposed methods take into account in-situ properties such as the
geothermal gradient, sedimentation and heating rates, elastic properties of
the diagenetic material (cement), cement distribution in the pore space, pore
aspect ratios, effective pressure and initial pre-stress conditions. All the
models, based on different physical assumptions, predict the same increasing
trend of velocity with burial history and similar values, indicating the
robustness of the different theories, either based on the critical-porosity
concept (e.g., Krief model) or an idealized shape of pores and grains to
obtain the dry-rock stiffness moduli.

## REFERENCES

> Allo, F., 2019, Consolidating rock-physics classics: A practical take on
> granular effective medium, The Leading Edge, 38(5), 334–340.
>
> Athy, L. F., 1930, Density, porosity, and compaction of sedimentary rocks,
> AAPG Bulletin, 14, 1–24
>
> Avseth, P., and Lehocki, I., 2016, Combining burial history and rock-physics
> modeling to constrain AVO analysis during exploration, The Leading Edge,
> 35(6), 528-534.
>
> Avseth, P., Mukerji, T., Mavko, G., and Dvorkin, J., 2010, Rock-physics
> diagnostics of depositional texture, diagenetic alterations, and reservoir
> heterogeneity in high-porosity siliciclastic sediments and rocks – A review
> of selected models and suggested work flows, Geophysics, 75, 75A31–75A47.
>
> Avseth, P., Skjei, N., and Mavko, G., 2016, Rock-physics modeling of stress
> sensitivity and 4D time shifts in patchy cemented sandstones – Application
> to the Visund Field, North Sea, The Leading Edge, 35(10), 868–878.
>
> Berryman , J. G., 1980, Long-wavelength propagation in composite elastic
> media II. Ellipsoidal inclusions, J. Acoust. Soc. Am., 68, 1820–1831.
>
> Carcione, J. M., 2014, Wave Fields in Real Media. Theory and numerical
> simulation of wave propagation in anisotropic, anelastic, porous and
> electromagnetic media, 3rd edition, revised and extended, Elsevier Science.
>
> Carcione, J. M., and Avseth, P., 2015, Rock-physics templates for clay-rich
> source rocks, Geophysics, 80, D481–D500.
>
> Carcione, J. M., and Gangi, A., 2000, Non-equilibrium compaction and
> abnormal pore-fluid pressures: effects on seismic attributes, Geophys.
> Prosp., 48, 521–537.
>
> Carcione, J. M., Farina, B., Poletto, F., Qadrouh, A. N., and Cheng, W.,
> 2020, Seismic attenuation in partially molten rocks, Physics of the Earth
> and Planetary Interiors, https://doi.org/10.1016/j.pepi.2020.106568
>
> Carcione, J. M., Helle, H. B., and Avseth, P., 2011, Source-rock seismic-
> velocity models, Gassmann versus Backus: Geophysics, 76, N37-N45.
>
> Carcione, J. M., Helle, H. B., Santos, J. E., and Ravazzoli, C. L., 2005, A
> constitutive equations and generalized Gassmann modulus for multi-mineral
> porous media, Geophysics, 70, N17-N26.
>
> Carcione, J. M., Landrø, M., Gangi, A. F., and Cavallini, F., 2007,
> Determining the dilation factor in 4D monitoring of compacting reservoirs by
> rock-physics models, Geophys. Prosp., 55, 793–804.
>
> Carcione, J. M., Poletto, F., and Farina, B., 2018, The Burgers/squirt-flow
> seismic model of the crust and mantle, Physics of the Earth and Planetary
> Interiors, 274, 14–22.
>
> Cheng, W., Carcione, J. M., Qadrouh, A., Alajmi, M., and Ba, J., 2020, Rock
> anelasticity, pore geometry and the Biot-Gardner effect, Rock Mechanics and
> Rock Engineering, https://doi.org/10.1007/s00603-020-02155-7.
>
> Draege, A., Jakobsen, M., and Johansen, T. A., 2006, Rock physics modelling
> of shale diagenesis, Petroleum Geoscience, 12, 49–57.
>
> Dvorkin, J., Berryman, J., and Nur, A, 1999, Elastic moduli of cemented
> sphere packs, Mechanics of Materials, 31, 461–469.
>
> Dvorkin, J., Nur, A., and Yin, H., 1994, Effective properties of cemented
> granular material, Mechanics of Materials 18, 351–366.
>
> Dvorkin, J., and Nur, A., 1996, Elasticity of high-porosity sandstones:
> Theory for two North Sea data sets, Geophysics, 61, 1363–1370.
>
> Gangi, A. F., and Carlson, R. L., 1996, An asperity-deformation model for
> effective pressure,Tectonophysics, 256, 241–251.
>
> Garcia, X., and Medina, E., 2007, Acoustic response of cemented granular
> sedimentary rocks: Molecular dynamics modeling, Physical Review E, 75,
> 061308 2007.
>
> Gei, D., and J. M. Carcione, 2003, Acoustic properties of sediments
> saturated with gas hydrate, free gas and water: Geophys. Prosp., 51,
> 141-157.
>
> Goffe, W. L., Ferrier, G. D., and Rogers, J., 1994, Global optimization of
> statistical functions with simulated annealing, Journal of Econometrics,
> 60(1-2), 65–99.
>
> Gurevich, B. and Carcione, J. M., 2000, Gassmann modeling of acoustic
> properties of sand/clay mixtures, Pure and Applied Geophysics, 157, 811–827.
>
> Hashin, Z., and Shtrikman, S., 1963, A variational approach to the theory of
> the elastic behaviour of multiphase materials, Journal of the Mechanics and
> Physics of Solids, 11, 127–140.
>
> Hertz, H., 1895, Gesammelte Werke, Band I, Schriften Vermischten Inhalts,
> Leipzig.
>
> Keller, W. D., Reynolds, R. C., and Inoue, A., 1986, Morphology of clay
> minerals in the smectite to illite conversion series by scanning electron
> microscopy, Clays & Clay Minerals, 34, 187–197.
>
> Krief, M., J. Garat, J. Stellingwerff, and J. Ventre, 1990, A petrophysical
> interpretation using the velocities of P and S waves (full waveform sonic):
> The log Analyst, 31, 355-369.
>
> Lander, R. H., and O. Walderhaug, 1999, Predicting porosity through
> simulating sandstone compaction and quartz cementation, AAPG Bulletin, 83,
> 433–449.
>
> Luo, X., Vasseur, G., 1996, Geopressuring mechanism of organic matter
> cracking: Numerical modeling, AAPG Bulletin, 80, 856–874
>
> Man, C.-S., and M. Huang, 2011, A simple explicit formula for the Voigt-
> Reuss-Hill average of elastic polycrystals with arbitrary crystal and
> texture symmetries, Journal of Elasticity, 105, 29–48.
>
> Mavko, G., T. Mukerji, and J. Dvorkin, 2009, The rock physics handbook:
> tools for seismic analysis in porous media, Cambridge Univ. Press.
>
> Mindlin, R. D., 1949, Compliance of elastic bodies in contact, J. Appl.
> Mech., 16, 259-268.
>
> Mondol, N. H., J. Jahren, K. Bjørlykke, and I. Brevik, 2008, Elastic
> properties of clay minerals, The Leading Edge, 27(6), 758-770.
>
> Paxton, S. T., Szabo, J. O., Ajdukiewicz, J. M., and Klimentidis, R. E.,
> 2002, Construction of an intergranular volume compaction curve for
> evaluating and predicting compaction and porosity loss in rigid-grain
> sandstone reservoirs, AAPG Bulletin, 86, 2047–2067.
>
> Perry, E., and Hower, J., 1970, Burial diagenesis in Gulf Coast pelitic
> sediments, Clays & Clay Minerals 18, 165–177
>
> Picotti, S., Carcione, J. M., and Ba, J., 2018, Rock-physics templates based
> on seismic $Q$, Geophysics, 84, MR13–MR23.
>
> Pinna, G., J. M. Carcione, and F. Poletto, 2011, Kerogen to oil conversion
> in source rocks. Pore-pressure build-up and effects on seismic velocities:
> J. Appl. Geophy., 74, 229–235.
>
> Pytte, A. M., and R. C. Reynolds, 1989, The thermal transformation of
> smectite to illite. In N. D. Naeser and T. H. McCulloh, editors, Thermal
> History of Sedimentary Basins: Methods and Case Histories, pages 133–140.
> Springer–Verlag.
>
> Qadrouh, A. N., Carcione, J. M., Alajmi, M., and Ba, J., 2020, Bounds and
> averages of seismic $Q$, Stud. Geophys. Geod.,
> https://doi.org/10.1007/s11200-019-1247-y.
>
> Qin, X., Han, D.-H., and Zhao, L., 2019, Elastic characteristics of
> overpressure due to smectite-to-illite transition based on micro-mechanism
> analysis, Geophysics, 84, 1JA-Z21.
>
> Rieke III, H. H., and Chilingarian, G. V., 1974, Compaction of argillaceous
> sediments. Elsevier, New York.
>
> Sadikh-Zadeh, L. A., 2006, Prediction of sandstone porosity through
> quantitative estimation of its mechanical compaction during lithogenesis
> (Bibieybat Field, South Caspian Basin), Natural Resources Research, 15, DOI:
> 10.1007/s11053-006-9008-3
>
> Scotchman, I. C., 1987, Clay diagenesis in the Kimmeridge Clay Formation,
> onshore UK, and its relation to organic maturation: Mineral. Mag., 51,
> 535-551.
>
> Sippel, R. F., 1968, Sandstone petrology, evidence from luminescence
> petrography, Journal of Sedimentary Petrology, 38, 530–554.
>
> Stephenson, L. P., Plumley, W. J., and Palciauskas, V. V., 1992, A model for
> sandstone compaction by grain interpenetration: Journal of Sedimentary
> Petrology, 62, 11–22.
>
> Thyberg, B. I., J. Jahren, T. Winje, K. Bjørlykke, and J. I. Faleide, 2009.
> From mud to shale: rock stiffening by micro-quartz cementation: First Break,
> 27, 27-33.
>
> Walderhaug, O.,1996, Kinetic modelling of quartz cementation and porosity
> loss in deeply buried sandstone reservoirs. AAPG Bulletin, 80, 731–745.
>
> Walderhaug, O., 2000, Modeling quartz cementation and porosity in Middle
> Jurassic Brent Group sandstones of the Kvitebjørn field, northern North Sea,
> AAPG Bulletin, 84, 1325–1339.
>
> Walton, K., 1987, The effective elastic moduli of a random packing of
> spheres, J. Mech. Phys. Solids, 35, 213–226.
>
> Wangen, M., 2000, Generation of overpressure by cementation of pore space in
> sedimentary rock, Geophys. J. Int., 143, 608–620.
>
> Weyl, P. K., 1959, Pressure solution and the force of crystallization – a
> phenomenological theory, J. Geophys. Res., 64, 2001–2025.

## Statements and Declarations

The authors declare that no funds, grants or other support were received
during the preparation of this manuscript, and have no relevant financial or
non-financial interests to disclose.

All authors contributed to the study conception and design. The manuscript was
written by J. M. Carcione, and D. Gei and S. Picotti checked the shale and
sandstone sections, respectively. A. N. Qadrouh, M. Alajmi and J. Ba
contributed to the review of the literature, data collection and analysis. All
authors commented on previous versions of the manuscript and approved the
final manuscript.

##  Appendix A Smectite/illite conversion kinetic model

The transformation of smectite to illite is probably the most important clay
mineral reaction in sedimentary rocks (Perry and Hower, 1970). Pytte and
Reynolds (1989) propose a model for the smectite/illite ratio $r$ based on the
$n$th-order Arrhenius-type reaction

$\frac{\partial r}{\partial t}=-r^{n}c\exp[-E/RT(t)],$ (A.1)

where $E$ is the activation energy, $R$ = 1.986 cal/( mol oK ) is the gas
constant, $c$ is a constant and $T(t)$ is the absolute temperature.

The illite/smectite ratio in percent is 100 ($1-r$). The solution of equation
A.1 is

$r(t)=m^{-1/m}\left\\{\frac{r_{0}^{-m}}{m}+\frac{c}{H}\left[\frac{E}{R}\left[{\rm
Ei}(x)-{\rm
Ei}(x_{0})\right]+T\exp(x)-T_{0}\exp(x_{0})\right]\right\\}^{-1/m},$ (A.2)

where $r_{0}$ is the initial ratio, $m=n-1$, Ei ($x$) is the exponential
integral,

$x=-\frac{E}{RT},\ \ \ x_{0}=-\frac{E}{RT_{0}},$ (A.3)

and the dependence on the deposition time is given in the absolute temperature
[see equation (2)]. To compute (A.2), we use the property Ei($x$) = $-$ ${\rm
E}_{1}(-x)$.

##  Appendix B Quartz cement precipitated in sandstones and compaction

It is assumed that dissolved silica is the source of cement (quartz
dissolution) at stylolites or at grain contacts containing clay or mica, and
diffuses short distances to the sandstone volume, so that no quartz
dissolution or grain interpenetration takes place within this volume.
Following Walderhaug (1996), the amount of quartz cement (cm3) precipitated in
a 1 cm3 of sandstone at time $t_{n+1}$ is

$\phi_{p}^{n+1}=\phi_{0}-(\phi_{0}-\phi_{p}^{n})\exp\left[-\frac{M_{q}aA_{0}}{\rho_{q}\phi_{0}bH\ln
10}\left(10^{bT_{n+1}}-10^{bT_{n}}\right)\right],\ \ \ T_{n}=Ht_{n}+T_{0},$
(B.1)

where $\phi_{p}^{n}$ is the amount of quartz cement present at time $t_{n}$,
$\phi_{0}$ is the initial porosity, $a$ and $b$ are constants, which have
units of mol/(cm2 s) and ∘C-1, respectively, $H$ is the heating rate,
$\rho_{q}$ is the quartz density, $M_{q}$ is the molar mass of quartz, and

$A_{0}=\frac{6fV}{D}$ (B.2)

is the initial quartz surface area, where $D$ is the grain diameter, $f$ is
the fraction of detrital quartz in the rock and $V$ is a unit volume (1 m3 if
$D$ is given in meters). Walderhaug (1996) shows examples with 65 and 50 %
detrital quartzite and the rest is feldspar.

In the first time step, $t_{0}$ = 0, $\phi_{p}^{0}$ = 0, $T_{1}=Hdt+T_{0}$ and

$\frac{\phi_{p}^{1}}{\phi_{0}}=1-\exp\left[-\frac{M_{q}aA_{0}}{\rho\phi_{0}bH\ln
10}(10^{bT_{1}}-10^{bT_{0}})\right],$ (B.3)

where $dt$ is the time step.

The porosity varies as

$\phi=\phi_{0}-\phi_{p},$ (B.4)

and the surface area as

$A=A_{0}\frac{\phi}{\phi_{0}}.$ (B.5)

The effect of clay coating is to reduce the precipitation and is modeled as

$A_{0}\rightarrow(1-{\cal C})A_{0},$ (B.6)

where ${\cal C}$ is a clay factor. ${\cal C}$ = 1 implies zero surface area
and no precipitation. The clay factor affects the surface area in the same
manner as the fraction of detrital quartz $f$ and reciprocally with respect to
the grain diameter $D$.

Equation (B.1) assumes a constant heating rate. An algorithm for a variable
heating rate is gieven by Eq. (5) in Walderhaug (1996).

Porosity reduction by compaction is based on the intergranular volume loss
index (IGV)

$\mbox{IGV}=\mbox{IGVi}+(\phi_{0}+m_{0}-{\rm IGVi})\exp(-\beta p_{e})$ (B.7)

(Lander and Walderhaug, 1999), where IGVi is the value at infinite effective
pressure, $m_{0}$ is the initial matrix volume fraction in the pore space, and
$\beta$ is a constant. If $p_{e}$ = 0, IGV = $\phi_{0}+m_{0}$, the maximum
value. For IGVi = $m_{0}$ = 0, we obtain Athy equation (7), while for IGVi =
$\phi_{0}$ and $m_{0}$ = 0, we obtain IGV = $\phi_{0}$. A detailed pictorial
explanation of the IGV index can be found in Paxton et al. (2002). IGVi and
$\beta$ are obtained by fitting experimental data.

In the absence of matrix material, the porosity is given by

$\phi=\mbox{IGV}-\phi_{p}$ (B.8)

(e.g., Sadikh-Zadeh, 2006). Equation (B.4) is obtained for $p_{e}$ = 0,
neglecting the pressure effects that lead to compaction. Actually, there is no
particular theory behind equation (B.7), which is an empirical equation of the
type used by Athy (1930) ninety years ago, of the form
$\phi=\phi_{0}\exp(-\beta^{\prime}z)$, where
$\beta^{\prime}=g(\bar{\rho}-\bar{\rho}_{w})\beta$.

##  Appendix C Hashin and Shtrikman bounds and averages

Let us denote the solid bulk and shear moduli by $K_{i}$ and $\mu_{i}$,
respectively. A two-solid composite, with no restriction on the shape of the
two phases, has stiffness bounds given by the Hashin and Shtrikman (1963)
equations,

$K_{\rm
HS}^{\pm}=K_{1}+\frac{\beta_{2}}{(K_{2}-K_{1})^{-1}+\beta_{1}\left(K_{1}+\displaystyle\frac{4}{3}\mu_{\beta}\right)^{-1}}$
(C.1)

and

$\mu_{\rm
HS}^{\pm}=\mu_{1}+\frac{\beta_{2}}{(\mu_{2}-\mu_{1})^{-1}+\beta_{1}\left[\mu_{1}+\displaystyle\frac{\mu_{\beta}}{6}\left(\displaystyle\frac{9K_{\beta}+8\mu_{\beta}}{K_{\beta}+2\mu_{\beta}}\right)\right]^{-1}},$
(C.2)

where $\beta_{1}$ and $\beta_{2}$ are the fractions of solid 1 and 2
($\beta_{1}+\beta_{2}$ = 1), and $\beta_{1}=\phi_{1}/(1-\phi)$. We obtain the
upper bounds when $K_{\beta}$ and $\mu_{\beta}$ are the maximum bulk and shear
moduli of the single components, and the lower bounds when these quantities
are the corresponding minimum moduli, i.e., we have the upper bound if 1 is
the stiffer medium and the lower bound is obtained if 1 is the softer medium
(Mavko et al., 2009).

The arithmetic averages of the bounds are frequently used to obtain the bulk
and shear moduli of a mineral mixture, i.e.,

$K_{s}=\frac{1}{2}(K_{\rm HS}^{+}+K_{\rm HS}^{-}),\ \ \
\mu_{s}=\frac{1}{2}(\mu_{\rm HS}^{+}+\mu_{\rm HS}^{-}).\ \ \ $ (C.3)

##  Appendix D The Hertz-Mindlin model

We model pressure dependence in Model 4 by using the Hertz-Mindlin (HM)
contact theory, which considers spherical grains. The stresses are calculated
in terms of the strains by considering the random packing of spheres as an
effective medium that exerts a mean-field force (as given by the contact
Hertzian theory) on a single representative grain (Hertz, 1895; Mindlin, 1949;
Walton, 1987; Mavko et al., 2009).

We modify the Hertz-Mindlin model by replacing the effective pressure $p_{e}$,
by the augmented value, $p_{e}+\delta$, following Gangi and Carlson (1996),
assuming that at the beginning of the burial process the grains are not
“floating” in the fluid, but there is an initial level of bonding between
grains determined by $\delta$ (a pre-stress condition). A HM model with
augmented effective pressure is given in Carcione et al. (2007) to model
bonded grains at $p_{e}$ = 0.

Then, the bulk and shear (uncemented) moduli at the critical porosity
$\phi_{0}$ are given by

$K_{u}=\left[\frac{{C}^{2}(1-\phi_{0})^{2}\mu_{1}^{2}(p_{e}+\delta)}{18\pi^{2}(1-\nu_{1})^{2}}\right]^{1/3}$
(D.1)

and

$\mu_{u}=\frac{3(5-4\nu_{1})}{5(2-\nu_{1})}K_{u},$ (D.2)

where $\mu_{1}$ is the shear modulus of the grains, $\nu_{1}$ is the Poisson
ratio of the grains and ${C}$ is the average number of contacts per spherical
grain.

##  5 Tables

Table 1. Material Properties.

property | smectite | illite | quartz | cement | water
---|---|---|---|---|---
$K$ (GPa) | 9 | 33 | 40 | 35 | 2.2
$\mu$ (GPa) | 6 | 28 | 39 | 30 | 0
$\rho$ (g/cm3) | 2.25 | 2.8 | 2.65 | 2.60 | 1.04

Figure 1:  Scanning electron micrographs (SEM) of a shale showing a texture
corresponding to a smectite/illite mixture (Sichuan Basin, China; depth: 3.3
km).  Figure 2:  Petrographic thin section of a sandstone showing the mineral
composition (Sichuan Basin, China; depth: 2 km).  Figure 3:  Clay diagenesis.
Smectite/illite conversion with release of bound water. Clays consist of bound
water at the time of deposition. Free water increases with burial with the
consequent release of bound water. Then, free water is squeezed out of the
original volume (modified from Rieke III and Chilingarian, 1974, Fig. 57).
Figure 4:  Sandstone diagenesis and cementation. Compaction reduces quartz
surface area by increasing grain contact area, as well as by the injection of
matrix material into the pore space. Cementation causes surface area reduction
when quartz grains are encased by pore-filling cements (modified from Lander
and Walderhaug, 1999).  Figure 5:  Model 4. Sandstone after cementation:
$\phi_{p0}$ is the maximum cement fraction bonded to the grains, $\phi_{p}$ is
the total volume fraction of cement, and $\phi_{p}-\phi_{p0}$ is the pore-
filling cement fraction. Figure 6:  A random packing of spheres corresponding
to the Hertz-Mindlin model, where the contact area between the grains depends
on the effective pressure. Figure 7:  Temperature-pressure-depth relation of
the linear basin modeling. Figure 8:  Illite-smectite fraction as a function
of depth for three values of the geothermal gradient. Figure 9:  Illite-
smectite fraction and density (a), and wave velocities of the mineral mixture
as a function of depth. Figure 10:  Shale porosity and bulk density as a
function of depth. Figure 11:  P-wave (a) and S-wave (b) velocities of the
shale as a function of depth, corresponding to Models 1, 2 and 3. Model 2
considers a pore aspect ratio $\gamma$ = 0.17. Figure 12:  P-wave velocities
of the shale as a function of depth, corresponding to Models 1, 2 and 3 for
two values of the geothermal gradient. Figure 13:  Cement fraction $\phi_{p}$,
porosity (no compaction) and porosity (with compaction) as a function of depth
of a sandstone body subject to precipitation and cementation. Figure 14:
Cement fraction $\phi_{p}$ and porosity $\phi$ as a function of depth for two
values of detrital quart in the sandstone $f$ = 0.4 and $f$ = 1. Figure 15:
Cement fraction $\phi_{p}$ as a function of depth for three values of the
heating rate. Figure 16:  P- and S-wave velocities of the sandstone as a
function of depth, corresponding to the four petro-elastical models. Figure
17:  P- and S-wave velocities of the sandstone as a function of depth,
corresponding to Model 2 (a) and Model 4 (b. for two values of the pore-aspect
ratio $\gamma$ and bonding parameter $\delta$, respectively. Figure 18:
P-wave velocities of the sandstone as a function of depth, corresponding to
Models 1 and 4 for two values of the geothermal gradient. Figure 19:  P- and
S-wave velocities of the sandstone as a function of depth, corresponding to
Model 4 and two different distributions of the contact cement, based on
equations (16) and (18) (cement at grain contacts and coating the grains,
respectively).


# Permanence and uniform asymptotic stability of positive solutions of SAIQH
models on time scales

Nedjoua Zine  University Mustapha Stambouli of Mascara, Algeria
<EMAIL_ADDRESS>,  Benaoumeur Bayour  University Mustapha
Stambouli of Mascara, Algeria<EMAIL_ADDRESS>and  Delfim F. M.
Torres  Center for Research and Development in Mathematics and Applications
(CIDMA), Department of Mathematics, University of Aveiro, 3810-193 Aveiro,
Portugal<EMAIL_ADDRESS>Research Center in Exact Sciences (CICE), Faculty of
Sciences and Technology (FCT), University of Cape Verde (Uni-CV), 7943-010
Praia, Cape Verde<EMAIL_ADDRESS>
(Date: Submitted January 31, 2024)

###### Abstract.

A Susceptible–Asymptomatic–Infectious–Quarantined–Hospitalized (SAIQH)
compartmental model on time scales is introduced and a suitable Lyapunov
function is defined. Main results include: the proof that the system is
permanent; proof of existence of solution; and sufficient conditions implying
the dynamic system to have a unique almost periodic solution that is uniformly
asymptotically stable. An example is presented supporting the obtained
results.

###### Key words and phrases:

SAIQH models; time scales; permanence; almost periodic solutions; uniform
asymptotic stability.

###### 2020 Mathematics Subject Classification:

34A12; 34N05

##  1\. Introduction

In 2019, the COVID-19 pandemic has appeared for the first time in Wuhan,
China, attracting many researchers to investigate outbreaks and the spread of
viruses [1, 2]. Some of the studies provided new mathematical compartmental
models, illustrating well the important contributions of Mathematics to fight
communicable diseases. Such models have been used to analyze the corresponding
dynamics and to supply useful techniques in disease epidemiology [3, 11, 13].

Kim et al. studied a SAIQH
(Susceptible–Asymptomatic–Infectious–Quarantined–Hospitalized) mathematical
model to analyze the transmission dynamics of MERS and to estimate
transmission rates [10]. Lemos-Paião et al. developed a SAIQH type model for
COVID-19, which was important to describe and understand the pandemic in
Portugal [11]. Similar to [11], here we also consider the $H_{IC}$ class of
hospitalized individuals in intensive care units.

Let the total living population under study, denoted by $N(t)$, $t\geq 0$, be
divided into six classes: (i) the susceptible individuals $S(t)$; (ii) the
infected individuals without (or with mild) symptoms $A(t)$ (the
Asymptomatic); (iii) infected individuals $I(t)$ with visible symptoms; (iv)
quarantined individuals $Q(t)$ in isolation at home; (v) hospitalized
individuals $H(t)$; (vi) and hospitalized individuals $H_{IC}(t)$ in intensive
care units. The mathematical model introduced and studied in [11] reads:

$\begin{cases}\dot{S}(t)=\Lambda+\omega nQ(t)-[\lambda(t)(1-p)+\phi
p+\gamma]S(t),\\\ \dot{A}(t)=\lambda(t)(1-p)S(t)-[q\nu+\gamma]A(t),\\\
\dot{I}(t)=q\nu A(t)-[\delta_{1}+\gamma]I(t),\\\ \dot{Q}(t)=\phi
pS(t)+\delta_{1}f_{1}I(t)+\delta_{2}(1-f_{2}-f_{3})H(t)-[\omega
m+\gamma]Q(t),\\\
\dot{H}(t)=\delta_{1}(1-f_{1})I(t)+\eta(1-k)H_{IC}(t)-\left[\delta_{2}(1-f_{2}-f_{3})+\delta_{2}f_{2}+\alpha_{1}f_{3}+\gamma\right]H(t),\\\
\dot{H_{IC}}(t)=\delta_{2}f_{2}H(t)-[\eta(1-k)+\alpha_{2}k+\gamma]H_{IC}(t),\end{cases}$
(1.1)

where, for all time $t\geq 0$,

$\lambda(t)=\dfrac{\beta\left(l_{A}A(t)+I(t)+l_{H}H(t)\right)}{N(t)}$

is a bounded positive function with

$N(t)=S(t)+A(t)+I(t)+Q(t)+H(t)+H_{IC}(t),$

$\beta,l_{A},l_{H}>0$, and all the other parameters in the model are
nonnegative. In addition,

$p,1-p,k,1-k,q,f_{1},1-f_{1},f_{2},f_{3},1-f_{2}-f_{3}\in[0,1].$

For more details on the mathematical model (1.1) we refer the reader to [11].

Nowadays, the study of dynamic equations on time scales is a leading topic of
research [4, 5, 6, 7, 9, 14, 15, 16, 17]. Motivated by the above mentioned
references and the 2019 paper of Khuddush and Prasad on a $n$-species
Lotka–Volterra system [9], here we extend (1.1) to a time scale $\mathbb{T}$,
studying the permanence and uniform asymptotic stability of the unique
positive almost periodic solution of the following SAIQH type model on time
scales:

$\begin{cases}S^{\Delta}(t)=\Lambda+\omega nQ(t)-[\lambda(t)(1-p)+\phi
p+\gamma]S^{\sigma}(t),\\\
A^{\Delta}(t)=\lambda(t)(1-p)S(t)-[q\nu+\gamma]A^{\sigma}(t),\\\
I^{\Delta}(t)=q\nu A(t)-[\delta_{1}+\gamma]I^{\sigma}(t),\\\
Q^{\Delta}(t)=\phi
pS(t)+\delta_{1}f_{1}I(t)+\delta_{2}(1-f_{2}-f_{3})H(t)-[\omega
m+\gamma]Q^{\sigma}(t),\\\
H^{\Delta}(t)=\delta_{1}(1-f_{1})I(t)+\eta(1-k)H_{IC}(t)-\left[\delta_{2}(1-f_{2}-f_{3})+\delta_{2}f_{2}+\alpha_{1}f_{3}+\gamma\right]H^{\sigma}(t),\\\
H_{IC}^{\Delta}(t)=\delta_{2}f_{2}H(t)-[\eta(1-k)+\alpha_{2}k+\gamma]H_{IC}^{\sigma}(t),\end{cases}$
(1.2)

where all parameters are real nonnegative numbers and $\lambda(t)$ is bounded
positive, as mentioned before. In the particular case
$\mathbb{T}=\mathbb{R}^{+}$, system (1.2) reduces to (1.1).

The text is organized as follows. In Section 2, we recall the main necessary
notions and results needed in the sequel. The original results are then
formulated and proved in Section 3: we show that system (1.2) is permanent
(Theorem 3); we give a sufficient condition for existence of a solution to
system (1.2) (Theorem 5); and we obtain conditions implying the dynamic system
(1.2) to have a unique almost periodic solution that is uniformly
asymptotically stable (Theorem 6).

For convenience, in the sequel we put $x_{1}=S$, $x_{2}=A$, $x_{3}=I$,
$x_{4}=Q$, $x_{5}=H$, and $x_{6}=H_{IC}$.

##  2\. Preliminaries

For a function $f(t)$ defined on $t\in\mathbb{T}^{+}$, we set

$f^{L}:=\inf\left\\{f(t):t\in\mathbb{T}^{+}\right\\},\;\;f^{U}:=\sup\left\\{f(t):\,t\in\mathbb{T}^{+}\right\\}.$

We begin by recalling some basic but fundamental concepts of the calculus on
time scales [4]. A time scale $\mathbb{T}$ is a nonempty closed subset of the
real numbers $\mathbb{R}$; on $\mathbb{T}$ one considers the topology
inherited from the real numbers with the standard topology. The jump operators
$\sigma,\rho:\,\mathbb{T}\longrightarrow\mathbb{T}$ and the graininess
function $\mu:\mathbb{T}\longrightarrow\mathbb{R}^{+}$ are defined,
respectively, by

$\sigma(t)=\inf\left\\{\tau\in\mathbb{T}:\tau>t\right\\},\;\rho(t)=\sup\left\\{\tau\in\mathbb{T}:\tau<t\right\\},\;\text{
and }\;\mu(t)=\sigma(t)-t,$

supplemented by $\inf\emptyset=\sup\mathbb{T}$ and
$\sup\emptyset=\inf\mathbb{T}$. A point $t\in\mathbb{T}$ is left-dense, left-
scattered, right-dense, or right-scattered when $\rho(t)=t$, $\rho(t)<t$,
$\sigma(t)=t$, and $\sigma(t)>t$, respectively. If $\mathbb{T}$ has a left-
scattered maximum $m$, then
$\mathbb{T}^{\kappa}:=\mathbb{T}\setminus\left\\{m\right\\}$; otherwise,
$\mathbb{T}^{\kappa}:=\mathbb{T}$. If $f:\mathbb{T}\rightarrow\mathbb{R}$,
then $f^{\sigma}:\mathbb{T}\rightarrow\mathbb{R}$ is given by
$f^{\sigma}(t)=f(\sigma(t))$ for all $t\in\mathbb{T}$.

######  Definition 1.

A function $f:\mathbb{T}\longrightarrow\mathbb{R}$ is called regressive
provided $1+\mu(t)f(t)\neq 0$ for all $t\in\mathbb{T}^{\kappa}$; a function
$g:\mathbb{T}\longrightarrow\mathbb{R}$ is called rd-continuous provided it is
continuous at right-dense points in $\mathbb{T}$ and its left-sided limits
exist (finite) at left-dense points in $\mathbb{T}$. The set of all regressive
and rd-continuous functions $f:\mathbb{T}\longrightarrow\mathbb{R}$ will be
denoted by $\mathcal{R}=\mathcal{R}(\mathbb{T},\mathbb{R})$. We define the set
$\mathcal{R}^{+}=\mathcal{R}^{+}(\mathbb{T},\mathbb{R})=\left\\{f\in\mathcal{R}:1+\mu(t)f(t)>0,\;\text{for}\,\text{all}\;t\in\mathbb{T}\right\\}$.

######  Lemma 1 (See [8]).

Assume that $\alpha>0$, $b>0$, and $-\alpha\in\mathcal{R}^{+}$. If

$y^{\Delta}(t)\geq(\leq)\,b-\alpha y^{\sigma}(t),\quad y(t)>0,\quad
t\in[t_{0},\infty)_{\mathbb{T}},$

then

$y(t)\geq(\leq)\,\dfrac{b}{\alpha}\left[1+\left(\dfrac{\alpha
y(t_{0})}{b}-1\right)e_{(-\alpha)}(t,t_{0})\right],\quad
t\in[t_{0},\infty)_{\mathbb{T}},$

where $e_{\cdot}(\cdot,\cdot)$ is the standard exponential function of the
time-scale calculus [4].

######  Definition 2 (See [12]).

A time scale $\mathbb{T}$ is called an almost periodic time scale if

$\Pi=\left\\{\tau\in\mathbb{R}:\,t+\tau\in\mathbb{T},\,\text{for}\,\text{all}\,\,t\in\mathbb{T}\right\\}\neq\left\\{0\right\\}.$

######  Definition 3 (See [12]).

Let $\mathbb{T}$ be an almost periodic time scale. A function $x\in
C(\mathbb{T},\mathbb{R}^{n})$ is called an almost periodic function if the
$\varepsilon$-translation set of $x$, that is,

$E\left\\{\varepsilon,x\right\\}=\left\\{\tau\in\Pi:\,\left|x(t+\tau)-x(t)\right|<\varepsilon\,\text{
for}\,\text{all}\,\,t\in\mathbb{T}\right\\},$

is a relatively dense set in $\mathbb{T}$ for all $\varepsilon>0$ and there
exists a constant $l(\varepsilon)>0$ such that each interval of length
$l(\varepsilon)$ contains a $\tau\in E\left\\{\varepsilon,x\right\\}$ for
which $\lvert x(t+\tau)-x(t)\lvert<\varepsilon$ for all $t\in\mathbb{T}$. The
value $\tau$ is known as the $\varepsilon$-translation number of $x$ and
$l(\varepsilon)$ is called the inclusion length of
$E\left\\{\varepsilon,x\right\\}$.

######  Definition 4 (See [12]).

Let $\mathbb{D}$ be an open set in $\mathbb{R}^{n}$ and let $\mathbb{T}$ be a
positive almost periodic time scale. A function $f\in
C(\mathbb{T}\times\mathbb{D},\mathbb{R}^{n})$ is called an almost periodic
function in $t\in\mathbb{T}$ uniformly for $x\in\mathbb{D}$ if the
$\varepsilon$-translation set of $f$,

$E\left\\{\varepsilon,f,\mathbb{S}\right\\}=\left\\{\tau\in\Pi:\mid
f(t+\tau)-f(t)\mid<\varepsilon,\,\text{ for }\,\text{ all
}\,\,(t,x)\in\mathbb{T}\times\mathbb{S}\right\\},$

is a relatively dense set in $\mathbb{T}$ for all $\varepsilon>0$ and, for
each compact subset $\mathbb{S}$ of $\mathbb{D}$, that is, for any given
$\varepsilon>0$ and each compact subset $\mathbb{S}$ of $\mathbb{D}$, there
exists a constant $l(\varepsilon,\mathbb{S})>0$ such that each interval of
length $l(\varepsilon,\mathbb{S})$ contains $\tau(\varepsilon,\mathbb{S})\in
E\left\\{\varepsilon,f,\mathbb{S}\right\\}$ for which

$\mid f(t+\tau,x)-f(t,x)\mid<\varepsilon,\,\,\text{ for }\,\text{ all
}\,\,(t,x)\in\mathbb{T}\times\mathbb{S}.$

Consider a system

$x^{\Delta}(t)=h(t,x),$ (2.1)

where
$h:\mathbb{T}^{+}\times\mathbb{S}_{B}\longrightarrow\mathbb{R}^{n},\;\mathbb{S}_{B}=\left\\{x\in\mathbb{R}^{n}:\,\|x\|<B\right\\}$
and $h(t,x)$ is almost periodic in $t$ uniformly for $x\in\mathbb{S}_{B}$ and
is continuous in $x$. In [17] the question of existence of a unique almost
periodic solution $f(t)\in\mathbb{S}$ of (2.1), which is uniformly
asymptotically stable, is investigated. For our model, we obtain from [17] the
following result.

######  Lemma 2 (cf. [17]).

Suppose that there exists a Lyapunov function $V(t,x,z)$ defined on
$\mathbb{T}^{+}\times\mathbb{S}_{B}\times\mathbb{S}_{B}$ satisfying the
following conditions:

(i) $\displaystyle a(\|x-z\|)\leq V(t,x,z)\leq b(\|x-z\|),\;\text{ where
}\;a,b\in\mathbb{K}$ $\displaystyle\text{ with }\mathbb{K}=\left\\{\alpha\in
C(\mathbb{R}^{+},\mathbb{R}^{+}):\alpha(0)=0\text{ and }\;\alpha\;\text{ is
increasing}\right\\};$ (ii) $\displaystyle|V(t,x,z)-V(t,x_{1},z_{1})|\leq
L(\|x-x_{1}\|+\|z-z_{1}\|),\;\text{ where }\;L>0\,\text{ is a constant};$
(iii) $\displaystyle D^{+}V^{\Delta}(t,x,z)\leq-cV(t,x,z),\;\text{ where
}\;c>0,\ -c\in\mathcal{R^{+}},$

and $D^{+}V^{\Delta}$ is the Dini derivative of $V$. Furthermore, if there
exists a solution $x(t)\in\mathbb{S}$ of system

$\begin{cases}x^{\Delta}_{1}(t)=\Lambda+\omega nx_{4}(t)-[\lambda(t)(1-p)+\phi
p+\gamma]x_{1}^{\sigma}(t),\\\
x^{\Delta}_{2}(t)=\lambda(t)(1-p)x_{1}(t)-[q\nu+\gamma]x_{2}^{\sigma}(t),\\\
x^{\Delta}_{3}(t)=q\nu x_{2}(t)-[\delta_{1}+\gamma]x_{3}^{\sigma}(t),\\\
x^{\Delta}_{4}(t)=\phi
px_{1}(t)+\delta_{1}f_{1}x_{3}(t)+\delta_{2}(1-f_{2}-f_{3})x_{5}(t)-[\omega
n+\gamma]x_{4}^{\sigma}(t),\\\
x^{\Delta}_{5}(t)=\delta_{1}(1-f_{1})x_{3}(t)+\eta(1-k)x_{6}(t)-[\delta_{2}(1-f_{2}-f_{3})+\delta_{2}f_{2}+\alpha_{1}f_{3}+\gamma]x_{5}^{\sigma}(t),\\\
x^{\Delta}_{6}(t)=\delta_{2}f_{2}x_{5}(t)-[\eta(1-k)+\alpha_{2}k+\gamma]x_{6}^{\sigma}(t),\end{cases}$
(2.2)

for $t\in\mathbb{T^{+}}$, where $\mathbb{S}\cup\mathbb{S}_{B}$ is a compact
set, then there exists a unique almost periodic solution $f(t)\in\mathbb{S}$
of system (2.2), which is uniformly asymptotically stable.

##  3\. Main Results

Let $t_{0}\in\mathbb{T}$ be a fixed positive initial time. Our main results
are: the proof that system (2.2) is permanent (Section 3.1); a sufficient
condition for existence of a solution to system (2.2) (Section 3.2); and
conditions that imply the dynamic system (2.2) to have a unique almost
periodic solution that is uniformly asymptotically stable (Section 3.3).

###  3.1. Permanence of solutions

We begin by introducing the notion of permanence of solutions.

######  Definition 5.

System (2.2) is said to be permanent if there exist positive constants $m$ and
$M$ such that
$m\leq\liminf_{t\longrightarrow\infty}x_{i}(t)\leq\limsup_{t\longrightarrow\infty}x_{i}(t)\leq
M$, $i=1,2,\ldots,6$, for any solution $(x_{1}(t),\ldots,x_{6}(t))$ of (2.2).

######  Theorem 3.

System (2.2) is permanent.

###### Proof.

Let $Z(t)=(x_{1}(t),\ldots,x_{6}(t))$ be a positive solution of system (2.2).
Then,

$\left\\{\begin{array}[]{l}x_{1}^{\Delta}(t)\geq\Lambda-[\lambda^{U}(1-p)+\phi
p+\gamma]x_{1}^{\sigma}(t),\\\
x_{2}^{\Delta}(t)\geq\lambda^{L}(1-p)m_{1}-(q\nu+\gamma)x_{2}^{\sigma}(t),\\\
x_{3}^{\Delta}(t)\geq q\nu m_{2}-(\delta_{1}+\gamma)x_{3}^{\sigma}(t),\\\
x_{4}^{\Delta}(t)\geq\phi pm_{1}+\delta_{1}f_{1}m_{3}-(\omega
n+\gamma)x_{4}^{\sigma}(t),\\\
x_{5}^{\Delta}(t)\geq\delta_{1}(1-f_{1})m_{3}-[\delta_{2}(1-f_{2}-f_{3})+\delta_{2}f_{2}+\alpha_{1}f_{3}+\gamma]x_{5}^{\sigma}(t),\\\
x_{6}^{\Delta}(t)\geq\delta_{2}f_{2}m_{5}-[\eta(1-k)+\alpha_{2}k+\gamma]x^{\sigma}_{6}(t).\end{array}\right.$

From Lemma 1, it follows that

$x_{1}(t)\geq\dfrac{\Lambda}{\lambda^{u}(1-p)+\phi
p+\gamma}\left[1+\left(\dfrac{\lambda^{U}(1-p)+\phi
p+\gamma}{\Lambda}x_{1}(t_{0})-1\right)e_{-(\lambda^{U}(1-p)+\phi
p+\gamma)}(t,t_{0})\right]$

and we have $e_{-(\lambda^{U}(1-p)+\phi p+\gamma)}(t,t_{0})\longrightarrow
0,\;\;as\;\;t\longrightarrow\infty$. Thus,

$\left\\{\begin{array}[]{ll}x_{1}(t)\geq
m_{1}:=\dfrac{\Lambda}{\lambda^{U}(1-p)+\phi p+\gamma},&\text{ for }\,t\geq
T_{1},\\\ x_{2}(t)\geq
m_{2}:=\dfrac{\lambda^{L}(1-p)m_{1}}{q\nu+\gamma},&\text{ for }\,t\geq
T_{2},\\\ x_{3}(t)\geq m_{3}:=\dfrac{q\nu m_{2}}{\delta_{1}+\gamma},&\text{
for }\,t\geq T_{3},\\\ x_{4}(t)\geq m_{4}:=\dfrac{\phi
pm_{1}+\delta_{1}f_{1}m_{3}}{\omega n+\gamma},&\text{ for }\,t\geq T_{4},\\\
x_{5}(t)\geq
m_{5}:=\dfrac{\delta_{1}(1-f_{1})m_{3}}{\delta_{2}(1-f_{2}-f_{3})+\delta_{2}f_{2}+\alpha_{1}f_{3}+\gamma},&\text{
for }\,t\geq T_{5},\\\ x_{6}(t)\geq
m_{6}:=\dfrac{\delta_{2}f_{2}m_{5}}{\eta(1-k)+\alpha_{2}k+\gamma},&\text{ for
}\,t\geq T_{6}.\end{array}\right.$

Let $m=\displaystyle\min_{1\leqslant i\leqslant 6}{m_{i}}$ and
$T=\displaystyle\max_{1\leqslant i\leqslant 6}{T_{i}}$. We can then write that
$x_{i}(t)\geq m$ for all $t>T$. ∎

###  3.2. Existence of solution

For system (2.2), we introduce the following assumption:

  * $(H_{1})$

$\lambda(t)$ is a bounded almost periodic function and satisfy
$0<\lambda^{L}\leq\lambda(t)\leq\lambda^{U}$.

To prove existence of solution, we first begin with a technical lemma.

######  Lemma 4.

If $(H1)$ holds, then, for any positive solution
$Z(t)=(x_{1}(t),\ldots,x_{6}(t))$ of system (2.2), there exist positive
constants $M$ and $T$ such that $x_{i}(t)<M$, $i=1,\ldots,6$, for all $t>T$.

###### Proof.

Let $Z(t)=(x_{1}(t),\ldots,x_{6}(t))$ be a positive solution of system (2.2).
Then,

$\left\\{\begin{array}[]{l}x_{1}^{\Delta}(t)\leqslant\Lambda+\omega
n\frac{\Lambda}{\gamma}-[\lambda^{L}(1-p)+\phi p+\gamma]x_{1}^{\sigma}(t),\\\
x_{2}^{\Delta}(t)\leqslant\lambda^{U}(1-p)M_{1}-(q\nu+\gamma)x_{2}^{\sigma}(t),\\\
x_{3}^{\Delta}(t)\leqslant q\nu M_{2}-(\delta_{1}+\gamma)x_{3}^{\sigma}(t),\\\
x_{4}^{\Delta}(t)\leqslant\phi
pM_{1}+\delta_{1}f_{1}M_{3}+\delta_{2}(1-f_{2}-f_{3})\frac{\Lambda}{\gamma}-(\omega
n+\gamma)x_{4}^{\sigma}(t),\\\
x_{5}^{\Delta}(t)\leqslant\delta_{1}(1-f_{1})M_{3}+\eta(1-k)\frac{\Lambda}{\gamma}-[\delta_{2}(1-f_{2}-f_{3})+\delta_{2}f_{2}+\alpha_{1}f_{3}+\gamma]x_{5}^{\sigma}(t),\\\
x_{6}^{\Delta}(t)\leqslant\delta_{2}f_{2}M_{5}-[\eta(1-k)+\alpha_{2}k+\gamma]x_{6}^{\sigma}(t).\end{array}\right.$

From Lemma 1,

$x_{1}(t)\leqslant\dfrac{\Lambda+\omega
n\frac{\Lambda}{\gamma}}{\lambda^{L}(1-p)+\phi
p+\gamma}\left[1+\left(\dfrac{\lambda^{L}(1-p)+\phi p+\gamma}{\Lambda+\omega
n\frac{\Lambda}{\gamma}}x_{1}(t_{0})-1\right)e_{-(\lambda^{L}(1-p)+\phi
p+\gamma)}(t,t_{0})\right]$

and we have $e_{-(\lambda^{L}(1-p)+\phi p+\gamma)}(t,t_{0})\longrightarrow
0,\;\;as\;\;t\longrightarrow\infty$. Then,

$x_{1}(t)\leqslant M_{1}:=\dfrac{\Lambda+\omega
n\frac{\Lambda}{\gamma}}{\lambda^{L}(1-p)+\phi p+\gamma},$
$\left\\{\begin{array}[]{ll}x_{1}(t)\leqslant M_{1}:=\dfrac{\Lambda+\omega
n\frac{\Lambda}{\gamma}}{\lambda^{L}(1-p)+\phi p+\gamma},&\text{ for }\,t\geq
T_{1},\\\ x_{2}(t)\leqslant
M_{2}:=\dfrac{\lambda^{U}(1-p)M_{1}}{q\nu+\gamma},&\text{ for }\,t\geq
T_{2},\\\ x_{3}(t)\leqslant M_{3}:=\dfrac{q\nu
M_{2}}{\delta_{1}+\gamma},&\text{ for }\,t\geq T_{3},\\\ x_{4}(t)\leqslant
M_{4}:=\dfrac{\phi
pM_{1}+\delta_{1}f_{1}M_{3}+\delta_{2}(1-f_{2}-f_{3})\frac{\Lambda}{\gamma}}{\omega
n+\gamma},&\text{ for }\,t\geq T_{4},\\\ x_{5}(t)\leqslant
M_{5}:=\dfrac{\delta_{1}(1-f_{1})M_{3}+\eta(1-k)\frac{\Lambda}{\gamma}}{\delta_{2}(1-f_{2}-f_{3})+\delta_{2}f_{2}+\alpha_{1}f_{3}+\gamma},&\text{
for }\,t\geq T_{5},\\\ x_{6}(t)\leqslant
M_{6}:=\dfrac{\delta_{2}f_{2}M_{5}}{\eta(1-k)+\alpha_{2}k+\gamma},&\text{ for
}\,t\geq T_{6}.\end{array}\right.$

Let $M=\displaystyle\max_{1\leqslant i\leqslant 6}{M_{i}}$ and
$T=\displaystyle\max_{1\leqslant i\leqslant 6}{T_{i}}$. Then
$x_{i}(t)\leqslant M$ for all $t>T$. ∎

Define

$\displaystyle\Omega$ $\displaystyle=$
$\displaystyle\\{(x_{1}(t),\ldots,x_{6}(t)):\,(x_{1}(t),\ldots,x_{6}(t))$
$\displaystyle\textit{ is a solution of (\ref{1.1}) and }\,0<m\leq x_{i}\leq
M,\;i=1,\ldots,6\\}.$

######  Theorem 5.

If $(H1)$ holds, then the set $\Omega$ is nonempty.

###### Proof.

The almost periodicity of $\lambda(t)$ implies that there is a sequence
$\\{\xi_{l}\\}\subseteq\mathbb{T}^{+}$ with $l\rightarrow\infty$ such that

$\lambda(t+\xi_{l})\rightarrow\lambda(t)\;\;\text{as}\;\;l\rightarrow\infty.$

From Theorem 3 and Lemma 4, for each sufficiently small $\varepsilon>0$, there
exists a $t_{1}\in\mathbb{T}^{+}$ such that

$m-\varepsilon\leq x_{i}(t)\leq M+\varepsilon,\,\text{ for all }\,t\geq
t_{1},\quad i=1,\ldots,6.$

Set $x_{il}(t)=x_{i}(t+\xi_{l})$ for $t\geq t_{1}-\xi_{l},\,l=1,2,\ldots$ For
any positive integer $k$, we obtain that there exists a sequence
$\\{x_{il}(t):\,l\geq k\\}$ such that the sequence $\\{x_{il}(t)\\}$ has a
subsequence, denoted by $\\{x_{il}^{*}(t)\\}$
$(x_{il}^{*}(t)=x_{il}(t+\xi_{l}^{*}))$, converging on any finite interval of
$\mathbb{T}^{+}$ as $l\rightarrow\infty$. So we have a sequence
$\\{y_{i}(t)\\}$ such that, for $t\in\mathbb{T}^{+}$,

$x_{il}^{*}(t)\longrightarrow y_{i}(t),\;\mbox{ as
}l\rightarrow\infty,\;i=1,\ldots,6.$ (3.1)

It is easy to see that the above sequence
$\\{\xi_{l}^{*}\\}\subseteq\mathbb{T}^{+}$ with $\xi_{l}^{*}\rightarrow\tau$
for $l\rightarrow\infty$ is such that

$\lambda(t+\xi_{l}^{*})\longrightarrow\lambda(t),\;\text{ as
}\;l\rightarrow\infty,$

which, together with (3.1) and

$\left\\{\begin{array}[]{l}x^{*\Delta}_{1l}(t)=\Lambda+\omega
nx^{*}_{4l}(t)-[\lambda(t+\xi_{l}^{*})(1-p)+\phi
p+\gamma]x_{1l}^{*\sigma}(t),\\\
x^{*\Delta}_{2l}(t)=\lambda(t+\xi_{l}^{*}))(1-p)x_{1l}^{*}(t)-[q\nu+\gamma]x_{2l}^{*\sigma}(t),\\\
x^{*\Delta}_{3l}(t)=q\nu
x_{2l}^{*}(t)-[\delta_{1}+\gamma]x_{3l}^{*\sigma}(t),\\\
x^{*\Delta}_{4l}(t)=\phi
px_{1l}^{*}(t)+\delta_{1}f_{1}x_{3l}^{*}(t)+\delta_{2}(1-f_{2}-f_{3})x_{5l}^{*}(t)-[\omega
n+\gamma]x^{*\sigma}_{4l}(t),\\\
x^{*\Delta}_{5l}(t)=\delta_{1}(1-f_{1})x_{3l}^{*}(t)+\eta(1-k)x_{6l}^{*}(t)-[\delta_{2}(1-f_{2}-f_{3})+\delta_{2}f_{2}+\alpha_{1}f_{3}+\gamma]x_{5l}^{*\sigma}(t),\\\
x^{*\Delta}_{6l}(t)=\delta_{2}f_{2}x_{5l}^{*}(t)-[\eta(1-k)+\alpha_{2}k+\gamma]x_{6l}^{*\sigma}(t),\end{array}\right.$

yields

$\left\\{\begin{array}[]{l}y^{\Delta}_{1}(t)=\Lambda+\omega
ny_{4}(t)-[\lambda(t)(1-p)+\phi p+\gamma]y_{1}^{\sigma}(t),\\\
y^{\Delta}_{2}(t)=\lambda(t)(1-p)y_{1}(t)-[q\nu+\gamma]y_{2}^{\sigma}(t),\\\
y^{\Delta}_{3}(t)=q\nu y_{2}(t)-[\delta_{1}+\gamma]y_{3}^{\sigma}(t),\\\
y^{\Delta}_{4}(t)=\phi
py_{1}(t)+\delta_{1}f_{1}y_{3}(t)+\delta_{2}(1-f_{2}-f_{3})y_{5}(t)-[\omega
n+\gamma]y^{\sigma}_{4}(t),\\\
y^{\Delta}_{5}(t)=\delta_{1}(1-f_{1})y_{3}(t)+\eta(1-k)y_{6}(t)-\left[\delta_{2}(1-f_{2}-f_{3})+\delta_{2}f_{2}+\alpha_{1}f_{3}+\gamma\right]y_{5}^{\sigma}(t),\\\
y^{\Delta}_{6}(t)=\delta_{2}f_{2}y_{5}(t)-\left[\eta(1-k)+\alpha_{2}k+\gamma\right]y_{6}^{\sigma}(t).\end{array}\right.$

It is clear that $(y_{1}(t),\ldots,y_{6}(t))$ is a solution of system (2.2)
and

$m-\varepsilon\leq y_{i}(t)\leq M+\varepsilon,\,\text{ for all
}\,t\in\mathbb{T}^{+},\quad i=1,\ldots,6.$

Since $\varepsilon$ is arbitrary, it follows that

$m\leq y_{i}(t)\leq M,\,\text{ for }\,t\in\mathbb{T}^{+},\quad i=1,\ldots,6.$

The proof is complete. ∎

###  3.3. Uniform asymptotic stability

Now, we establish sufficient conditions for the existence of a unique positive
almost periodic solution to system (2.2) that is uniform asymptotically
stable. We introduce some more notations. Let

$\displaystyle A_{1}$ $\displaystyle:=$ $\displaystyle\lambda^{L}(1-p)+\phi
p+\gamma;\quad A_{2}:=q\nu+\gamma;$ $\displaystyle A_{3}$ $\displaystyle:=$
$\displaystyle\delta_{1}+\gamma;\quad A_{4}:=\omega n+\gamma;$ $\displaystyle
A_{5}$ $\displaystyle:=$
$\displaystyle\delta_{2}(1-f_{3})+\alpha_{1}f_{3}+\gamma;\quad
A_{6}:=\eta(1-k)+\alpha_{2}k+\gamma;$ $\displaystyle B_{1}$ $\displaystyle:=$
$\displaystyle\lambda^{U}(1-p)+\phi p;\quad B_{2}:=q\nu+2\frac{\gamma\beta
l_{A}(1-p)M}{\Lambda};$ $\displaystyle B_{3}$ $\displaystyle:=$
$\displaystyle\delta_{1}+2\frac{\gamma\beta(1-p)M}{\Lambda};\quad
B_{4}:=\omega n;$ $\displaystyle B_{5}$ $\displaystyle:=$
$\displaystyle\delta_{2}(1-f_{3})+2\frac{\gamma\beta
l_{H}(1-p)M}{\Lambda};\quad B_{6}:=\eta(1-k).$

Moreover, let $A:=\displaystyle\min_{1\leq i\leq 6}A_{i}$ and
$B:=\displaystyle\max_{1\leq i\leq 6}B_{i}$.

In our next result (Theorem 6), we assume the following additional hypothesis:

  * (H2)

$B<A$.

######  Theorem 6.

If (H1) and (H2) hold, then the dynamic system (2.2) has a unique almost
periodic solution $Z(t)=(x_{1}(t),\ldots,x_{6}(t))\in\Omega$ that is uniformly
asymptotically stable.

###### Proof.

According to Theorem 3, every solution $Z(t)=(x_{1}(t),\ldots,x_{6}(t))$ of
system (2.2) satisfies $x_{i}^{L}\leq x_{i}\leq x_{i}^{U}$. Hence,
$\left|x_{i}(t)\right|\leq K_{i},\;\;i=1,\ldots,6$. Suppose that
$Z(t)=(x_{1}(t),\ldots,x_{6}(t))$ and
$\widehat{Z}(t)=(\widehat{x}_{1}(t),\ldots,\widehat{x}_{6}(t))$ are two
positive solutions of system (2.2). We have

$\left\\{\begin{array}[]{l}x^{\Delta}_{1}(t)=\Lambda+\omega
nx_{4}(t)-[\lambda(t)(1-p)+\phi p+\gamma]x_{1}^{\sigma}(t),\\\
x^{\Delta}_{2}(t)=\lambda(t)(1-p)x_{1}(t)-[q\nu+\gamma]x_{2}^{\sigma}(t),\\\
x^{\Delta}_{3}(t)=q\nu x_{2}(t)-[\delta_{1}+\gamma]x_{3}^{\sigma}(t),\\\
x^{\Delta}_{4}(t)=\phi
px_{1}(t)+\delta_{1}f_{1}x_{3}(t)+\delta_{2}(1-f_{2}-f_{3})x_{5}(t)-[\omega
n+\gamma]x^{\sigma}_{4}(t),\\\
x^{\Delta}_{5}(t)=\delta_{1}(1-f_{1})x_{3}(t)+\eta(1-k)x_{6}(t)-[\delta_{2}(1-f_{2}-f_{3})+\delta_{2}f_{2}+\alpha_{1}f_{3}+\gamma]x_{5}^{\sigma}(t),\\\
x^{\Delta}_{6}(t)=\delta_{2}f_{2}x_{5}(t)-[\eta(1-k)+\alpha_{2}k+\gamma]x_{6}^{\sigma}(t),\end{array}\right.$

and

$\left\\{\begin{array}[]{l}\widehat{x}^{\Delta}_{1}(t)=\Lambda+\omega
n\widehat{x}_{4}(t)-[\widehat{\lambda}(t)(1-p)+\phi
p+\gamma]\widehat{x}_{1}^{\sigma}(t),\\\
\widehat{x}^{\Delta}_{2}(t)=\widehat{\lambda}(t)(1-p)\widehat{x}_{1}(t)-[q\nu+\gamma]\widehat{x}_{2}^{\sigma}(t),\\\
\widehat{x}^{\Delta}_{3}(t)=q\nu\widehat{x}_{2}(t)-[\delta_{1}+\gamma]\widehat{x}^{\sigma}_{3}(t),\\\
\widehat{x}^{\Delta}_{4}(t)=\phi
p\widehat{x}_{1}(t)+\delta_{1}f_{1}\widehat{x}_{3}(t)+\delta_{2}(1-f_{2}-f_{3})\widehat{x}_{5}(t)-[\omega
n+\gamma]\widehat{x}_{4}^{\sigma}(t),\\\
\widehat{x}^{\Delta}_{5}(t)=\delta_{1}(1-f_{1})\widehat{x}_{3}(t)+\eta(1-k)\widehat{x}_{6}(t)-[\delta_{2}(1-f_{2}-f_{3})+\delta_{2}f_{2}+\alpha_{1}f_{3}+\gamma]\widehat{x}_{5}^{\sigma}(t),\\\
\widehat{x}^{\Delta}_{6}(t)=\delta_{2}f_{2}\widehat{x}_{5}(t)-[\eta(1-k)+\alpha_{2}k+\gamma]\widehat{x}_{6}^{\sigma}(t).\end{array}\right.$

Denote

$\left\|Z\right\|=\left\|(x_{1}(t),\ldots,x_{6}(t))\right\|=\displaystyle\sup_{t\in\mathbb{T}^{+}}\displaystyle\sum_{i=1}^{6}\left|x_{i}(t)\right|.$

Then $\left\|Z\right\|\leq K$ and $\left\|\widehat{Z}\right\|\leq K$ where
$K=\displaystyle\sum_{i=1}^{6}K_{i}$. Define the Lyapunov function
$V(t,Z,\widehat{Z})$ on $\mathbb{T}^{+}\times\Omega\times\Omega$ as

$V(t,Z,\widehat{Z})=\displaystyle\sum^{6}_{i=1}\left|x_{i}(t)-\widehat{x}_{i}(t)\right|=\displaystyle\sum^{6}_{i=1}V_{i}(t),$
(3.2)

where $V_{i}(t)=\left|x_{i}(t)-\widehat{x}_{i}(t)\right|$. Then the two norms

$\left\|Z-\widehat{Z}\right\|=\displaystyle\sum^{6}_{i=1}\left|x_{i}(t)-\widehat{x}_{i}(t)\right|$

and

$\left\|Z-\widehat{Z}\right\|_{*}=\displaystyle\sup_{t\in\mathbb{R}_{+}}\left[\displaystyle\sum^{6}_{i=1}\left(x_{i}(t)-\widehat{x}_{i}(t)\right)^{2}\right]^{\frac{1}{2}}$

are equivalent, that is, there exist two constants $\eta_{1},\,\eta_{2}>0$
such that

$\eta_{1}\left\|Z-\widehat{Z}\right\|_{*}\leq\left\|Z-\widehat{Z}\right\|\leq\eta_{2}\left\|Z-\widehat{Z}\right\|_{*}.$

Hence,

$\eta_{1}\left\|Z-\widehat{Z}\right\|_{*}\leq
V(t,Z,\widehat{Z})\leq\eta_{2}\left\|Z-\widehat{Z}\right\|_{*}.$

Let $a,\,b\in C(\mathbb{R}^{+},\mathbb{R}^{+})$, $a(x)=\eta_{1}x$, and
$b(x)=\eta_{2}x$. Then the assumption (i) of Lemma 2 is satisfied. On the
other hand, we have

$\displaystyle\left|V(t,Z(t),\widehat{Z}(t))-V(t,Z^{*}(t),\widehat{Z}^{*}(t))\right|$
$\displaystyle=$
$\displaystyle\left|\displaystyle\sum^{6}_{i=1}\left|x_{i}(t)-\widehat{x}_{i}(t)\right|-\displaystyle\sum^{6}_{i=1}\left|x_{i}^{*}(t)-\widehat{x}_{i}^{*}(t)\right|\right|$
$\displaystyle\leq$
$\displaystyle\displaystyle\sum^{6}_{i=1}\left|\left|x_{i}(t)-\widehat{x}_{i}(t)\right|-\left|x_{i}^{*}(t)-\widehat{x}_{i}^{*}(t)\right|\right|$
$\displaystyle\leq$
$\displaystyle\displaystyle\sum^{6}_{i=1}\left|x_{i}(t)-x^{*}_{i}(t)\right|+\displaystyle\sum^{6}_{i=1}\left|\widehat{x}_{i}(t)-\widehat{x}_{i}^{*}(t)\right|$
$\displaystyle\leq$
$\displaystyle\left\|Z-Z^{*}\right\|+\left\|\widehat{Z}-\widehat{Z}^{*}\right\|,$

where $L=1$, so condition (ii) of Lemma 2 is also satisfied.

Now, using Lemma 4.2 of [18], it follows that

$D^{+}V^{\Delta}_{i}(t)\leq
sign(x_{i}^{\sigma}(t)-\widehat{x}_{i}^{\sigma}(t))(x_{i}^{\Delta}(t)-\widehat{x}_{i}^{\Delta}(t)),\;\;i=1,\ldots,6,$

where $D^{+}V^{\Delta}_{i}$ is the Dini derivative of $V_{i}$. For $i=1$,

$\displaystyle D^{+}V^{\Delta}_{1}(t)$ $\displaystyle\leq$ $\displaystyle
sign(x_{1}^{\sigma}(t)-\widehat{x}_{1}^{\sigma}(t))(x_{1}^{\Delta}(t)-\widehat{x}_{1}^{\Delta}(t))$
$\displaystyle=$ $\displaystyle
sign(x_{1}^{\sigma}(t)-\widehat{x}_{1}^{\sigma}(t))[\omega
n(x_{4}(t)-\widehat{x}_{4}(t))-\lambda(t)(1-p)x_{1}^{\sigma}(t)$
$\displaystyle-[\phi
p+\gamma](x_{1}^{\sigma}(t)-\widehat{x}^{\sigma}_{1}(t))+\widehat{\lambda}(t)(1-p)\widehat{x}^{\sigma}_{1}(t)]$
$\displaystyle=$ $\displaystyle
sign(x_{1}^{\sigma}(t)-\widehat{x}_{1}^{\sigma}(t))[\omega
n(x_{4}(t)-\widehat{x}_{4}(t))-\lambda(t)(1-p)x_{1}^{\sigma}(t)$
$\displaystyle-[\phi
p+\gamma](x_{1}^{\sigma}(t)-\widehat{x}_{1}^{\sigma}(t))+\widehat{\lambda}(t)(1-p)\widehat{x}^{\sigma}_{1}(t)+\lambda(t)(1-p)\widehat{x}^{\sigma}_{1}(t)-\lambda(t)(1-p)\widehat{x}_{1}^{\sigma}(t)]$
$\displaystyle=$ $\displaystyle
sign(x_{1}^{\sigma}(t)-\widehat{x}_{1}^{\sigma}(t))[\omega
n(x_{4}(t)-\widehat{x}_{4}(t))$ $\displaystyle-[\lambda(t)(1-p)+\phi
p+\gamma](x_{1}^{\sigma}(t)-\widehat{x}_{1}^{\sigma}(t))$
$\displaystyle+(\widehat{\lambda}(t)-\lambda(t))(1-p)\widehat{x}_{1}^{\sigma}(t)]$
$\displaystyle=$ $\displaystyle
sign(x_{1}^{\sigma}(t)-\widehat{x}_{1}^{\sigma}(t))[\omega
n(x_{4}(t)-\widehat{x}_{4}(t))-[\lambda(t)(1-p)+\phi
p+\gamma](x_{1}^{\sigma}(t)-\widehat{x}_{1}^{\sigma}(t))$
$\displaystyle+(\widehat{\lambda}(t)-\lambda(t))(1-p)\widehat{x}^{\sigma}_{1}(t)]$
$\displaystyle\leq$ $\displaystyle\omega
n\left|x_{4}(t)-\widehat{x}_{4}(t)\right|-[\lambda^{L}(1-p)+\phi
p+\gamma]\left|x_{1}^{\sigma}(t)-\widehat{x}_{1}^{\sigma}(t)\right|$
$\displaystyle+\frac{\gamma\beta
l_{A}(1-p)M}{\Lambda}\left|x_{2}(t)-\widehat{x}_{2}(t)\right|+\frac{\gamma\beta(1-p)M}{\Lambda}\left|x_{3}(t)-\widehat{x}_{3}(t)\right|$
$\displaystyle+\frac{\gamma\beta
l_{H}(1-p)M}{\Lambda}\left|x_{5}(t)-\widehat{x}_{5}(t)\right|;$

for $i=2$,

$\displaystyle D^{+}V^{\Delta}_{2}(t)$ $\displaystyle\leq$ $\displaystyle
sign(x_{2}^{\sigma}(t)-\widehat{x}_{2}^{\sigma}(t))(x_{2}^{\Delta}(t)-\widehat{x}_{2}^{\Delta}(t))$
$\displaystyle=$ $\displaystyle
sign(x_{2}^{\sigma}(t)-\widehat{x}_{2}^{\sigma}(t))[\lambda(t)(1-p)x_{1}(t)-[q\nu+\gamma](x_{2}^{\sigma}(t)-\widehat{x}_{2}^{\sigma}(t))-\widehat{\lambda}(t)(1-p)\widehat{x}_{1}(t)]$
$\displaystyle=$ $\displaystyle
sign(x_{2}^{\sigma}(t)-\widehat{x}_{2}^{\sigma}(t))[\lambda(t)(1-p)x_{1}(t)-[q\nu+\gamma](x_{2}^{\sigma}(t)-\widehat{x}_{2}^{\sigma}(t))$
$\displaystyle-\widehat{\lambda}(t)(1-p)\widehat{x}_{1}(t)+\lambda(t)(1-p)\widehat{x}_{1}(t)-\lambda(t)(1-p)\widehat{x}_{1}(t)]$
$\displaystyle=$ $\displaystyle
sign(x_{2}^{\sigma}(t)-\widehat{x}_{2}^{\sigma}(t))[-[q\nu+\gamma](x_{2}^{\sigma}(t)-\widehat{x}_{2}^{\sigma}(t))$
$\displaystyle+\lambda(t)(1-p)(x_{1}(t)-\widehat{x}_{1}(t))-(\widehat{\lambda}(t)-\lambda(t))(1-p)\widehat{x}_{1}(t)]$
$\displaystyle\leq$
$\displaystyle\lambda^{U}(1-p)\left|x_{1}(t)-\widehat{x}_{1}(t)\right|-[q\nu+\gamma]\left|x_{2}^{\sigma}(t)-\widehat{x}_{2}^{\sigma}(t)\right|$
$\displaystyle+(1-p)M\left|\lambda(t)-\widehat{\lambda}(t)\right|$
$\displaystyle\leq$
$\displaystyle\lambda^{U}(1-p)\left|x_{1}(t)-\widehat{x}_{1}(t)\right|-[q\nu+\gamma]\left|x_{2}^{\sigma}(t)-\widehat{x}_{2}^{\sigma}(t)\right|$
$\displaystyle+\frac{\gamma\beta
l_{A}(1-p)M}{\Lambda}\left|x_{2}(t)-\widehat{x}_{2}(t)\right|$
$\displaystyle+\frac{\gamma\beta(1-p)M}{\Lambda}\left|x_{3}(t)-\widehat{x}_{3}(t)\right|+\frac{\gamma\beta
l_{H}(1-p)M}{\Lambda}\left|x_{5}(t)-\widehat{x}_{5}(t)\right|;$

for $i=3$,

$\displaystyle D^{+}V^{\Delta}_{3}(t)$ $\displaystyle\leq$ $\displaystyle
sign(x_{3}^{\sigma}(t)-\widehat{x}_{3}^{\sigma}(t))(x_{3}^{\Delta}(t)-\widehat{x}_{3}^{\Delta}(t))$
$\displaystyle=$ $\displaystyle
sign(x_{3}^{\sigma}(t)-\widehat{x}_{3}^{\sigma}(t))[q\nu(x_{2}(t)-\widehat{x}_{2}(t))$
$\displaystyle-[\delta_{1}+\gamma](x_{3}^{\sigma}(t)-\widehat{x}_{3}^{\sigma}(t))]$
$\displaystyle\leq$ $\displaystyle
q\nu\left|x_{2}(t)-\widehat{x}_{2}(t)\right|-[\delta_{1}+\gamma]\left|x_{3}^{\sigma}(t)-\widehat{x}_{3}^{\sigma}(t)\right|;$

for $i=4$,

$\displaystyle D^{+}V^{\Delta}_{4}(t)$ $\displaystyle\leq$ $\displaystyle
sign(x_{4}^{\sigma}(t)-\widehat{x}_{4}^{\sigma}(t))(x_{4}^{\Delta}(t)-\widehat{x}_{4}^{\Delta}(t))$
$\displaystyle=$ $\displaystyle
sign(x_{4}^{\sigma}(t)-\widehat{x}_{4}^{\sigma}(t))[\phi
p(x_{1}(t)-\widehat{x}_{1}(t))+\delta_{1}f_{1}(x_{3}(t)-\widehat{x}_{3}(t))$
$\displaystyle+\delta_{2}(1-f_{2}-f_{3})(x_{5}(t)-\widehat{x}_{5}(t))-[\omega
n+\gamma](x_{4}^{\sigma}(t)-\widehat{x}_{4}^{\sigma}(t))]$ $\displaystyle\leq$
$\displaystyle\phi
p\left|x_{1}(t)-\widehat{x}_{1}(t)\right|+\delta_{1}f_{1}\left|x_{3}(t)-\widehat{x}_{3}(t)\right|$
$\displaystyle+\delta_{2}(1-f_{2}-f_{3})\left|x_{5}(t)-\widehat{x}_{5}(t)\right|-[\omega
n+\gamma]\left|x_{4}^{\sigma}(t)-\widehat{x}_{4}^{\sigma}(t)\right|;$

for $i=5$,

$\displaystyle D^{+}V^{\Delta}_{5}(t)$ $\displaystyle\leq$ $\displaystyle
sign(x_{5}^{\sigma}(t)-\widehat{x}_{5}^{\sigma}(t))(x_{5}^{\Delta}(t)-\widehat{x}_{5}^{\Delta}(t))$
$\displaystyle=$ $\displaystyle
sign(x_{5}^{\sigma}(t)-\widehat{x}_{5}^{\sigma}(t))[\delta_{1}(1-f_{1})(x_{3}(t)-\widehat{x}_{3}(t))+\eta(1-k)(x_{6}(t)-\widehat{x}_{6}(t))$
$\displaystyle-[\delta_{2}(1-f_{2}-f_{3})+\delta_{1}f_{1}+\alpha_{1}f_{3}+\gamma](x_{5}^{\sigma}(t)-\widehat{x}_{5}^{\sigma}(t))]$
$\displaystyle\leq$
$\displaystyle\delta_{1}(1-f_{1})\left|x_{3}(t)-\widehat{x}_{3}(t)\right|+\eta(1-k)\left|x_{6}(t)-\widehat{x}_{6}(t)\right|$
$\displaystyle-[\delta_{2}(1-f_{2}-f_{3})+\delta_{1}f_{1}+\alpha_{1}f_{3}+\gamma]\left|x_{5}^{\sigma}(t)-\widehat{x}_{5}^{\sigma}(t)\right|;$

and, finally, for $i=6$,

$\displaystyle D^{+}V^{\Delta}_{6}(t)$ $\displaystyle\leq$ $\displaystyle
sign(x_{6}^{\sigma}(t)-\widehat{x}_{6}^{\sigma}(t))(x_{6}^{\Delta}(t)-\widehat{x}_{6}^{\Delta}(t))$
$\displaystyle=$ $\displaystyle
sign(x_{6}^{\sigma}(t)-\widehat{x}_{6}^{\sigma}(t))[\delta_{2}f_{2}(x_{5}(t)-\widehat{x}_{5}(t))$
$\displaystyle-[\eta(1-k)+\alpha_{2}k+\gamma](x_{6}^{\sigma}(t)-\widehat{x}_{6}^{\sigma}(t)$
$\displaystyle\leq$
$\displaystyle\delta_{2}f_{2}\left|x_{5}(t)-\widehat{x}_{5}(t)\right|-[\eta(1-k)+\alpha_{2}k+\gamma]\left|x_{6}^{\sigma}(t)-\widehat{x}_{6}^{\sigma}(t)\right|.$

It follows that

$\displaystyle D^{+}V^{\Delta}(t)$ $\displaystyle\leq$ $\displaystyle\omega
n\left|x_{4}(t)-\widehat{x}_{4}(t)\right|-[\lambda^{L}(1-p)+\phi
p+\gamma]\left|x_{1}^{\sigma}(t)-\widehat{x}_{1}^{\sigma}(t)\right|$
$\displaystyle+\frac{\gamma\beta
l_{A}(1-p)M}{\Lambda}\left|x_{2}(t)-\widehat{x}_{2}(t)\right|+\frac{\gamma\beta(1-p)M}{\Lambda}\left|x_{3}(t)-\widehat{x}_{3}(t)\right|$
$\displaystyle+\frac{\gamma\beta
l_{H}(1-p)M}{\Lambda}\left|x_{5}(t)-\widehat{x}_{5}(t)\right|+\lambda^{U}(1-p)\left|x_{1}(t)-\widehat{x}_{1}(t)\right|$
$\displaystyle-[q\nu+\gamma]\left|x_{2}^{\sigma}(t)-\widehat{x}_{2}^{\sigma}(t)\right|+\frac{\gamma\beta
l_{A}(1-p)M}{\Lambda}\left|x_{2}(t)-\widehat{x}_{2}(t)\right|$
$\displaystyle+\frac{\gamma\beta(1-p)M}{\Lambda}\left|x_{3}(t)-\widehat{x}_{3}(t)\right|+\frac{\gamma\beta
l_{H}(1-p)M}{\Lambda}\left|x_{5}(t)-\widehat{x}_{5}(t)\right|$
$\displaystyle+q\nu\left|x_{2}(t)-\widehat{x}_{2}(t)\right|-[\delta_{1}+\gamma]\left|x_{3}^{\sigma}(t)-\widehat{x}_{3}^{\sigma}(t)\right|$
$\displaystyle+\phi
p\left|x_{1}(t)-\widehat{x}_{1}(t)\right|+\delta_{1}f_{1}\left|x_{3}(t)-\widehat{x}_{3}(t)\right|$
$\displaystyle+\delta_{2}(1-f_{2}-f_{3})\left|x_{5}(t)-\widehat{x}_{5}(t)\right|-[\omega
n+\gamma]\left|x_{4}^{\sigma}(t)-\widehat{x}_{4}^{\sigma}(t)\right|$
$\displaystyle+\delta_{1}(1-f_{1})\left|x_{3}(t)-\widehat{x}_{3}(t)\right|+\eta(1-k)\left|x_{6}(t)-\widehat{x}_{6}(t)\right|$
$\displaystyle-[\delta_{2}(1-f_{2}-f_{3})+\delta_{1}f_{1}+\alpha_{1}f_{3}+\gamma]\left|x_{5}^{\sigma}(t)-\widehat{x}_{5}^{\sigma}(t)\right|$
$\displaystyle+\delta_{2}f_{2}\left|x_{5}(t)-\widehat{x}_{5}(t)\right|-[\eta(1-k)+\alpha_{2}k+\gamma]\left|x_{6}^{\sigma}(t)-\widehat{x}_{6}^{\sigma}(t)\right|$
$\displaystyle=$ $\displaystyle-[\lambda(t)(1-p)+\phi
p+\gamma]\left|x_{1}^{\sigma}(t)-\widehat{x}_{1}^{\sigma}(t)\right|-[q\nu+\gamma]\left|x_{2}^{\sigma}(t)-\widehat{x}_{2}^{\sigma}(t)\right|$
$\displaystyle-[\delta_{1}+\gamma]\left|x_{3}^{\sigma}(t)-\widehat{x}_{3}^{\sigma}(t)\right|-[\omega
n+\gamma]\left|x_{4}^{\sigma}(t)-\widehat{x}_{4}^{\sigma}(t)\right|$
$\displaystyle-[\delta_{2}(1-f_{2}-f_{3})+\delta_{1}f_{1}+\alpha_{1}f_{3}+\gamma]\left|x_{5}^{\sigma}(t)-\widehat{x}_{5}^{\sigma}(t)\right|$
$\displaystyle-[\eta(1-k)+\alpha_{2}k+\gamma]\left|x_{6}^{\sigma}(t)-\widehat{x}_{6}^{\sigma}(t)\right|$
$\displaystyle+\left\\{\lambda^{U}(1-p)+\phi
p\right\\}\left|x_{1}(t)-\widehat{x}_{1}(t)\right|$
$\displaystyle+\left\\{q\nu+\frac{\gamma\beta
l_{A}(1-p)M}{\Lambda}+\frac{\gamma\beta
l_{A}(1-p)M}{\Lambda}\right\\}\left|x_{2}(t)-\widehat{x}_{2}(t)\right|$
$\displaystyle+\left\\{\delta_{1}f_{1}+\delta_{1}(1-f_{1})+\frac{\gamma\beta(1-p)M}{\Lambda}+\frac{\gamma\beta(1-p)M}{\Lambda}\right\\}\left|x_{3}(t)-\widehat{x}_{3}(t)\right|$
$\displaystyle+\omega
n\left|x_{4}(t)-\widehat{x}_{4}(t)\right|+\left\\{\delta_{2}(1-f_{2}-f_{3})+\delta_{2}f_{2}+\frac{\gamma\beta
l_{H}(1-p)M}{\Lambda}\right.$ $\displaystyle\left.+\frac{\gamma\beta
l_{H}(1-p)M}{\Lambda}\right\\}\left|x_{5}(t)-\widehat{x}_{5}(t)\right|+\eta(1-k)\left|x_{6}(t)-\widehat{x}_{6}(t)\right|$
$\displaystyle=$ $\displaystyle-
A_{1}\left|x_{1}^{\sigma}(t)-\widehat{x}_{1}^{\sigma}(t)\right|-A_{2}\left|x_{2}^{\sigma}(t)-\widehat{x}_{2}^{\sigma}(t)\right|$
$\displaystyle-
A_{3}\left|x_{3}^{\sigma}(t)-\widehat{x}_{3}^{\sigma}(t)\right|-A_{4}\left|x_{4}^{\sigma}(t)-\widehat{x}_{4}^{\sigma}(t)\right|$
$\displaystyle-
A_{5}\left|x_{5}^{\sigma}(t)-\widehat{x}_{5}^{\sigma}(t)\right|-A_{6}\left|x_{6}^{\sigma}(t)-\widehat{x}_{6}^{\sigma}(t)\right|$
$\displaystyle+B_{1}\left|x_{1}(t)-\widehat{x}_{1}(t)\right|+B_{2}\left|x_{2}(t)-\widehat{x}_{2}(t)\right|$
$\displaystyle+B_{3}\left|x_{3}(t)-\widehat{x}_{3}(t)\right|+B_{4}\left|x_{4}(t)-\widehat{x}_{4}(t)\right|$
$\displaystyle+B_{5}\left|x_{5}(t)-\widehat{x}_{5}(t)\right|+B_{6}\left|x_{6}(t)-\widehat{x}_{6}(t)\right|$
$\displaystyle=$ $\displaystyle-AV(\sigma(t))+BV(t)$ $\displaystyle=$
$\displaystyle(B-A)V(t)-A\mu(t)D^{+}V^{\Delta}(t)$

and $D^{+}V^{\Delta}(t)\leq\frac{B-A}{1+A\mu(t)}V(t)\leq-\psi(t)V(t)$ with
$\psi=\frac{A-B}{1+A\mu^{U}}$. By $(H2)$, we have
$\psi(t)=\frac{A-B}{1+A\mu^{U}}>0$ and
$1-\psi\mu(t)=1+A(\mu^{U}-\mu(t))+\mu(t)B>0$. Hence,
$-\psi\in\mathcal{R}^{+}$. Thus, the assumption (iii) of Lemma 2 is satisfied
and it follows from Lemma 2 that there exists a unique almost periodic
solution $Z(t)=(x_{1}(t),\ldots,x_{6}(t))$ of the dynamic system (2.2) that is
uniformly asymptotically stable with $Z(t)\in\Omega$. ∎

We illustrate our results with an example.

######  Example 7.

Based on [11], let us consider the following system on the time scale
$\mathbb{T}=\mathbb{Z}_{0}^{+}$:

$\left\\{\begin{array}[]{l}x^{\Delta}_{1}(t)=\Lambda+\omega
nx_{4}(t)-[\lambda(t)(1-p)+\phi p+\gamma]x_{1}^{\sigma}(t),\\\
x^{\Delta}_{2}(t)=\lambda(t)(1-p)x_{1}(t)-[q\nu+\gamma]x_{2}^{\sigma}(t),\\\
x^{\Delta}_{3}(t)=q\nu x_{2}(t)-[\delta_{1}+\gamma]x_{3}^{\sigma}(t),\\\
x^{\Delta}_{4}(t)=\phi
px_{1}(t)+\delta_{1}f_{1}x_{3}(t)+\delta_{2}(1-f_{2}-f_{3})x_{5}(t)-[\omega
n+\gamma]x^{\sigma}_{4}(t),\\\
x^{\Delta}_{5}(t)=\delta_{1}(1-f_{1})x_{3}(t)+\eta(1-k)x_{6}(t)-[\delta_{2}(1-f_{2}-f_{3})+\delta_{2}f_{2}+\alpha_{1}f_{3}+\gamma]x_{5}^{\sigma}(t),\\\
x^{\Delta}_{6}(t)=\delta_{2}f_{2}x_{5}(t)-[\eta(1-k)+\alpha_{2}k+\gamma]x_{6}^{\sigma}(t),\end{array}\right.$
(3.3)

subject to

$x_{1}(0)=10283785,\quad x_{2}(0)=13,\quad x_{3}(0)=2,\quad x_{4}(0)=0,\quad
x_{5}(0)=0,\quad x_{6}(0)=0,$

where $\Lambda=\frac{22614}{53}$, $\omega=1/31$, $n=0.075$,

$\lambda(t)=\dfrac{\beta\left(l_{A}x_{2}(t)+x_{3}(t)+l_{H}x_{5}(t)\right)}{N(t)}$

with $\beta=1.93$, $l_{A}=1$, $l_{H}=0.1$, and $N(t)=\sum_{i=1}^{6}x_{i}(t)$,
$p=0.68$, $\phi=1/12$, $\gamma=\frac{47833615}{N_{0}}$ with $N_{0}=N(0)$,
$q=0.15$, $\nu=1/15$, $\delta_{1}=1/3$, $\delta_{2}=1/3$, $f_{1}=0.96$,
$f_{2}=0.21$, $f_{3}=0.03$, $\eta=1/7$, $k=0.03$, $\alpha_{1}=1/7$, and
$\alpha_{2}=1/15$.

System (3.3) is permanent with $\lambda^{L}=1.876738171\times 10^{-7}$,
$\lambda^{U}=1.93$, $M_{1}=65.83271997$, $M_{2}=8.722416333$,
$M_{3}=0.017498412509$, $M_{4}=4.428264471$, $M_{6}=0.02788991356$,
$M=\max_{i=1,\ldots,6}(M_{i})=M_{1}$, $m_{1}=58.16800031$, $m_{2}=7.7068897$,
$m_{3}=0.0154611$, $m_{4}=0.7093454$, $m_{5}=0.0000414$,
$m_{6}=6.0482208\times 10^{-7}$, and $m=\min_{i=1,\ldots,6}(m_{i})=m_{6}$. In
addition, the conditions of Theorem 6 are verified and we have

$4.653775371=A>B=4.148857053,\quad\psi=0.0505839,\quad
1-\psi\mu(t)=0.94941603>0.$

We conclude that system (3.3) has a unique positive almost periodic solution,
which is uniformly asymptotic stable. This is illustrated in Figure 1.

(a) Evolution of $x_{1}(t)$

(b) Evolution of $x_{2}(t)$

(c) Evolution of $x_{3}(t)$

(d) Evolution of $x_{4}(t)$

(e) Evolution of $x_{5}(t)$

(f) Evolution of $x_{6}(t)$

Figure 1. Example 7: solution of (3.3) during 7 days.

## Acknowledgments

This work is part of first authors’ PhD project. Torres is funded by the
Fundação para a Ciência e a Tecnologia, I.P. (FCT, Funder ID = 50110000187)
under grants UIDB/04106/2020 and UIDP/04106/2020 and the project
2022.03091.PTDC “Mathematical Modelling of Multi-scale Control Systems:
applications to human diseases” (CoSysM3), financially supported by national
funds (OE) through FCT/MCTES.

## References

  * [1] P. Agarwal, J. J. Nieto, M. Ruzhansky and D. F. M. Torres, Analysis of Infectious Disease Problems (Covid-19) and Their Global Impact, Infosys Science Foundation Series in Mathematical Sciences, Springer, Singapore, 2021.
  * [2] P. Agarwal, J. J. Nieto and D. F. M. Torres, Mathematical Analysis of Infectious Diseases, Academic Press, London, UK, 2022.
  * [3] D. I. Baleanu, H. Mohammadi and S. Rezapour, A fractional differential equation model for the COVID-19 transmission by using the Caputo-Fabrizio derivative, Adv. Difference Equ. 2020, Paper No. 299, 27 pp.
  * [4] M. Bohner and A. Peterson, Dynamic equations on time scales, Birkhäuser Boston, Boston, MA, 2001.
  * [5] M. Bohner and A. Peterson, Advances in dynamic equations on time scales, Birkhäuser Boston, Inc., Boston, MA, 2003.
  * [6] M. J. Bohner, S. H. Streipert and D. F. M. Torres, Exact solution to a dynamic SIR model, Nonlinear Anal. Hybrid Syst. 32 (2019), 228–238.
  * [7] M. J. Bohner, J. Heim and A. L. Liu, Solow models on time scales, Cubo 15 (2013), no. 1, 13–31.
  * [8] M. Hu and L. Wang, Dynamic inequalities on time scales with applications in permanence of predator-prey system, Discrete Dyn. Nat. Soc. 2012, Art. ID 281052, 15 pp.
  * [9] M. Khuddush and K. R. Prasad, Permanence and stability of multi-species nonautonomous Lotka-Volterra competitive systems with delays and feedback controls on time scales, Khayyam J. Math. 7 (2021), no. 2, 241–256.
  * [10] Y. Kim, S. Lee, C. Chu, S. Choe, S. Hong and Y. Shin, The characteristics of middle eastern respiratory syndrome coronavirus transmission dynamics in South Korea, Osong Public Health Res. Perspect. 7 (2016), no. 1, 49–55.
  * [11] A. P. Lemos-Paião, C. J. Silva and D. F. M. Torres, A new compartmental epidemiological model for COVID-19 with a case study of Portugal, Ecological Complexity 44 (2020), Art. 100885, 8 pp.
  * [12] W.-T. Li and L. L. Wang, Existence and global attractivity of positive periodic solutions of functional differential equations with feedback control, J. Comput. Appl. Math. 180 (2005), no. 2, 293–309.
  * [13] K. Shah, M. Arfan, W. Deebani, M. Shutaywi and D. Baleanu, Investigation of COVID-19 mathematical model under fractional order derivative, Math. Model. Nat. Phenom. 16 (2021), Paper No. 50, 14 pp.
  * [14] K. R. Prasad and M. Khuddush, Stability of positive almost periodic solutions for a fishing model with multiple time varying variable delays on time scales, Bull. Int. Math. Virtual Inst. 9 (2019), no. 3, 521–533.
  * [15] K. R. Prasad and M. Khuddush, Existence and uniform asymptotic stability of positive almost periodic solutions for three-species Lotka-Volterra competitive system on time scales, Asian-Eur. J. Math. 13 (2020), no. 3, 2050058, 24 pp.
  * [16] C. C. Tisdell and A. H. Zaidi, Basic qualitative and quantitative results for solutions to nonlinear, dynamic equations on time scales with an application to economic modelling, Nonlinear Anal. 68 (2008), no. 11, 3504–3524.
  * [17] H. T. Zhang and Y. Li, Almost periodic solutions to dynamic equations on time scales, J. Egyptian Math. Soc. 21 (2013), no. 1, 3–10.
  * [18] Y. Li and C. Wang, Almost periodic solutions to dynamic equations on time scales and applications, J. Appl. Math. 2012, Art. ID 463913, 19 pp.