# Bibliography management: natbib package

Luke S. Baker  School of Mathematical and Statistical Sciences, Arizona State
University, Tempe, AZ  Applied Mathematics and Plasma Physics, Los Alamos
National Laboratory, Los Alamos, NM<EMAIL_ADDRESS>Saif R.
Kazi  Applied Mathematics and Plasma Physics, Los Alamos National Laboratory,
Los Alamos, NM  Center for Nonlinear Studies, Los Alamos National Laboratory,
Los Alamos, NM  Anatoly Zlotnik  Applied Mathematics and Plasma Physics, Los
Alamos National Laboratory, Los Alamos, NM

(August 27, 2024)

#  Transitions from Monotonicity to Chaos in Gas Mixture Dynamics
in Pipeline Networks

Luke S. Baker  School of Mathematical and Statistical Sciences, Arizona State
University, Tempe, AZ  Applied Mathematics and Plasma Physics, Los Alamos
National Laboratory, Los Alamos, NM<EMAIL_ADDRESS>Saif R.
Kazi  Applied Mathematics and Plasma Physics, Los Alamos National Laboratory,
Los Alamos, NM  Center for Nonlinear Studies, Los Alamos National Laboratory,
Los Alamos, NM  Anatoly Zlotnik  Applied Mathematics and Plasma Physics, Los
Alamos National Laboratory, Los Alamos, NM

(August 27, 2024)

###### Abstract

The blending of hydrogen generated using clean energy into natural gas
pipeline networks is proposed in order to utilize existing energy systems for
their planned lifetimes while reducing their reliance on fossil fuels. We
formulate a system of partial differential equations (PDEs) that govern the
flow dynamics of mixtures of gases in pipeline networks under the influence of
time-varying compressor and regulator control actions. The formulation is
derived for general gas networks that can inject or withdraw arbitrary time-
varying mixtures of gases into or from the network at arbitrarily specified
nodes. The PDE formulation is discretized in space to form a nonlinear control
system that is used to prove that homogeneous mixtures are well-behaved and
heterogeneous mixtures may be ill-behaved in the sense of monotone-ordering of
solutions. We use numerical simulations to compute interfaces in the parameter
region of sinusoidal boundary conditions that delimit monotonic, periodic, and
chaotic system responses. The interfaces suggest that any solution in the
monotonic response region is not chaotic and will eventually approach a
periodic orbit. The results are demonstrated using examples for a single
pipeline and a small test network.

monotonicity; periodicity; network flow; phase transition

##  I Introduction

Although natural gas is projected to be a primary fuel source through the year
2050 [1], societies worldwide are investing intensively to transition from
fossil fuels such as natural gas and coal to more sustainable and cleaner
resources. Hydrogen is an energy carrier that can be cleanly produced [2] and
can address climate change because it does not result in carbon dioxide
emissions or other harmful emissions when it is burned. Several qualities of
hydrogen make it an attractive fuel option for a variety of applications that
include transportation and high temperature manufacturing. Hydrogen can also
be used to power turbines, which can potentially be used for aviation and
electric power production. Hydrogen can be produced directly from fossil
fuels, biomass, or direct electrolysis, by splitting water into its
constituent components of hydrogen and oxygen. After hydrogen is produced, it
can be transported to end users economically by dedicated pipeline systems.

Recent studies have proposed that natural gas pipelines can safely transport
mixtures of up to 20% hydrogen or more by volume [3, 4]. Hydrogen could be
transported through the existing infrastructure and then separated, or the
mixture could be burned directly as an end-use fuel. Because the physical and
chemical properties of hydrogen and natural gas (primarily methane) differ
significantly, the mass and energy transport dynamics of inhomogeneous
mixtures of these constituent gases are considerably more complex than for a
homogeneous gas [5]. The mathematical modeling of such mixtures is also
considerably more challenging than what has traditionally been done for gas
pipelines [6]. The introduction of substantial proportions of much lighter
hydrogen into natural gas pipelines requires much closer spacing of gas
compressors, and this relationship has been characterized in an empirical
study [7]. Additionally, the pressure and flow dynamics in gas networks have
been proven to satisfy certain physically intuitive and conceptually valuable
monotonicity properties [8], which must be re-examined in the presence of
inhomogeneous gas mixing.

The physical complexities of blending hydrogen in natural gas pipelines
present several mathematical challenges. Additional state variables are needed
to account for changes in mass fraction, which affect total density, energy
content, and flow dynamics. Modeling the flow of a homogeneous gas on a
network requires partial differential equations (PDEs) for mass and momentum
conservation on each pipe, and a linear mass flow balance equation at each
network junction. Injection of a second gas into the network requires the
addition of another PDE on each pipe and a bilinear nodal balance equation at
each junction to account for conservation of composition. This more than
doubles the state space of the continuous model. Moreover, the faster wave
speed corresponding to the lower density of hydrogen worsens the numerical
ill-conditioning of the dynamic model. Such issues have been highlighted by
the numerical simulations of hydrogen and natural gas flows in pipelines [9,
10, 11, 12, 13, 14, 15, 7, 16].

One recent study has demonstrated conditions under which pipeline pressures
may exceed allowable upper limits, and that the likelihood of this occurrence
increases proportionally with increasing hydrogen concentration [13]. Another
study examined the effects of hydrogen blending on the detection and
estimation of leaks [16], and demonstrated that the amount of leak discharge
increases as the concentration of hydrogen increases. A moving grid method and
an implicit backward difference method for tracking gas concentration were
both shown to perform well for numerical simulations but the implicit
difference method may lose some finer detail due to numerical diffusion [10].
The method of characteristics was also applied for the numerical simulation of
transient flows on cyclic networks with homogeneous flow mixtures [12].
Modeling networks of pipelines with composition tracking was the focus of
another recent study [7], although this model does not include control actions
of compressor units. In general, these models demonstrate a simulation
capability or sensitivity study for a specific network. Addressing challenging
design, operational, and economic issues in the pipeline transport of gas
mixtures will require minimal and generalizable mathematical models that
adequately describe the relevant physics, in addition to complex simulations
with comprehensive characterizations of pipeline flows.

The scope of the present study is threefold. First, we extend general control
system models for gas pipeline networks [17] to account for heterogeneous
mixtures of hydrogen and natural gas. Similar lumped parameter modeling has
been well-studied for pipeline simulation since at least a decade ago [18,
19]. The state variables are flows, partial densities, and pressures
throughout the network, and the control variables are the actions of
compressor and regulator units. Control actions may be designed to minimize
fuel consumption [20, 21, 22] or maximize economic value [23]. The PDE control
system of the mixture is discretized in space using a finite volume method
[24] and written in matrix form as a finite-dimensional control system of
nonlinear ordinary differential equations (ODEs). Second, we prove that
solutions to initial boundary value problems (IBVPs) of a gas mixture have
certain monotone ordering properties if the concentration is homogeneous but,
in general, do not have these properties if the concentration is
heterogeneous. The homogeneous monotonicity result generalizes the pure
natural gas monotonicity result for obtaining control formulations that are
robust to uncertainty in pressure and withdrawal profiles [25, 8]. Third, we
demonstrate that the solution of an IBVP may be irregular, in the sense of
generating a continuous distribution of harmonic modes, and may also be
chaotic, in the sense of being sensitive to initial conditions. Numerical
simulations are used to characterize flow solution behavior in a phase space
of periodic forcing functions and to identify boundaries between the regions
of monotonic, periodic, and chaotic solution behavior. Transitions through
such fluid mixing phase regions were observed in oceanic wind bursts [26, 27]
and in flame combustion of hydrogen and air mixtures [28, 29]. Inspection of
the response interfaces that we compute suggest that any solution in the
monotonic response region is not chaotic and will eventually approach a
periodic orbit. Simulation-based analyses such as those presented here could
be used to evaluate appropriate limitations on blending of hydrogen into
existing natural gas pipeline networks.

The rest of this paper is organized as follows. The PDEs that govern
heterogeneous mixtures of hydrogen and natural gas are presented in Section
II. In Section III, the PDE system is discretized in space to obtain a system
of ODEs. Section IV presents a derivation of equivalent ODE systems in terms
of other state variables of interest. Section V contains a proof that each of
the equivalent systems have monotonic solutions if the concentration is
homogeneous, as well as a proof that the solutions are, in general, non-
monotonic if the concentration is heterogeneous. In Section VI, we illustrate
non-monotonic system responses using numerical simulations of flows through a
small test network that contains a loop, and which was examined in a previous
study [30]. Moreover, that section illustrates that certain types of
equivalent systems may have more desirable monotone system behavior than
others in certain response regimes. Sections VII, VIII, and IX describe
techniques to compute interfaces, in the region of boundary condition
parameters for flow in a single pipeline, between regions that do and do not
exhibit monotonic, periodic, and chaotic properties, respectively. We provide
concluding remarks and an outlook for future work in Section X.

##  II Gas Network Modeling

A gas transport network is modeled as a connected and directed graph
$(\mathcal{E},\mathcal{V})$ consisting of edges $\mathcal{E}=\\{1,\dots,E\\}$
and nodes $\mathcal{V}=\\{1,\dots,V\\}$, where $E$ and $V$ denote the numbers
of edges and nodes, respectively. It is assumed that the elements of these
sets are ordered according to their integer labels. The edges represent
pipelines and the nodes represent junctions or stations where gas can be
injected into or withdrawn from the network. The symbol $k$ is reserved for
indexing edges in $\mathcal{E}$ and the symbols $i$ and $j$ are reserved for
indexing nodes in $\mathcal{V}$. The graph is directed by assigning a positive
flow direction along each edge. It is assumed that gas physically flows in
only the direction of positive flow, so that the mass flow and velocity values
of the gas are positive quantities everywhere in the network. The notation
$k:i\mapsto j$ means that edge $k\in\mathcal{E}$ is directed from node
$i\in\mathcal{V}$ to node $j\in\mathcal{V}$. For each node $j\in\mathcal{V}$,
we define (potentially empty) incoming and outgoing sets of pipelines by
${}_{\mapsto}j=\\{k\in\mathcal{E}|k:i\mapsto j\\}$ and
$j_{\mapsto}=\\{k\in\mathcal{E}|k:j\mapsto i\\}$, respectively. All
nomenclature is listed in Appendix C.

###  II.1 Modeling Physical Flow in a Pipe

Compressible flow of a homogeneous ideal gas through a pipe is described using
the one-dimensional isothermal Euler equations [31],

$\displaystyle\partial_{t}\rho+\partial_{x}(\rho u)$ $\displaystyle=0,$ (1a)
$\displaystyle\partial_{t}(\rho u)+\partial_{x}(p+\rho u^{2})$
$\displaystyle=-\frac{\lambda}{2D}\rho u|u|-\rho g\frac{\partial h}{\partial
x},$ (1b) $\displaystyle p=\rho ZR\mathbf{T}$ $\displaystyle=\sigma^{2}\rho,$
(1c)

where the variables $u(t,x)$, $p(t,x)$, and $\rho(t,x)$ represent velocity,
pressure, and density of the gas, respectively. Here, $t\in[0,T]$ and
$x\in[0,\ell]$, where $T$ denotes the time horizon and $\ell$ denotes the
length of the pipe. The symbols $\partial_{t}$ and $\partial_{x}$ denote the
differential operators with respect to time $t$ and axial location $x$,
respectively. The above system describes mass conservation (1a), momentum
conservation (1b), and the gas equation of state (1c). The variable $h$
represents the elevation of the pipe. The dominant term in the momentum
equation (1b) is the phenomenological Darcy-Weisbach term that models momentum
loss caused by turbulent friction, and is scaled by a dimensionless parameter
$\lambda$ called the friction factor. The remaining parameters are the
internal pipe diameter $D$, the wave (sound) speed
$\sigma=\sqrt{ZR\mathbf{T}}$ in the gas, and the gravitational acceleration
$g$, where $Z$, $R$, and $\mathbf{T}$ are the gas compressibility factor,
specific gas constant, and absolute temperature, respectively. Here, we assume
that gas pressure $p$ and gas density $\rho$ satisfy the ideal gas equation of
state (1c) with wave speed $\sigma$. While non-ideal modeling is necessary in
practice to correctly quantify flows at pressures used in large gas transport
pipelines, ideal gas modeling still qualitatively captures the flow
phenomenology, so we use it for simplicity of exposition. Extension to non-
ideal gas modeling can be made by applying appropriate nonlinear transforms
[30].

It is standard to use the per area mass flux $\varphi=\rho u$, and assume that
gas flow is an isothermal process, that flow is turbulent and has high
Reynolds number, and that the flow is adiabatic, i.e. there is no heat
exchange with the ground [32]. For slowly varying boundary conditions, the
kinetic energy term $\partial_{x}(\rho u^{2})$ and the inertia term
$\partial_{t}(\rho u)$ in equation (1b) may be omitted [31]. With these
assumptions, and given no elevation changes, Eq. (1) can be reduced to

$\displaystyle\partial_{t}\rho+\partial_{x}\varphi$ $\displaystyle=0,$ (2a)
$\displaystyle\partial_{x}(\sigma^{2}\rho)$
$\displaystyle=-\frac{\lambda}{2D}\frac{\varphi|\varphi|}{\rho}.$ (2b) where
$\rho$ and $\varphi$ denote density and mass flux (in per-area units). The
above set of equations have been used in several previous studies [33, 8], and
we refer the reader there for further justifications. Here, we extend these
equations to the case of a mixture of two constituent gases, whose partial
pressures, partial densities, partial fluxes, and mass fractions are denoted
by $p^{(m)}$, $\rho^{(m)}$, $\varphi^{(m)}$, and $\eta^{(m)}$, respectively,
where $m=1$ and $m=2$ are used to identify the two distinct gases. The
fraction of mass of each gas is related to the partial density variables by
$\eta^{(m)}=\rho^{(m)}/(\rho^{(1)}+\rho^{(2)})$. The propagation of either
mass fraction quantity $\eta^{(m)}$ can be modeled by the convection-diffusion
equation with diffusion terms omitted [10], i.e.,
$\displaystyle\partial_{t}\eta^{(m)}+\frac{\varphi}{\rho}\partial_{x}\eta^{(m)}=0.$
(2c)

It follows from the relation $\eta^{(1)}=(1-\eta^{(2)})$ that only one mass
fraction variable should be modeled if there are exactly two distinct gases in
the mixture. It can be shown that $\eta^{(1)}$ solves Eq. (2c) if and only if
$\eta^{(2)}$ does. Partial densities $\rho^{(m)}$ and partial fluxes
$\varphi^{(m)}$ are proportional to their total counterparts
$\rho=(\rho^{(1)}+\rho^{(2)})$ and $\varphi=(\varphi^{(1)}+\varphi^{(2)})$
with the mass fraction being the coefficient of proportionality. In
particular, $\rho^{(m)}=\eta^{(m)}\rho$ and $\varphi^{(m)}=\eta^{(m)}\varphi$.
However, proportionality according to mass fraction does not generally hold
true for pressure. The ideal equation of state for each constituent gas is
defined by $p^{(m)}=\sigma_{m}^{2}\rho^{(m)}$, where $\sigma_{1}$ and
$\sigma_{2}$ are the wave speeds of the two gases. Using Dalton’s law [34],
the total pressure of the mixture is defined to be the summation of partial
pressures given by

$p=p^{(1)}+p^{(2)}=\sigma_{1}^{2}\rho^{(1)}+\sigma_{2}^{2}\rho^{(2)}=\sigma^{2}\rho,$

where $\sigma^{2}=(\sigma_{1}^{2}\eta^{(1)}+\sigma_{2}^{2}\eta^{(2)})$. It
follows that the local wave speed $\sigma$ of the mixture depends on the local
mass fraction of the gases. Moreover, the total pressure and partial pressure
variables are related through the volumetric concentration defined by
$\nu^{(m)}=\sigma_{m}^{2}\rho^{(m)}/(\sigma_{1}^{2}\rho^{(1)}+\sigma_{2}^{2}\rho^{(2)})$.
In particular, $p^{(m)}=\nu^{(m)}p$. From here onward, we use the terms mass
fraction and concentration interchangeably, and specifically refer to
_volumetric_ concentration where that quantity is examined.

In our application of energy transport, superscripts “(1)” and “(2)” will
henceforth identify correspondence of variables to natural gas and hydrogen,
respectively. The wave speeds of natural gas and hydrogen are defined,
respectively, by $\sigma_{1}$ and $\sigma_{2}$. We note that natural gas is
itself a mixture composed of several distinct gases, and its composition has
historically been modeled as constant and homogeneous in the majority of
academic pipeline simulation studies. This simplification is accepted because
the composition of processed pipeline quality gas is at least 90% methane and
typically 95% methane by molar mass, and with ethane comprising the majority
of the other constituents [35]. The calorific content of methane from
different sources may vary somewhat [36], and numerical methods for pipeline
simulation that account for variation in gas composition are well developed
[10]. In contrast, the calorific values and molecular masses of hydrogen and
natural gas both differ by an order of magnitude. In this study, we suppose
that natural gas has a nominal homogeneous and constant composition, and focus
on the phenomenology of dynamic pipeline response to injection of much less
dense hydrogen gas with variation in time and by location. In general,
although our focus here is on blending two gases in a network of pipelines,
the models and theory here can be extended to mixtures of more than two gases
with additional equations and variables corresponding to superscripts “($m$)”
for $m=1,\dots,M$, where $M$ would denote the total number of gas components
of the mixture.

###  II.2 Gas Mixture Dynamics on a Network

With the above assumptions, the flow dynamics through the horizontal pipeline
of index $k\in\mathcal{E}$ is modeled with the friction-dominated PDEs

$\displaystyle\partial_{t}\rho_{k}^{(m)}+\partial_{x}\left(\frac{\rho_{k}^{(m)}}{\rho_{k}^{(1)}+\rho_{k}^{(2)}}\varphi_{k}\right)$
$\displaystyle=$ $\displaystyle 0,$ (3)
$\displaystyle\partial_{x}\left(\sigma^{2}_{1}\rho_{k}^{(1)}+\sigma^{2}_{2}\rho_{k}^{(2)}\right)$
$\displaystyle=$
$\displaystyle-\frac{\lambda_{k}}{2D_{k}}\frac{\varphi_{k}|\varphi_{k}|}{\rho_{k}^{(1)}+\rho_{k}^{(2)}},\quad$
(4)

where Eq. (3) is defined for both $m=1$ and $m=2$. We leave it to the reader
to verify that Eqs. (3)-(4) defined in terms of partial densities $\rho^{(1)}$
and $\rho^{(2)}$ and total flow $\varphi$ are equivalent to Eq. (2) defined in
terms of total density $\rho$, total flow $\varphi$, and one concentration
variable $\eta^{(2)}$. The wave speeds, $\sigma_{1}$ and $\sigma_{2}$, and the
diameter $D_{k}$, length $\ell_{k}$, and friction factor $\lambda_{k}$
associated with each pipeline $k\in\mathcal{E}$ are the parameters of the
system.

Compressor and regulator stations are critical components that actuate the
flow of gas through the network and reduce pressure in the direction of flow,
respectively. For convenience, we assume that a compressor is located at the
inlet and a regulator is located at the outlet of each pipeline, where inlet
and outlet are defined with respect to the oriented positive flow direction.
For each pipeline $k\in\mathcal{E}$, compression and regulation are modeled
with multiplicative control variables $\underline{\mu}_{k}(t)\geq 1$ and
$\overline{\mu}_{k}(t)\geq 1$, respectively. That is, discharge pressure and
density are $\underline{\mu}_{k}(t)$ times larger than the suction pressure
and density.

The boundary conditions for a mixture of gases allow for more degrees of
freedom than those for a single gas, and are formulated here to enable
definition of a range of potential scenarios. All of the flow quantities
defined in this paragraph are, in general, time-varying, but we suppress time-
dependence for readability. The network nodes are partitioned into slack nodes
$\mathcal{V}_{s}\subset\mathcal{V}$ and non-slack nodes
$\mathcal{V}_{d}\subset\mathcal{V}$. Slack (pressure) nodes are typically used
to represent large sources of gas, at which pressure and concentration are
specified and inflow to the network is a dependent variable that is determined
by solving the network flow equations defined below. At non-slack (flow)
nodes, the withdrawal flow from the network is specified, and the pressure and
mass fraction is determined by solving the initial boundary value problem.
Alternatively, injection can be specified at a non-slack (flow) node if
concentration is also provided. Slack nodes are assumed to be ordered in
$\mathcal{V}$ before non-slack nodes, so that $i<j$ for all
$i\in\mathcal{V}_{s}$ and $j\in\mathcal{V}_{d}$. A mixture of gas is injected
into the network at each slack node $i\in\mathcal{V}_{s}$. The boundary
conditions at the slack nodes $i\in\mathcal{V}_{s}$ are defined by specifying
individual densities $\bm{s}_{i}^{(1)}$ and $\bm{s}_{i}^{(2)}$. Alternatively,
pressure $(\bm{p}_{s})_{i}$ and concentration $\bm{\alpha}_{i}^{(m)}$ may be
specified at slack nodes $i\in\mathcal{V}_{s}$. The relations
$(\bm{p}_{s})_{i}=(\sigma_{1}^{2}\bm{s}_{i}^{(1)}+\sigma_{2}^{2}\bm{s}_{i}^{(2)})$
and
$\bm{\alpha}_{i}^{(m)}=\bm{s}_{i}^{(m)}/(\bm{s}_{i}^{(1)}+\bm{s}_{i}^{(2)})$
can then be used to determine the corresponding partial densities that will
achieve the specified pressures and concentrations. Non-slack nodes are
partitioned into injection nodes $\mathcal{V}_{q}\subset\mathcal{V}_{d}$ and
withdrawal nodes $\mathcal{V}_{w}\subset\mathcal{V}_{d}$. We order the non-
slack nodes $\mathcal{V}_{d}$ with injection nodes enumerated before
withdrawal nodes, so that $i<j$ for all $i\in\mathcal{V}_{q}$ and
$j\in\mathcal{V}_{w}$. A mixture is withdrawn from the network at each
withdrawal node $j\in\mathcal{V}_{w}$ with boundary conditions specified by
mass outflow $\bm{w}_{j}\geq 0$. At each injection node $j\in\mathcal{V}_{q}$,
a mixture is injected into the network with boundary conditions specified by
both the mass inflow $\bm{q}_{j}$, with $\bm{q}_{j}\geq 0$, and the
concentration $\bm{\beta}_{j}^{(m)}$. Although a mass inflow is specified at
each injection node $j\in\mathcal{V}_{q}$ with concentration
$\bm{\beta}_{j}^{(m)}$, this does not, in general, imply that the
concentration flowing from node $j$ into outgoing edges is equal to
$\bm{\beta}_{j}^{(m)}$, because the nodal concentration is a mixture of flows
entering node $j$ either by injection or from incoming pipelines. Boundary
condition designations are illustrated for a small example network in Fig. 1.

Figure 1: Configuration of the boundary conditions. Here,
$j\in\mathcal{V}_{w}$, $k_{1}:i_{1}\mapsto j$ with $i_{1}\in\mathcal{V}_{s}$,
and $k_{2}:i_{2}\mapsto j$ with $i_{2}\in\mathcal{V}_{q}$.

For each non-slack node $j\in\mathcal{V}_{d}$, the nodal partial density
variables are denoted by $\bm{\rho}_{j}^{(m)}$ and the dependent concentration
variables are denoted by
$\bm{\eta}_{j}^{(m)}=\bm{\rho}_{j}^{(m)}/(\bm{\rho}_{j}^{(1)}+\bm{\rho}_{j}^{(2)})$
for $m=1$ and $m=2$. All of the nodal quantities in this study are identified
with bold symbols. Inlet and outlet edge variables are defined by attaching
underlines below and overlines above the associated edge variables,
respectively. For example, $\underline{\varphi}_{k}(t)=\varphi_{k}(t,0)$ and
$\overline{\varphi}_{k}(t)=\varphi_{k}(t,\ell_{k})$. Let us denote the cross-
sectional area of edge $k\in\mathcal{E}$ by $\chi_{k}=\pi D_{k}^{2}/4$. The
boundary conditions for the flow of the mixture are defined for $m=1$ and
$m=2$ by

$\displaystyle\underline{\rho}_{k}^{(m)}$ $\displaystyle=$
$\displaystyle\underline{\mu}_{k}\bm{s}_{i}^{(m)},\qquad\overline{\rho}_{k}^{(m)}=\overline{\mu}_{k}\bm{\rho}_{j}^{(m)},$
(5) $\displaystyle\underline{\rho}_{k}^{(m)}$ $\displaystyle=$
$\displaystyle\underline{\mu}_{k}\bm{\rho}_{i}^{(m)},\qquad\overline{\rho}_{k}^{(m)}=\overline{\mu}_{k}\bm{\rho}_{j}^{(m)},$
(6) $\displaystyle\bm{\gamma}_{j}^{(m)}\bm{d}_{j}$ $\displaystyle=$
$\displaystyle\sum_{k\in_{\mapsto}j}\chi_{k}\overline{\eta}_{k}^{(m)}\overline{\varphi}_{k}-\sum_{k\in
j_{\mapsto}}\chi_{k}\underline{\eta}_{k}^{(m)}\underline{\varphi}_{k},$ (7)

where Eq. (5) is defined for $k:i\mapsto j$ with $i\in\mathcal{V}_{s}$, Eq.
(6) is defined for $k:i\mapsto j$ with $i,j\in\mathcal{V}_{d}$, and Eq. (7) is
defined for $j\in\mathcal{V}_{d}$ with the condition that
$\bm{\gamma}_{j}^{(m)}\bm{d}_{j}=\bm{\eta}_{j}^{(m)}\bm{w}_{j}$ if
$j\in\mathcal{V}_{w}$ and
$\bm{\gamma}_{j}^{(m)}\bm{d}_{j}=-\bm{\beta}_{j}^{(m)}\bm{q}_{j}$ if
$j\in\mathcal{V}_{q}$. The conditions in Eqs. (5)-(6) are multiplicative
relations between the partial densities at the boundaries of pipelines and the
auxiliary “nodal” variables for partial densities that are internal to the
nodes at each end (see Fig. 1). The multiplicative factors are used as control
variables that represent gas compression and regulation. The condition for the
conservation of mass through each non-slack node in Eq. (7) depends on whether
the node is an injection node or a withdrawal node. For injection, we specify
both the concentration of the mixture and the mass flow being injected into
the node, but for withdrawal, we specify only the mass outflow. The
concentration of the withdrawn gas is determined by solving the network flow
equations.

The initial conditions of partial density are assumed to be a steady-state
solution given for all $k\in\mathcal{E}$ and $x\in[0,\ell_{k}]$ by

$\displaystyle\rho_{k}^{(m)}(0,x)=\varrho_{k}^{(m)}(x).$ (8)

The steady-state configuration is defined to be the solution of the system in
Eqs. (3)-(7) when the boundary condition profiles are time-invariant (i.e.
equal to the initial values of the time-varying boundary profiles). More
details on the initial condition for the discretized system are provided in
the following section. Mass flux is not initially specified because it is
uniquely determined from the density. We assume standard conditions for well-
posedness [37], and specifically that the boundary conditions are smooth,
slowly-varying, bounded in their respective domains, and compatible with the
initial conditions to ensure the existence of a smooth, slowly-varying,
bounded solution. The flow of the mixture of gases in the network is defined
by the initial-boundary value system of PDEs defined by equations (3)-(8).

##  III Spatial Discretization

To analyze the system of PDEs (3)-(8) on the graph
$(\mathcal{E},\mathcal{V})$, we have developed a process of discretization,
which includes a _refinement_ of the graph, _approximation_ of the PDE system
by an ODE system using a finite volume approach, and a _reformulation_ in
terms of variable vectors and parameter matrices. The vectors include
variables that represent the states and boundary condition profiles, and the
matrices incorporate network model parameters, the incidence structure of the
graph, and the control values.

Graph Refinement: A refinement $(\hat{\mathcal{E}},\hat{\mathcal{V}})$ of the
graph $(\mathcal{E},\mathcal{V})$ is created by adding auxiliary nodes to
$\mathcal{V}$ in order to subdivide the edges of $\mathcal{E}$ so that
$\ell_{k}\leq\ell$ for all $k\in\hat{\mathcal{E}}$, where $\ell$ is
sufficiently small [18]. Henceforth, we assume that $\ell\leq 1$ (km), and
will use that threshold for computational studies as well. The refined graph
inherits the prescribed orientation of the parent graph. Assuming sufficiently
fine network refinement, the relative difference of the density variables of
adjacent nodes in the solution to the IBVP (3)-(8) can be made arbitrarily
small in magnitude because of continuity of the solution to the system given
well-posed conditions [37]. We assume for all $k\in\hat{\mathcal{E}}$ that

$\frac{\left|\overline{\rho}_{k}^{(m)}-\underline{\rho}_{k}^{(m)}\right|}{\underline{\rho}_{k}^{(m)}}<\epsilon,\qquad\frac{\left|\overline{\rho}_{k}^{(m)}-\underline{\rho}_{k}^{(m)}\right|}{\overline{\rho}_{k}^{(m)}}<\epsilon,$
(9)

where $0\leq\epsilon\ll 1$. The proofs that follow only require $\epsilon\leq
1$. We assume that the graph has been sufficiently refined to satisfy Eq. (9)
and that the hats may be omitted moving forward.

Finite Volume Approximation: The system of ODEs is obtained by integrating the
dynamics in Eqs. (3)-(4) along the length of each refined pipeline segment so
that

$\displaystyle\int_{0}^{\ell}\partial_{t}\rho^{(m)}dx=-\int_{0}^{\ell}\partial_{x}\left(\frac{\rho^{(m)}}{\rho^{(1)}+\rho^{(2)}}\varphi\right)dx,$
$\displaystyle\int_{0}^{\ell}\partial_{x}\left(\sigma^{2}_{1}\rho^{(1)}+\sigma^{2}_{2}\rho^{(2)}\right)dx=-\frac{\lambda}{2D}\int_{0}^{\ell}\frac{\varphi|\varphi|}{\rho^{(1)}+\rho^{(2)}}dx,$

where edge subscripts have been removed for readability. The above integrals
of space derivatives are evaluated using the fundamental theorem of calculus.
The remaining integrals are evaluated by approximating pipeline density by
outlet density and pipeline flux by inlet flux. These approximations are
independent of $x$ and may be factored out of the integrals. The above
equations become

$\displaystyle\ell\dot{\overline{\rho}}^{(m)}$ $\displaystyle=$
$\displaystyle\underline{\eta}^{(m)}\underline{\varphi}-\overline{\eta}^{(m)}\overline{\varphi},$
(10)
$\displaystyle\sum_{n=1}^{2}\sigma_{n}^{2}\left(\overline{\rho}^{(n)}-\underline{\rho}^{(n)}\right)$
$\displaystyle=$
$\displaystyle-\frac{\lambda\ell}{2D}\frac{\underline{\varphi}\left|\underline{\varphi}\right|}{\overline{\rho}^{(1)}+\overline{\rho}^{(2)}},$
(11)

where a dot above a variable represents the time-derivative of the variable.

Matrix Form: We now write the discretized system in matrix-vector form. Define
$E\times E$ diagonal matrices $L$ and $X$ with diagonal entries
$L_{kk}=\ell_{k}$ and $X_{kk}=\chi_{k}$. Define the time-varying (transposed)
incidence matrix $M$ of size $E\times V$ componentwise by

$\displaystyle M_{kj}$
$\displaystyle=\begin{cases}\overline{\mu}_{k}(t),&\text{ edge
$k\in_{\mapsto}j$ enters node $j$},\\\ -\underline{\mu}_{k}(t),&\text{ edge
$k\in j_{\mapsto}$ leaves node $j$,}\\\ 0,&\text{else.}\end{cases}$ (12)

Define the $E\times V_{s}$ submatrix $M_{s}$ of $M$ by the removal of columns
$i\in\mathcal{V}_{d}$, the $E\times(V-V_{s})$ submatrix $M_{d}$ of $M$ by the
removal of columns $i\in\mathcal{V}_{s}$, and the positive and negative parts
of $M_{d}$ by $\overline{M}_{d}$ and $\underline{M}_{d}$ so that
$M_{d}=(\overline{M}_{d}+\underline{M}_{d})/2$ and
$|M_{d}|=(\overline{M}_{d}-\underline{M}_{d})/2$, where $V_{s}$ denotes the
number of slack nodes and $|A|$ denotes the componentwise absolute value of a
matrix $A$. Define the signed matrices $Q_{d}=$sign$(M_{d})$,
$\overline{Q}_{d}=$sign$(\overline{M}_{d})$,
$\underline{Q}_{d}=$sign$(\underline{M}_{d})$, and similarly for $M_{s}$.
These signed matrices are well-defined by the lower-bound constraints on
compression and regulation. The incidence matrix provides a compact
representation of the network structure that can be used to specify the
dynamical equations for all of the refined edges. Moreover, the signed
incidence matrix is used to produce the first order finite difference
approximation of the spatial derivative of each component of the vector on
which it acts. The absolute values of the positive and negative parts of the
weighted incidence matrix apply the multiplicative boost in compression and
reduction in regulation, respectively, to the vectors on which they act.

We define the $V_{d}\times V_{d}$ identity matrix $I$, the $V_{d}\times V_{q}$
submatrix $I_{q}$ of $I$ by the removal of columns $j\in\mathcal{V}_{w}$, and
the $V_{d}\times V_{d}$ matrix $I_{w}$ by replacing columns
$j\in\mathcal{V}_{q}$ of $I$ with the zero vector. Here, $V_{d}$ and $V_{q}$
denote the numbers of non-slack nodes and non-slack injection nodes,
respectively. We also define inlet and outlet edge mass flux vectors by
$\underline{\varphi}=(\underline{\varphi}_{1},\dots,\underline{\varphi}_{E})^{T}$
and
$\overline{\varphi}=(\overline{\varphi}_{1},\dots,\overline{\varphi}_{E})^{T}$,
and similarly for inlet and outlet edge concentrations $\underline{\eta}$ and
$\overline{\eta}$. Then, let us define the vectors
$\bm{\rho}^{(m)}=(\bm{\rho}_{V_{s}+1}^{(m)},\dots,\bm{\rho}_{V_{d}}^{(m)})^{T}$,
$\bm{\alpha}^{(m)}=(\bm{\alpha}_{1}^{(m)},\dots,\bm{\alpha}_{V_{s}}^{(m)})^{T}$,
and
$\bm{\beta}^{(m)}=(\bm{\beta}^{(m)}_{V_{s}+1},\dots,\bm{\beta}^{(m)}_{V_{q}})^{T}$,
where the subscripts of the vector components are indexed according to the
node labels in $\mathcal{V}$. Similarly, define the vectors
$\bm{\eta}^{(m)}=(\bm{\eta}^{(m)}_{V_{s}+1},\dots,\bm{\eta}^{(m)}_{V_{d}})^{T}$
and $\bm{d}=(\bm{d}_{V_{s}+1},\dots,\bm{d}_{V_{d}})^{T}$. Recall that the
components of $\bm{d}$ are positive for those corresponding to non-slack
withdrawal nodes and negative for non-slack injection nodes. Define the
function $f:\mathbb{R}^{E}\times\mathbb{R}^{E}\rightarrow\mathbb{R}^{E}$
component-wise for $k\in\mathcal{E}$ by

$\displaystyle f_{k}(y,z)$
$\displaystyle=-\text{sign}(z_{k})\Lambda_{k}\left|y_{k}z_{k}\right|^{1/2},$
(13)

where $\Lambda_{k}=\sqrt{2D_{k}/(\lambda_{k}\ell_{k})}$. This function is used
to express $\underline{\varphi}$ in Eq. (11) in terms of density and its
spatial derivative so that we may eliminate flux from the dynamic equations.
Using the function in Eq. (13), the discretized flow in Eqs. (10)-(11)
together with the boundary conditions in Eqs. (5)-(7) may be expressed in
matrix-vector form as

$\displaystyle\\!\\!\\!\\!L\overline{M}_{d}\dot{\bm{\rho}}^{(m)}$
$\displaystyle=$ $\displaystyle\underline{\eta}^{(m)}\odot
F-\overline{\eta}^{(m)}\odot\overline{\varphi},$ (14)
$\displaystyle\\!\\!\\!\\!\bm{\gamma}^{(m)}\odot\bm{d}$ $\displaystyle=$
$\displaystyle\overline{Q}_{d}^{T}X\left(\overline{\eta}^{(m)}\odot\overline{\varphi}\right)+\underline{Q}_{d}^{T}X\left(\underline{\eta}^{(m)}\odot
F\right)\\!,\quad$ (15)

where $\odot$ is the Hadamard product and

$\displaystyle\\!\\!F\\!$
$\displaystyle=\\!f\\!\left(\\!\overline{M}_{d}(\bm{\rho}^{(1)}\\!+\\!\bm{\rho}^{(2)}),\sum_{m}\\!\sigma_{m}^{2}(M_{s}\bm{s}^{(m)}\\!+\\!M_{d}\bm{\rho}^{(m)})\\!\right)\\!.$
(16)

We suppose that regulators vary slowly so the time derivative of
$\overline{M}_{d}$ is insignificant, justifying its removal from Eq. (14).

The interested reader is encouraged to sketch a small network such as the one
in Fig. 1; label the edges, nodes, compressors, and regulators; and construct
the associated weighted incidence matrix, its submatrices, and their signed
correspondences. Upon substituting these matrices into Eq. (14) and performing
the matrix multiplications, we should observe that each component of the
resulting vector equation is simply Eq. (10) corresponding to the edge
component. Likewise, the components of Eq. (16) reduce to Eq. (11) after some
algebraic rearrangements. Multiplying both sides of Eq. (14) on the left by
$\overline{Q}_{d}^{T}X$ and using Eq. (15), we may combine Eq. (14) and Eq.
(15) to form the equation
$\overline{Q}_{d}^{T}XL\overline{M}_{d}\dot{\bm{\rho}}^{(m)}=[Q_{d}^{T}X(\underline{\eta}^{(m)}\odot
F)-\bm{\gamma}^{(m)}\odot\bm{d}]$, where we have used
$Q_{d}=(\underline{Q}_{d}+\overline{Q}_{d})$. By writing edge concentrations
in terms of nodal concentrations, and nodal concentrations in terms of
concentrations of flows into the nodes, the system in Eqs. (14)-(16) may be
written for both $m=1$ and $m=2$ as

$\displaystyle R\dot{\bm{\rho}}^{(m)}$ $\displaystyle=$ $\displaystyle
Q_{d}^{T}X\left(\left(|\underline{Q}_{s}|\bm{\alpha}^{(m)}+|\underline{Q}_{d}|\bm{\eta}^{(m)}\right)\odot
F\right)$ (17)
$\displaystyle\qquad-\left(I_{q}\bm{\beta}^{(m)}+I_{w}\bm{\eta}^{(m)}\right)\odot\bm{d},\quad$

where $R=\overline{Q}_{d}^{T}XL\overline{M}_{d}$. The system in Eq. (17) will
be called the partial density system of ODEs. Each row $k$ of
$\overline{M}_{d}$ contains exactly one nonzero component given by
$\overline{M}_{kj}=\overline{\mu}_{k}$ for $k\in$${}_{\mapsto}j$. Using the
additional fact that $X$ and $L$ are diagonal, it can be shown that the mass
matrix $R$ on the left-hand-side of Eq. (17) is diagonal with positive
diagonal components given by
$\bm{r}_{j}=\sum_{k\in_{\mapsto}j}\chi_{k}\ell_{k}\overline{\mu}_{k}$ for
$j\in\mathcal{V}_{d}$. Therefore, the matrix $R$ may readily be inverted to
obtain a nonlinear control system in the usual, although complicated, ODE
form. The initial condition in Eq. (8), sampled at the refined nodes of the
network, is the time-invariant solution of the system in Eq. (17) with
$\bm{d}=\bm{d}(0)$, $\bm{\alpha}^{(m)}=\bm{\alpha}^{(m)}(0)$, and
$\bm{\beta}^{(m)}=\bm{\beta}^{(m)}(0)$. We assume that this steady-state
solution is the initial condition of the partial density system.

##  IV Equivalent Systems

The system in Eq. (17) is expressed in terms of partial densities at non-slack
nodes. Equivalent systems expressed in terms of other variables of interest
may be derived from Eq. (17) using appropriate transformations, such as that
performed in the continuous case going from Eqs. (2a)-(2c) to Eqs. (3)-(4). In
fact, such transformations exist even for homogeneous gas systems. For
example, the equations of natural gas flow may be expressed in terms of
pressure and velocity, in terms of density and mass flux, or in terms of their
dimensionless quantities. Define vectors $\bm{\rho}$, $\bm{p}$,
$\bm{\nu}^{(m)}$, and $\bm{E}$ of nodal values for density, pressure,
volumetric concentration, and energy, respectively, at non-slack nodes by

$\displaystyle\bm{\rho}$ $\displaystyle=$
$\displaystyle\bm{\rho}^{(1)}+\bm{\rho}^{(2)},$ (18) $\displaystyle\bm{p}$
$\displaystyle=$
$\displaystyle\sigma_{1}^{2}\bm{\rho}^{(1)}+\sigma_{2}^{2}\bm{\rho}^{(2)},$
(19) $\displaystyle\bm{\nu}^{(m)}$ $\displaystyle=$
$\displaystyle\frac{\sigma_{m}^{2}\bm{\rho}^{(m)}}{\sigma_{1}^{2}\bm{\rho}^{(1)}+\sigma_{2}^{2}\bm{\rho}^{(2)}},$
(20) $\displaystyle\bm{E}$ $\displaystyle=$
$\displaystyle(|\overline{Q}_{d}^{T}|X\underline{\varphi})\odot\left(r^{(1)}\bm{\eta}^{(1)}+r^{(2)}\bm{\eta}^{(2)}\right),$
(21)

where $r^{(1)}=44.2$ (MJ/kg) and $r^{(2)}=141.8$ (MJ/kg). Equivalent systems
may be expressed in terms of any two vector variables from the set
$\\{\bm{\rho}^{(m)},\bm{\eta}^{(m)},\bm{\nu}^{(m)},\bm{\rho},\bm{p},\bm{E}\\}$,
excluding pairs from the subset
$\\{\bm{\eta}^{(1)},\bm{\eta}^{(2)},\bm{\nu}^{(1)},\bm{\nu}^{(2)}\\}$ because
variables in the latter subset would reduce to constant vectors in the case of
homogeneous mixtures. The choice of which equivalent system to use may depend
on the sought application, although some systems have better conditioning with
fewer nonlinear operations than others. Define (potentially time-varying)
localized wave speed vectors
$\bm{a}=(\sigma_{1}^{2}\bm{\alpha}^{(1)}+\alpha_{2}^{2}\bm{\alpha}^{(2)})^{1/2}$
and
$\bm{b}=(\sigma_{1}^{2}\bm{\beta}^{(1)}+\sigma_{2}^{2}\bm{\beta}^{(2)})^{1/2}$,
where the square-root is applied component-wise. The transformation from
partial densities to _total density and pressure_ is obtained by superimposing
Eq. (17) for $m=1,2$ to obtain an equation for $\dot{\bm{\rho}}$ and linearly
combining Eq. (17) for $m=1,2$ with coefficients $\sigma_{1}^{2}$ and
$\sigma_{2}^{2}$ to obtain an equation for $\dot{\bm{p}}$. This transformation
produces the system

$\displaystyle R\dot{\bm{\rho}}$ $\displaystyle=$ $\displaystyle
Q_{d}^{T}XF-\bm{d},$ (22) $\displaystyle R\dot{\bm{p}}$ $\displaystyle=$
$\displaystyle
Q_{d}^{T}X\left(\left(|\underline{Q}_{s}|\bm{a}^{2}+|\underline{Q}_{d}|\frac{\bm{p}}{\bm{\rho}}\right)\odot
F\right)$ (23)
$\displaystyle\qquad\quad-\left(I_{q}\bm{b}^{2}+I_{w}\frac{\bm{p}}{\bm{\rho}}\right)\odot\bm{d},\quad$

where $F=f(\overline{M}_{d}\bm{\rho},M_{s}\bm{p}_{s}+M_{d}\bm{p})$. The system
in Eqs. (22)-(23) will be called the total density and pressure system of
ODEs. We do not derive other equivalent systems. Instead, we compute the
solution of the partial density system of ODEs numerically, and, thereafter,
obtain the other variables of interest by subsequently applying the
appropriate transformations.

If $\bm{\eta}^{(m)}$ is a constant vector, then the system of total density
and pressure decouples into two isolated subsystems that are equivalent to one
another because, for constant concentration, we have
$\bm{p}=\bm{c}^{2}\odot\bm{\rho}$, where
$\bm{c}=(\sigma_{1}^{2}\bm{\eta}^{(1)}+\sigma_{2}^{2}\bm{\eta}^{(2)})^{1/2}$
is a constant vector. In particular, the total density and pressure system in
Eqs. (22)-(23) reduces by half its dimension to the isolated system

$\displaystyle R\dot{\bm{p}}$ $\displaystyle=$ $\displaystyle
Q_{d}^{T}X\left(\left(|\underline{Q}_{s}|\bm{a}^{2}+|\underline{Q}_{d}|\bm{c}^{2}\right)\right.$
(24)
$\displaystyle\qquad\quad\odot\left.f\left(\overline{M}_{d}\frac{\bm{p}}{\bm{c}^{2}},M_{s}\bm{p}_{s}+M_{d}\bm{p}\right)\right)$
$\displaystyle-\left(I_{q}\bm{b}^{2}+I_{w}\bm{c}^{2}\right)\odot\bm{d}.\quad$

The system in Eq. (24) is called the isolated total pressure system of ODEs.
Equivalent isolated subsystems expressed in terms of one vector variable from
the set $\\{\bm{\rho}^{(m)},\bm{\rho},\bm{E}\\}$ may be derived. Each isolated
subsystem is only applicable if the concentration vector $\bm{\eta}^{(m)}$ is
constant. Rigorous definitions and proofs of conditions on
$\bm{\alpha}^{(m)}$, $\bm{\beta}^{(m)}$, $\bm{q}$, $\bm{w}$, and network
topology that would result in $\bm{\eta}^{(m)}$ being constant are outside the
scope of this study.

##  V Monotonicity

The monotonicity properties of solutions to initial boundary value problems
for flows of a homogeneous gas through an actuated transport network were
examined as a means to reduce the complexity of optimization and optimal
control of natural gas networks in the presence of uncertainty [8]. Here, we
examine how such concepts can be extended to the transport of inhomogeneous
gas mixtures, and specifically to characterize the acceptable extent and
variability of hydrogen blending into a natural gas pipeline. We first present
some analytical results before proceeding with numerical simulations in the
next section.

A nonlinear input-to-state initial-value system of ODEs may be generally
expressed as

$\dot{x}=g(x,u,d),\qquad x(0)=y,$ (25)

where $x(t)\in\mathcal{X}\subset\mathbb{R}^{n}$ is the state vector,
$u(t)\in\mathcal{U}\subset\mathbb{R}^{m}$ is the control input vector, and
$d(t)\in\mathcal{D}\subset\mathbb{R}^{r}$ is the parameter input vector
defined for $t\in[0,T]$. It is assumed that the subsets $\mathcal{X}$,
$\mathcal{U}$, and $\mathcal{D}$ are compact and convex and that the function
$g:\mathcal{X}\times\mathcal{U}\times\mathcal{D}\rightarrow\mathcal{X}$ is
Lipschitz in $\mathcal{X}\times\mathcal{U}\times\mathcal{D}$.

Definitions: Suppose that two independent state solutions
$\\{x_{1}(t),x_{2}(t)\\}\subset\mathcal{X}$ exist (and are thus unique because
$g$ is Lipschitz) with initial conditions
$\\{y_{1},y_{2}\\}\subset\mathcal{X}$, and which correspond to the piecewise-
continuous control inputs $\\{u_{1}(t),u_{2}(t)\\}\subset\mathcal{U}$ and
piecewise-continuous parameter inputs
$\\{d_{1}(t),d_{2}(t)\\}\subset\mathcal{D}$ for $t\in[0,T]$. For the given set
of control inputs, the system in Eq. (25) is said to be monotone-ordered with
respect to $d(t)$ if $x_{1}(t)\leq x_{2}(t)$ for $t\in[0,T]$ whenever
$y_{1}\leq y_{2}$ and $d_{1}(t)\leq d_{2}(t)$, where inequalities for vectors
are taken componentwise. In this case, the solution states $x_{1}$ and $x_{2}$
are said to be monotone-ordered. For simplicity, we say that a monotone-
ordered system and a set of monotone-ordered solutions are monotone,
monotonic, and have the property of monotonicity. An $n\times n$ matrix $A$ is
called Metzler if all of its off-diagonal elements are non-negative, i.e.
$A_{ij}\geq 0$ for all $i\not=j$. An $n\times m$ matrix is called non-negative
if all of its entries are non-negative.

Theorem 1 (Monotonicity) [38, 39]: The nonlinear system in Eq. (25) is
monotone if and only if the Jacobian matrices $\partial g/\partial x$ and
$\partial g/\partial d$ are, respectively, Metzler and non-negative almost
everywhere in $\mathcal{X}\times\mathcal{U}\times\mathcal{D}$.

###  V.1 Homogeneous Concentration

The equivalent systems described in Section IV are first reformulated in terms
of the monotone system definitions above. In steady-state [8], the pressure
$\bm{p}$ increases componentwise with decreasing withdrawal $\bm{w}\geq 0$ and
with increasing injection $-\bm{q}\leq 0$. In reference to Eq. (25), we assume
that all non-slack nodes are injection nodes and define the input parameter by
$d=\\{\bm{p}_{s},\bm{d}\\}=\\{\bm{p}_{s},-\bm{q}\\}$.

Proposition 1 (Monotonicity of Total Pressure and Density): Assume that i) all
non-slack nodes are injection nodes; ii) gas flows only in the positive
direction through each edge according to its orientation in the network graph;
iii) pressure is positive in each node; and iv) Eq. (9) is satisfied. Suppose
that the concentration vector $\bm{\eta}^{(2)}$ is constant and that there
exist two state solutions $\bm{p}_{1}$, $\bm{p}_{2}$ of the system in Eq. (24)
with respective initial conditions $\bm{\pi}_{1}$, $\bm{\pi}_{2}$, slack node
pressures $(\bm{p}_{s})_{1}$, $(\bm{p}_{s})_{2}$, and non-slack injection
flows $\bm{q}_{1}$, $\bm{q}_{2}$ for a given fixed set of control inputs
$\\{\underline{\mu},\overline{\mu}\\}$. Here, the vector subscripts denote the
first and second solutions and not the refined nodes. If
$\bm{\pi}_{1}\leq\bm{\pi}_{2}$, $(\bm{p}_{s})_{1}(t)\leq(\bm{p}_{s})_{2}(t)$,
and $\bm{q}_{1}(t)\geq\bm{q}_{2}(t)$ componentwise for all $t\in[0,T]$, then
$\bm{p}_{1}(t)\leq\bm{p}_{2}(t)$. Consequently,
$\bm{\rho}_{1}(t)\leq\bm{\rho}_{2}(t)$, where $\bm{\rho}_{1}$ and
$\bm{\rho}_{2}$ are the total densities of the two solutions.

Proof: Throughout this proof, the state and input subscripts correspond to the
nodes of the refined graph. Because flow is in the positive oriented
direction, it follows from Eq. (4) that
$\underline{\mu}_{k}\bm{p}_{i}(t)>\overline{\mu}_{k}\bm{p}_{j}(t)$ for all
$i,j\in\mathcal{V}$ with $k:i\mapsto j$. Thus, the sign and absolute value
operations in Eq. (13) are unnecessary. The $j$-th state dynamics in Eq. (24)
for $j\in\mathcal{V}_{d}$ may be written as

$\displaystyle\bm{r}_{j}\dot{\bm{p}}_{j}$ $\displaystyle=\sum_{k:i\mapsto
j}\frac{\bm{\sigma}^{2}_{i}\chi_{k}\Lambda_{k}}{\bm{c}_{j}}\left(\overline{\mu}_{k}\bm{p}_{j}\left(\underline{\mu}_{k}\bm{p}_{i}-\overline{\mu}_{k}\bm{p}_{j}\right)\right)^{1/2}$
$\displaystyle-\sum_{k:j\mapsto
i}\frac{\bm{c}^{2}_{j}\chi_{k}\Lambda_{k}}{\bm{c}_{i}}\left(\overline{\mu}_{k}\bm{p}_{i}\left(\underline{\mu}_{k}\bm{p}_{j}-\overline{\mu}_{k}\bm{p}_{i}\right)\right)^{1/2}$
(26) $\displaystyle\qquad\qquad+\bm{b}_{j}^{2}\bm{q}_{j},$

where $\bm{p}_{i}=(\bm{p}_{s})_{i}$ and $\bm{\sigma}_{i}^{2}=\bm{a}_{i}^{2}$
if $i\in\mathcal{V}_{s}$, whereas $\bm{\sigma}_{i}^{2}=\bm{c}_{i}^{2}$ if
$i\in\mathcal{V}_{d}$. It is clear from this expanded form that the function
on the right-hand-side of Eq. (24) is continuously differentiable (hence
Lipschitz) in the state and input variables over the domain of positive flow
and pressure. In reference to Theorem 1, we first show that the state Jacobian
matrix is Metzler, i.e., $\partial\dot{\bm{p}_{j}}/\partial\bm{p}_{i}$ is non-
negative for all $i,\,j\in\mathcal{V}_{d}$ with $i\not=j$. If $i$ and $j$ are
non-adjacent with $i\not=j$, then clearly
$\partial\dot{\bm{p}_{j}}/\partial\bm{p}_{i}=0$. Suppose that $i$ and $j$ are
adjacent with $k:j\mapsto i$. Substituting Eqs. (5)-(6) into Eq. (9) and using
the relation between pressure and partial densities, we can show that
$(\overline{\mu}_{k}\bm{p}_{i}-\underline{\mu}_{k}\bm{p}_{j})>-\overline{\mu}_{k}\bm{p}_{i}$.
Thus, the Jacobian component

$\displaystyle\frac{\partial\dot{\bm{p}_{j}}}{\partial\bm{p}_{i}}$
$\displaystyle=$
$\displaystyle\frac{\bm{c}_{j}^{2}\chi_{k}\Lambda_{k}\overline{\mu}_{k}(2\overline{\mu}_{k}\bm{p}_{i}-\underline{\mu}_{k}\bm{p}_{j})}{2\bm{r}_{j}\bm{c}_{i}(\overline{\mu}_{k}\bm{p}_{i}(\underline{\mu}_{k}\bm{p}_{j}-\overline{\mu}_{k}\bm{p}_{i}))^{1/2}}$
(27)

is positive. Suppose that $i$ and $j$ are adjacent with $k:i\mapsto j$. Then

$\displaystyle\frac{\partial\dot{\bm{p}_{j}}}{\partial\bm{p}_{i}}$
$\displaystyle=$
$\displaystyle\frac{\bm{\sigma}_{i}^{2}\chi_{k}\Lambda_{k}\underline{\mu}_{k}\overline{\mu}_{k}\bm{p}_{j}}{2\bm{r}_{j}\bm{c}_{j}(\overline{\mu}_{k}\bm{p}_{j}(\underline{\mu}_{k}\bm{p}_{i}-\overline{\mu}_{k}\bm{p}_{j}))^{1/2}}>0.$
(28)

Because $j\in\mathcal{V}_{d}$ is arbitrary, it follows that the state Jacobian
matrix is Metzler. We now show that the parameter Jacobian matrix is non-
negative. The above computation can be extended to show that
$\partial\dot{\bm{p}_{j}}/\partial(\bm{p}_{s})_{i}$ is non-negative for
$i\in\mathcal{V}_{s}$. With respect to mass inflow parameters, the Jacobian
components
$\partial\dot{\bm{p}_{j}}/\partial\bm{q}_{i}=\bm{b}_{j}^{2}/\bm{r}_{j}\delta_{i,j}$
are non-negative ($\delta_{i,j}$ is the Kronecker delta). We conclude from
Theorem 1 that the system in Eq. (24) is monotone. Because
$\bm{p}_{j}=\bm{c}_{j}^{2}\bm{\rho}_{j}$ for $j\in\mathcal{V}_{d}$, it follows
that the isolated total density system is monotone as well. $\square$

Corollary 1 (Monotonicity of Equivalent Systems): Assume that the conditions
hold from Proposition 1. Then
$\bm{\rho}_{1}^{(m)}(t)\leq\bm{\rho}_{2}^{(m)}(t)$ componentwise for all
$t\in[0,T]$, where $\bm{\rho}_{1}^{(m)}$ and $\bm{\rho}_{2}^{(m)}$ are the
partial densities of the two solutions.

Proof: The mass fraction $\bm{\eta}^{(m)}$ is constant, therefore it follows
from Proposition 1 that
$\bm{\rho}_{1}^{(m)}=\bm{\eta}^{(m)}\odot\bm{\rho}_{1}\leq\bm{\eta}^{(m)}\odot\bm{\rho}_{2}=\bm{\rho}_{2}^{(m)}$.
$\square$

###  V.2 Heterogeneous Concentration

Proposition 2 (Non-Monotonicity of Total Pressure and Density): Assume that i)
all non-slack nodes are injection nodes; ii) gas flows only in the positive
direction through each edge according to its orientation in the network graph;
and iii) pressure and density are positive in each node. Suppose that, for a
given fixed set of control inputs $\\{\underline{\mu},\overline{\mu}\\}$,
there exist two state solutions $(\bm{\rho},\bm{p})_{1}$,
$(\bm{\rho},\bm{p})_{2}$ of the system in Eqs. (22)-(23) with respective
initial conditions $(\bm{\varrho},\bm{\pi})_{1}$,
$(\bm{\varrho},\bm{\pi})_{2}$, slack inputs $(\bm{\rho}_{s},\bm{p}_{s})_{1}$,
$(\bm{\rho}_{s},\bm{p}_{s})_{2}$, and non-slack mass inflows $\bm{q}_{1}$,
$\bm{q}_{2}$ that satisfy
$(\bm{\varrho},\bm{\pi})_{1}\leq(\bm{\varrho},\bm{\pi})_{2}$,
$(\bm{\rho}_{s}(t),\bm{p}_{s}(t))_{1}\leq(\bm{\rho}_{s}(t),\bm{p}_{s}(t))_{2}$,
and $\bm{q}_{1}(t)\geq\bm{q}_{2}(t)$ componentwise for all $t\in[0,T]$. If
$\bm{\eta}^{(m)}(t)$ is time-varying, then, in general,
$(\bm{\rho}(t),\bm{p}(t))_{1}\not\leq(\bm{\rho}(t),\bm{p}(t))_{2}$ component-
wise for all $t\in[0,T]$.

Proof: Throughout this proof, the state and input subscripts correspond to the
nodes of the refined graph. From Theorem 1, it suffices to show that one
component of the state Jacobian matrix is negative. The $j$-th nodal pressure
dynamics in Eq. (23) may be written as

$\displaystyle\bm{r}_{j}\dot{\bm{p}}_{j}$ $\displaystyle=\sum_{k:i\mapsto
j}\bm{\sigma}^{2}_{i}\chi_{k}\Lambda_{k}\left(\overline{\mu}_{k}\bm{\rho}_{j}\left(\underline{\mu}_{k}\bm{p}_{i}-\overline{\mu}_{k}\bm{p}_{j}\right)\right)^{1/2}$
$\displaystyle-\sum_{k:j\mapsto
i}\frac{\bm{p}_{j}}{\bm{\rho}_{j}}\chi_{k}\Lambda_{k}\left(\overline{\mu}_{k}\bm{\rho}_{i}\left(\underline{\mu}_{k}\bm{p}_{j}-\overline{\mu}_{k}\bm{p}_{i}\right)\right)^{1/2}$
(29) $\displaystyle\qquad\qquad+\bm{b}_{j}^{2}\bm{q}_{j},$

where $\bm{p}_{i}=(\bm{p}_{s})_{i}$, $\bm{\rho}_{i}=(\bm{\rho}_{s})_{i}$, and
$\bm{\sigma}_{i}^{2}=\bm{a}_{i}^{2}$ if $i\in\mathcal{V}_{s}$, and
$\bm{\sigma}_{i}^{2}=\bm{p}_{i}/\bm{\rho}_{i}$ if $i\in\mathcal{V}_{d}$. By
adding a refined edge to the graph if necessary, we assume that there is an
edge $k^{\prime}:i^{\prime}\mapsto j$ with $i^{\prime}\in\mathcal{V}_{d}$. The
Jacobian component corresponding to $\bm{\rho}_{i^{\prime}}$ is given by

$\displaystyle\frac{\partial\dot{\bm{p}_{j}}}{\partial\bm{\rho}_{i^{\prime}}}$
$\displaystyle=$
$\displaystyle-\frac{\chi_{k^{\prime}}\Lambda_{k^{\prime}}}{\bm{r}_{j}}\frac{\bm{p}_{i^{\prime}}}{\bm{\rho}_{i^{\prime}}^{2}}\left(\overline{\mu}_{k^{\prime}}\bm{\rho}_{j}\left(\underline{\mu}_{k^{\prime}}\bm{p}_{i^{\prime}}-\overline{\mu}_{k^{\prime}}\bm{p}_{j}\right)\right)^{1/2},$

which is negative. It follows from Theorem 1 that the system in Eqs. (22)-(23)
is not monotone, regardlesss of Eq. (9). $\square$

##  VI Network Case Study

Figure 2: Network configuration (not to scale). The triangles represent
compressor stations. Pipeline dimensions between nodes: blue to black (20 km),
black to green (70 km), green to maroon (10 km), black to maroon (60 km),
maroon to cyan (80 km). The pipelines have uniform diameter (0.9144 m) and
friction factor (0.01), except for the black to maroon pipeline that has
diameter (0.635 m) and friction factor (0.015). Figure 3: Two solutions (solid
lines vs. dots) at the color-coordinated network nodes in Fig. 2. From top to
bottom, the depicted nodal solutions are total pressure (MPa), total density
(kg/m3), total energy (GJ/s), concentration by mass, and concentration by
volume. The boundary conditions for both solutions are
$(\bm{p}_{s})_{\text{blue}}=5$ MPa,
$\bm{\alpha}_{\text{blue}}^{(2)}(t)=0.01(1+\sin(4\pi t/T))$,
$\bm{\beta}_{\text{green}}^{(2)}(t)=0.125(1+\sin(12\pi t/T))$,
$\bm{q}_{\text{green}}(t)=3$ (kg/s), $\bm{w}_{\text{black}}(t)=60(1-\sin(6\pi
t/T))$ (kg/s), $\mu_{\text{red}}=1.0678$,
$\underline{\mu}_{\text{yellow}}=1.0140$, and
$\underline{\mu}_{\text{purple}}=1.0734$, where $T=60$ (hr). The boundary
condition that differs between the two solutions is
$\bm{w}_{\text{cyan}}(t)=110$ (kg/s) (solid lines) and
$\bm{w}_{\text{cyan}}(t)=130$ (kg/s) (dots).  Figure 4: Same boundary
conditions as in Fig. 3 except for $(\bm{p}_{s})_{\text{blue}}=10$ (MPa).
Figure 5: The boundary conditions for the two solutions are
$(\bm{p}_{s})_{\text{blue}}=5$ (MPa), $\bm{\alpha}_{\text{blue}}^{(2)}(t)=0$,
$\bm{\beta}_{\text{green}}^{(2)}(t)=1$,
$\bm{q}_{\text{green}}(t)=9(1+\sin(6\pi t/T))$ (kg/s),
$\bm{w}_{\text{black}}(t)=100(1-\sin(6\pi t/T))$ (kg/s),
$\mu_{\text{red}}=1.1096$, $\underline{\mu}_{\text{yellow}}=1.0057$, and
$\underline{\mu}_{\text{purple}}=1.1301$, where $T=80$ (hr). The other
boundary condition is $\bm{w}_{\text{cyan}}(t)=130$ (kg/s) (solid lines) and
$\bm{w}_{\text{cyan}}(t)=150$ (kg/s) (dots). Figure 6: Same boundary
conditions as in Fig. 5 except for $(\bm{p}_{s})_{\text{blue}}=10$ (MPa).

We use numerical simulations to examine the effects of time-varying
heterogeneity of a transported mixture on the flow dynamics of equivalent
system variables throughout a network. The simulations are performed for a
test network that was used in a previous study [30], in which the authors
presented a staggered grid discretization method for the numerical solution of
homogeneous natural gas pipeline flow. We refer the reader to Appendix A in
which we show the results of our implementation of the IBVP posed in the
former study in order to verify that we obtain the same solution when no
hydrogen is present. The configuration and dimensions of the network are shown
in Fig. 2. The dark blue node is a slack node at which pressure and
concentration are specified, the black, maroon, and cyan nodes are non-slack
withdrawal nodes, and the green node is a non-slack injection node. The wave
speeds are chosen to be $\sigma_{1}=377$ (m/s) and $\sigma_{2}=2.8\sigma_{1}$.
We simulate several examples to illustrate that some physical quantities may
exhibit fewer crossings than others in certain response regimes, given ordered
boundary parameters. These examples provide insight into which equivalent
system may be more useful for simulation and optimization of gas mixture
dynamics [40] for a particular response regime. Figs. 3-7 show the solutions
of five different examples. Each example computes two solutions that
correspond to slightly different but monotone-ordered boundary conditions. For
each example, the two solutions are depicted in the figures distinctly with
solid lines and dots and with colors coordinated to match the colors of the
nodes of the network diagram in Fig. 2. We now describe the simulation results
for each example.

The first example in Fig. 3 considers time-varying sinusoidal forcing in the
concentration of the gas mixtures being injected into the slack node (blue
node) and the non-slack injection node (green node). In this figure, the total
pressure, density, and energy solutions at the non-slack nodes do not cross,
but the mass and volumetric concentrations do show crossings. The solutions in
Fig. 4 have the same boundary conditions as those in Fig. 3 except for the
slack node pressure. If the slack node pressure is doubled, then the total
density trajectories will exhibit crossings at each non-slack node but the
pressure and energy trajectories still do not cross. In Figs. 5 and 6, the
blue node injects pure natural gas and the green node injects pure hydrogen
with a varying mass inflow profile. As seen in Fig. 5, the pressure and energy
solutions at each node do not show crossings. However, a close examination
shows that the density solutions exhibit crossovers at every node upstream
from the point of hydrogen injection. Moreover, the concentration solutions
show crossings at only the cyan node. The solutions in Fig. 6 have the same
boundary conditions as those in Fig. 5 except for the slack node pressure. If
the slack node pressure is doubled, then the resulting pressure, density, and
energy trajectories will cross at all of the non-slack nodes. Moreover, the
concentration trajectories in this example cross at every node upstream from
the node of hydrogen injection. At nodes downstream the injection of hydrogen,
the concentration of hydrogen is zero, as it ought to be. We note that the
solutions in Figs. 5 and 6 may not be realistic in the current operation of
natural gas pipelines because the concentration of hydrogen reaches very high
levels. However, these figures indicate that flow dynamics may exhibit
similarly rapid transients in pipelines that are upgraded to deliver
significant amounts of hydrogen. All of the solution variables show large
gradient surges in small time intervals.

From the simulations presented thus far, it may appear that the concentration
variables are the most likely of the equivalent system variables to violate
monotonicity. The simulation shown in Fig. 7 demonstrates a counterexample to
this conjecture in which pressure, density, and energy trajectories all
exhibit crossings even though the concentration solutions do not. However,
strictly speaking, density trajectories cross only upstream from the point of
hydrogen injection. The difference between the solid line and dotted solutions
in Fig. 7 is that the solid line represents the solution for homogeneous
natural gas and the dotted solution represents a slightly perturbed solution
that results from a small injection of hydrogen made at the green non-slack
node. The concentration variables associated with the two solutions may be
identical at times in certain network nodes, but there cannot be strict
crossings in the concentration trajectories for this example. The reason is
because the homogeneous natural gas solution corresponds to zero hydrogen
concentration and this is the lower bound that the concentration variables can
achieve.

In each of the examples above, the five edges of the network are discretized
into 240 refined edges with $\ell_{k}=1$ (km) for all $k\in\hat{\mathcal{E}}$
and the simulations are performed using the partial density system of ODEs in
Eq. (17). The equivalent system variables are computed using the
transformations presented in Section IV. Although one kilometer discretization
size is sufficient to demonstrate non-monotonicity for slowly-varying
concentrations, a much smaller discretization size is required to accurately
simulate rapidly-varying concentrations. We note that even the slowly-varying
solutions in Figs. 3-7 show noticeable convergence as the discretization size
is decreased from 1 (km) to 100 (m). For small discretization lengths
($\ell_{k}\leq 100$ (m)), the crossings of the solutions in Figs. 3-7 may be
more pronounced. We will return to this point later on in Sections VIII and IX
when we simulate pipeline flows of highly heterogeneous gas mixtures.

Figure 7: The boundary conditions for the two solutions are
$(\bm{p}_{s})_{\text{blue}}=11$ (MPa), $\bm{\alpha}_{\text{blue}}^{(2)}(t)=0$,
$\bm{\beta}_{\text{green}}^{(2)}(t)=1$,
$\bm{w}_{\text{black}}(t)=100(1-\sin(6\pi t/T))$ (kg/s),
$\bm{w}_{\text{cyan}}(t)=130(1+0.5\sin(6\pi t/T))$ (kg/s),
$\mu_{\text{red}}=1.0240$, $\underline{\mu}_{\text{yellow}}=1.0029$, and
$\underline{\mu}_{\text{purple}}=1.0199$, where $T=80$ (hr). The other
boundary condition is $\bm{q}_{\text{green}}(t)=0$ (kg/s) (solid lines) and
$\bm{q}_{\text{green}}(t)=2$ (kg/s) (dots).

##  VII Monotonic Interface

Proposition 2 shows that the total pressure and density system of ODEs is not
monotone-ordered over the entire boundary condition parameter region
$\mathcal{D}=\\{\bm{p}_{s},\bm{d},\bm{\alpha}^{(m)},\bm{\beta}^{(m)}\\}$.
However, by Proposition 1 and the continuity of solutions with respect to
initial conditions and inputs [41], the non-isolated total pressure and
density system of ODEs is expected to be monotone-ordered over a certain
subregion of $\mathcal{D}$ that consists of concentration vectors that are
uniformly close to a constant concentration vector. Again by continuity,
monotonicity is also expected to hold for slow variations in concentration
with large amplitudes. This suggests that there may be a nontrivial monotonic
interface (MI) that partitions $\mathcal{D}$ into monotonic and non-monotonic
phase regions for each equivalent system variable. Moreover, the simulation
results from Section VI suggest that the MIs for each equivalent system
variable may be significantly different from one another. We focus on
partitioning the subregion of $\mathcal{D}$ that consists of only
concentration boundary parameters, because the concentration variable is the
only factor that leads to conditions in which monotonicity does not hold.

We analyze the MI numerically for a single pipeline with concentration and
pressure specified at the inlet of the pipeline (node 1) and with mass outflow
specified at the outlet (node 2). The parameters and boundary conditions that
do not change are pipeline length $\ell=50$ (km), diameter $D=0.5$ (m),
friction factor $\lambda=0.11$, and constant slack node pressure
$\bm{p}_{s}=7$ (MPa). We denote the concentration of hydrogen at the inlet
slack node by $\bm{\alpha}_{1}(t)=\bm{\alpha}_{1}^{(2)}(t)$ and specify it to
be

$\bm{\alpha}_{1}(t)=\alpha_{1}\left(1+\kappa\sin(2\pi\omega_{*}t)\right),$
(30)

where $\kappa$ is the amplitude factor of the sinusoid, $\omega_{*}$ is the
forcing frequency in cycles per hour, and $\alpha_{1}$ is the mean
concentration profile around which the sinusoid oscillates. Here, the
subscript is with respect to the node number. The sub-region of $\mathcal{D}$
that we consider consists of all pairs $(\omega_{*},\kappa)$ with
$0\leq\omega_{*}\leq 2$ and $0\leq\kappa\leq 1$. We note that the pipeline
system in Eq. (17) with the boundary conditions specified as above may be
written as an autonomous system by extending the state space by two dimensions
and writing $\bm{\alpha}_{1}$ as a state solution of the harmonic oscillator.
We note this extended formulation so that our subsequent results can be
interpreted in the context of autonomous dynamical systems theory. However,
the extension is not necessary for the analysis, so we omit the details.

Figure 8: Monotonic interfaces $(\omega_{*},\kappa_{*}(\omega_{*}))$ in the
$(\omega_{*},\kappa)$ plane for each of the equivalent system variables listed
in the legend. Two solutions corresponding to two ordered mass outflows and
the concentration pair $(\omega_{*},\kappa)$ will not cross one another if
$\kappa<\kappa_{*}(\omega_{*})$.

In addition to the boundary conditions and parameters specified above, we use
the parameter values $\sigma_{1}=377$ (m/s), $\sigma_{2}=2.8\sigma_{1}$, and
$\alpha_{1}=0.02$ in the analyses in this section. We now describe our process
of computing the MI. For each $(\omega_{*},\kappa)$ in Eq. (30), we compute
three solutions corresponding to three monotone-ordered constant mass outflows
$\bm{w}_{2}=\overline{\varphi}\pi(D/2)^{2}$ (kg/s), where
$\overline{\varphi}=120$, 140, and 160 (kg/m2s). The region in the
$(\omega_{*},\kappa)$ plane defined by $0\leq\omega_{*}\leq 2$ and
$0\leq\kappa\leq 1$ is discretized into a $21\times 41$ grid of discrete
pairs. For each discrete $\omega_{*}$, we compute the three solutions for each
discrete $\kappa$ starting from $\kappa=0$ and increasing $\kappa$ until we
achieve the lower bound $\kappa=\kappa_{*}(\omega_{*})$ with which at least
two of the three solutions exhibit crossings at some point in time. Thus, the
three solutions corresponding to the pair $(\omega_{*},\kappa)$ will not cross
if $\kappa<\kappa_{*}(\omega_{*})$ and at least two of the three will cross if
$\kappa\geq\kappa_{*}(\omega_{*})$. The MI is defined by the set of pairs
$(\omega_{*},\kappa_{*}(\omega_{*}))$ in the plane, and each of the equivalent
system variables has its own MI associated with it. The cubic spline
interpolated MI curves for several equivalent system variables are depicted in
Fig. 8. The region below the MI curve is called the monotone response region
(MRR). Fig. 8 shows that the MRRs for hydrogen density, natural gas density,
total density, energy, and pressure are nested increasing sets where the MRRs
for hydrogen density and pressure are the smallest and largest sets,
respectively.

For time-varying concentration profiles, Fig. 8 suggests that the pressure and
energy equivalent system should be used if monotonicity properties are
important to the formulation. This is the conclusion that we arrive at in
Section VI. Of the five examples from that section, the only examples that
consider similar sinusoidal forcing in concentration are those that correspond
to Figs. 3 and 4. Both of these examples introduce a combination of two
sinusoidal forcing frequencies to the network, $0.1$ (cyc/hr) and 0.033
(cyc/hr). Recall from Figs. 3 and 4 that only the pressure and energy
solutions did not exhibit crossings for these two examples. This observation
agrees with the MIs in Fig. 8, where the operating points
$(\omega_{*},\kappa)=(0.1,1)$ and $(\omega_{*},\kappa)=(0.033,1)$ are above
all of the MIs except for the pressure and energy MIs. It is important to
mention that the MIs computed here for a single pipeline do not necessarily
match the MIs that would result for the network topology from Section VI, even
if the slack node pressure, total outflow, and other boundary condition
parameters of the network were equal to those used for the single pipeline.

Observe that the MRRs for the pressure and energy variables are significantly
larger in area than the MRRs for density, as seen in Fig. 8. This is inherent
because of the manner in which the boundary conditions have been specified.
Recall that the boundary condition at the slack node of the pipeline presented
here has been specified to maintain constant pressure. Although the model has
been derived more generally, constant slack node pressure is a common
specification for pipeline simulation, so we have followed this convention.
However, if the pressure is constant and the concentration varies at the slack
node, then the density will be forced to vary at the slack node, because of
the equation of state of the mixture. It follows that the density will
oscillate if the concentration oscillates, and there will be some propagation
of density oscillations at speeds related to the mass outflows. Because we
specify different outflows (that are monotone-ordered) to compute the MIs, the
density waves corresponding to the different outflows will typically be out of
phase and the solutions will exhibit crossings. Therefore, we expect the
density variables to be the most sensitive to ordering properties under such
boundary conditions, and this is apparent in Fig. 8.

As $\omega_{*}$ increases from $\omega_{*}=0$ to $\omega_{*}=2$ (cyc/hr), the
MI curves qualitatively decrease from unity to a lower bound, flatten out, and
then increase. The fact that the amplitude factor generally increases along
the MI as $\omega_{*}$ increases beyond $\omega_{*}=0.75$ is a robustness
feature of monotonicity to high frequency uncertainty. This property appears
to be a consequence of wave attenuation in strongly dissipative gas pipeline
flow [42]. In particular, the gas pipeline demonstrates low-pass filtering
characteristics with which the amplitudes of high frequency travelling waves
tend to be significantly attenuated over short distances, and, therefore, the
likelihood of conditions in which monotonicity does not hold decreases as the
frequency of the high frequency oscillation increases. If the concentration of
hydrogen injected into the network contains a small variation of high
frequency uncertainty, then the MIs demonstrate that this uncertainty
typically will not cause an otherwise theoretically monotonic operation to
become non-monotonic.

##  VIII Periodic Interface

We demonstrate that non-periodic solutions can arise from sinusoidal forcing
in concentration. To numerically study periodicity and the breakdown thereof,
we must simulate the solutions over long time intervals that span hundreds of
hours. In addition, we will consider large and fast variations in
concentration. As we have mentioned at the end of Section VI, fast variations
require an extremely fine spatial discretization size for the finite volume
discretization method presented in Section III. The small spatial
discretization size creates a large ODE system, which is difficult to
implement over a long time interval on a digital computer. Therefore in our
study of periodicity, instead of using the finite volume method, we discretize
the pipeline at the (translated) nodes of Chebyshev polynomials for which
exponential convergence properties are obtained (e.g., see [43]). We briefly
outline the method in Appendix B. The analysis presented in this section is
performed in the single 50 (km) pipeline that was used in the previous section
to study the MI, with $D=0.5$ (m), $\lambda=0.11$, and constant slack node
pressure $\bm{p}_{s}=7$ (MPa). However, in addition to the other parameters
specified in that section, we now use $\sigma_{1}=338.38$ (m/s),
$\sigma_{2}=4\sigma_{1}$, and $\alpha_{1}=0.2$, but all of the other
parameters remain the same.

To introduce our analysis on periodicity, we show three examples in Figs. 9 to
11 that share the same boundary conditions with one another except for the
different frequencies $\omega_{*}$ and amplitude factors $\kappa$ of the
sinusoidal concentration profile in Eq. (30). The top of the three figures
depict the pressure solutions at the outlet of the pipeline for $t\in[rT,T]$
with $0.7\leq r\leq 0.95$, where $T=400$ (hr). The tail-ends of the solutions
are used so that initial transient responses do not affect the analysis of
periodic orbits. The bottom left-hand-sides of Figs. 9 to 11 show the phase
space diagrams of outlet density and outlet pressure during the later stages
of the simulations. We see that the solutions in Figs. 9 and 10 approach
periodic orbits and that the solution in Fig. 11 does not appear to do so.
However, even the two periodic responses in Figs. 9 to 10 have certain
properties that are not observed in homogeneous natural gas simulations [42].
Particularly, in Fig. 9, multiple local minima in the pressure appear for
every local minimum of the sinusoidal forcing over the time interval
$[0.95T,T]$. The additional local minima correspond to the inner loop of the
periodic orbit. The pressure in Fig. 10 has the same number of local minima as
the forcing sinusoid over the interval $[0.75T,T]$, but has twice the period.
These examples demonstrate that periodic solutions arising from heterogeneous
mixtures of gases may be irregular in the following sense. From the laws of
fluid dynamics, gas pressure should decrease with decreasing density under
constant temperature and volume. However, because of the oscillating gas
composition, the phase space diagram in Fig. 9 contains sequences of four
small time intervals during which density decreases while pressure increases,
and the phase space diagram in Fig. 10 contains two such time intervals. Note
that the phase space diagram of outlet pressure and density for flow of a
single ideal gas in the transient regime is simply a line with positive slope,
which the solution traverses.

Figure 9: Pipeline solution with boundary conditions
$\bm{w}_{2}(t)=75\pi(D/2)^{2}$ (kg/s), $\omega_{*}=0.25$, and $\kappa=1.0$.
The periodicity measure in Eq. (32) is $\mathcal{P}=0.30$. Figure 10: Pipeline
solution with boundary conditions $\bm{w}_{2}(t)=75\pi(D/2)^{2}$ (kg/s),
$\omega_{*}=0.1$, and $\kappa=0.98$. The periodicity measure in Eq. (32) is
$\mathcal{P}=0.62$. Figure 11: Pipeline solution with boundary conditions
$\bm{w}_{2}(t)=75\pi(D/2)^{2}$ (kg/s), $\omega_{*}=0.5$, and $\kappa=0.9$. The
periodicity measure in Eq. (32) is $\mathcal{P}=1.19$.

The frequency responses of the outlet pressures are depicted on the bottom
right-hand-sides of Figs. 9 to 11 using the discrete Fourier transform [44]
defined below in Eq. (31). The dominant frequency mode in the solution appears
at the forcing frequency $\omega_{n}=\omega_{*}$ in Figs. 9 to 11. The
generated harmonic modes in Fig. 9 appear at integer multiples of
$\omega_{*}$. This behavior is typical for homogeneous natural gas pipeline
flow [42]. The generated harmonic modes in Fig. 10 appear at half the values
of the integer multiples of $\omega_{*}$. This behavior indicates period-
doubling bifurcations [45] at the forcing frequency $\omega_{*}=0.1$ as the
amplitude factor $\kappa$ increases. The pressure dynamics in Fig. 11 appear
to be comprised by a continuous distribution of generated harmonic modes.
These observations inspire a quantitative measure of periodicity in terms of
the frequency response of the solution. This is the approach taken in [26] for
the transition to chaotic responses in oceanic wind bursts. We define a
sequence of evenly-spaced samples of the tail-end of the outlet pressure by
$\bm{p}_{2}[k]=\bm{p}_{2}((0.6+k/N)T)$ for $k=0,\dots,0.4N$, where $N$ is
equal to the number of time samples of the numerical solution over the
interval $[0,T]$. For such a sampled sequence $\bm{\psi}[k]$, the normalized
discrete Fourier transform (DFT) is defined as

$\\{\mathcal{F}\bm{\psi}\\}[\omega_{n}]=\frac{\sum_{k=0}^{0.4N}\bm{\psi}[k]e^{-\bm{j}2\pi\omega_{n}k}}{\max_{\omega_{n}}\left|\sum_{k=0}^{0.4N}\bm{\psi}[k]e^{-\bm{j}2\pi\omega_{n}k}\right|},$
(31)

where $\bm{j}$ is the imaginary unit and $\omega_{n}=n/(0.4T)$ (cyc/hr) are
the sampling frequencies for $n=0,\dots,0.4N$. Periodicity is measured with
the average power spectrum defined by

$\displaystyle\mathcal{P}=\frac{1}{0.4N+1}\sum_{n=0}^{0.4N}\left|\\{\mathcal{F}(\bm{p}_{2}-\bm{\pi}_{2})\\}[\omega_{n}]\right|^{2}\times
100,$ (32)

where $\bm{\pi}_{2}=\bm{p}_{2}(0)$ is the initial steady-state value of
pressure at the outlet of the pipeline. The shifted pressure in the power
spectrum is used to suppress the zero frequency component of the initial
state.

The power spectrum $\mathcal{P}$ is depicted in the form of a color map as a
function of $(\omega_{*},\kappa)$ in Fig. 12, where $\omega_{*}$ is the
forcing frequency and $\kappa$ is its amplitude factor given in Eq. (30). This
figure has been obtained numerically as follows. Similarly to the way that we
have computed the MIs, the region in the $(\omega_{*},\kappa)$ plane defined
by $0\leq\omega_{*}\leq 2$ and $0.5\leq\kappa\leq 1$ is discretized into a
$21\times 25$ grid of discrete pairs. For each frequency and amplitude factor
of the forcing concentration on this grid, we numerically simulate the
solution in the pipeline for 400 hours. We then compute the normalized DFT and
the average power spectrum of the tail-end of the sampled solution as defined
above. These computations provide the discrete set of quantified values
depicted in Fig. 12. The periodic interface (PI) in Fig. 12 partitions the
$(\omega_{*},\kappa)$ plane into periodic and non-periodic response
subregions. The computation of the PI is performed as follows. For each
$\omega_{*}$, the parameter $\kappa$ is increased from $\kappa=0$ to
$\kappa=\kappa^{*}(\omega_{*})$, where $\kappa^{*}(\omega_{*})$ is the upper
bound on $\kappa$ below which the DFT of the outlet pressure consists of
countably many pulses. From numerical simulations, the average power spectrum
corresponding to the forcing pair $(\omega_{*},\kappa^{*}(\omega_{*}))$ is
typically around $\mathcal{P}=0.3$ (which depends on parameters and the number
of time samples of the numerical solution). Therefore, we quantitatively
define $\kappa^{*}(\omega_{*})$ to be the upper bound on $\kappa$ below which
$\mathcal{P}<0.3$. The periodic response region (PRR) is defined to be the set
of boundary condition parameter pairs $(\omega_{*},\kappa)$ with
$\kappa<\kappa^{*}(\omega_{*})$.

Figure 12: Color map of the power spectrum $\mathcal{P}$ in Eq. (32) as a
function of $(\omega_{*},\kappa)$ in Eq. (30). The boundary conditions are
$\alpha_{1}=0.2$ and $\bm{w}_{2}(t)=75$. In this figure, we plot the minimum
between 1 and $\mathcal{P}$ in Eq. (32).

##  IX Chaotic Interface

Figure 13: Example of the difference $\Psi[n]=\Psi(t_{n})$ between two outlet
pressure solutions and corresponding outlet density solutions in the pipeline
from Section VIII. The boundary condition parameters are
$(\omega_{*},\kappa)=(0.5,0.95)$, $\alpha_{1}=0.2$, and the initial conditions
of the two solutions correspond to steady mass flows $\bm{w}_{2}(0)=75(\pi
D^{2}/4)$ and $\bm{w}_{2}(0)=75.0001(\pi D^{2}/4)$.

The analysis presented in the previous section demonstrates that heterogeneous
mixtures of gases may behave irregularly for periodic boundary conditions with
frequency and amplitude $(\omega_{*},\kappa)$ outside of the periodic response
region. In this section, we demonstrate that the solutions that exhibit such
behavior are also chaotic in the sense of being highly sensitive to initial
conditions [46]. The extent of chaos in a finite-dimensional system can be
quantified by the largest Lyapunov exponent of the system [47]. This measure
provides an estimate on the exponential rate of divergence between two
solutions that begin their trajectories close to one another. Several
approaches have been developed to approximate the Lyapunov exponents of an
unknown finite-dimensional system from numerical data by using the time series
associated with a numerical solution of the system [48, 49, 50]. These methods
require an appropriate embedding dimension on which the calculation of the
Lyapunov exponent depends. Although there are methods that approximate the
optimal embedding dimension of time-series observations of a finite-
dimensional system [51], there may not always be a consistent approximation
for an observable generated by an underlying infinite-dimensional PDE system.
We measure the extent of chaos in our system by using a rate of divergence
between two specific solutions. It is important to note that the measure we
present does not necessarily provide an estimate of the largest Lyapunov
exponent of our system.

Our analysis will be performed in the single pipeline with parameters and
boundary conditions defined in the previous section. Consider two scalar-
valued time series solutions $\psi_{1}[n]$ and $\psi_{2}[n]$ with
$|\psi_{2}[0]-\psi_{1}[0]|<\delta$, where $\delta$ is small relative to
$\psi_{1}[0]$. Here, $\psi[n]=\psi(t_{n})$ represents any one of the dynamic
variables evaluated at the outlet of the pipe at discrete times $t_{n}=(n/N)T$
for $n=0,\dots,N$, where $(N+1)$ is the number of time samples of the
numerical solution over the time interval $[0,T]$. The rate of divergence
between the two solutions is defined to be the slope of a linear estimation of
the difference function

$\Psi[n]=\log\left|\frac{\psi_{2}[n]-\psi_{1}[n]}{\psi_{2}[0]-\psi_{1}[0]}\right|$

over an interval $I_{0}=[n_{0},n_{1}]$ during which the two solutions diverge
exponentially (assuming that they diverge). Fig. 13 shows an example of how
the difference function may increase with time for two divergent solutions
that began their trajectories close to one another. Of course, boundedness of
the physical system limits the divergence to within a bounded region in the
state space. For the example in Fig. 13, the rate of divergence may be
computed as the slope of a linear estimation of $\Psi[n]$ over the interval
$[0.5N,0.55N]$ corresponding to $t_{n}\in[100,110]$. One such linear
estimation is the least squares regression line.

The interval $I_{0}$ typically varies with respect to parameters and could be
difficult to construct for two general solutions that may or may not diverge.
For two solutions that do not diverge, the slope of the linear estimation of
$\Psi[n]$ on the interval $[0,N]$ would be at most approximately zero.
Moreover, there are usually several subintervals $I_{0}\subset[0,N]$ over
which the slope of the linear estimation of $\Psi$ changes sign, particularly
for two solutions that do not diverge. Consequently, the linear estimation of
$\Psi[n]$ over an interval $I_{0}$ may not always reflect the actual mean
growth of $|\psi_{2}[n]-\psi_{1}[n]|$ over the entire interval $[0,N]$. In an
attempt to moderate these difficulties, the amount of divergence between two
solutions $\psi_{1}$ and $\psi_{2}$ is defined by the measure

$\displaystyle\mathcal{C}:$ $\displaystyle=$
$\displaystyle\frac{1}{n_{2}-n_{1}}\left(\frac{1}{|I_{T}|}\sum_{n\in
I_{T}}\Psi[n]-\frac{1}{|I_{0}|}\sum_{n\in I_{0}}\Psi[n]\right)$
$\displaystyle=$ $\displaystyle\frac{1}{(n_{2}-n_{1})|I_{T}|}\sum_{n\in
I_{T}}\log|\psi_{2}[n]-\psi_{1}[n]|$
$\displaystyle\qquad-\frac{1}{(n_{2}-n_{1})|I_{0}|}\sum_{n\in
I_{0}}\log|\psi_{2}[n]-\psi_{1}[n]|,\qquad$

with $I_{T}=[n_{2},n_{3}]$ where $n_{1}<n_{2}$. The function $\mathcal{C}$ in
Eq. (IX) will be referred to as the chaos measure. This measure is interpreted
as the slope of the local mean of $\Psi$ or as the slope of the local mean
logarithmic difference of the two solutions computed over an initial interval
$I_{0}$ and a final interval $I_{T}$. Large and positive $\mathcal{C}$
indicates exponential divergence between the two time series (over some
interval). The means over the initial and final time intervals are used to
suppress fluctuations and estimate the expected values of the associated
quantities over the intervals. The interval $I_{0}$ must span a range that
precedes any exponential divergence and should also not include initial
transients. Moreover, the interval $I_{T}$ must follow any exponential
divergence.

Figure 14: Color map of the measure $\mathcal{C}$ in Eq. (IX) as a function of
$(\omega_{*},\kappa)$ in Eq. (30). The boundary conditions are the same as
those in Fig. 12 and the initial conditions of the two solutions correspond to
$\bm{w}_{2}(0)=75(\pi D^{2}/4)$ and $\bm{w}_{2}(0)=75.1(\pi D^{2}/4)$.

We analyze the measure $\mathcal{C}$ for the outlet pressure
$\psi[n]=p(t_{n},\ell)$ of the single pipeline used previously in Section VIII
with $\alpha_{1}=0.2$, but we note that equivalent system variables share
similar divergence properties. As for the computations of the power spectrum
shown in Fig. 12, we discretize the region in the $(\omega_{*},\kappa)$ plane
defined by $0\leq\omega_{*}\leq 2$ and $0.5\leq\kappa\leq 1$ into a $21\times
25$ grid of discrete pairs. For each frequency and amplitude factor of the
forcing concentration on this grid, we numerically simulate two solutions in
the pipeline for 100 hours whose initial conditions correspond to withdrawal
rates $\bm{w}_{2}(0)=75(\pi D^{2}/4)$ and $\bm{w}_{2}(0)=75.1(\pi D^{2}/4)$,
where the subscript refers to the outlet node. We then compute the value
$\mathcal{C}$ for each discrete pair on the grid using the intervals
$I_{0}=[0.08N,0.15N]$ and $I_{T}=[0.5N,0.8N]$ with $N=10,000$ and then depict
the collected values as a color map in Fig. 14. While we suppose that the
intervals $I_{0}$ and $I_{T}$ are sufficient for the specified initial
conditions, these intervals may not be adequate for the example in Fig. 13
because the solutions in that example are initially closer to one another and
require a larger time to begin diverging. For each frequency, the values of
$\mathcal{C}$ in Fig. 14 are more scattered than the associated values of
$\mathcal{P}$ in Fig. 12. In the latter case, the values of the power spectrum
either increase or decrease almost monotonically for each frequency as the
values of the amplitude factor increase. However, for a fixed frequency, the
values of the chaos measure may oscillate sporadically around $\mathcal{C}=0$
as $\kappa$ increases from $\kappa=0$ to $\kappa=1$.

Figure 15: Values of $\mathcal{C}$ as a function of $\kappa$ for the fixed
frequencies $\omega_{*}=0.5$ and $\omega_{*}=0.7$.

The chaotic interface (CI) is defined to be the set of pairs
$(\omega_{*},\kappa^{*})$, where $\kappa^{*}(\omega_{*})$ is the lower bound
that satisfies $\mathcal{C}>0$ for all $\kappa>\kappa^{*}(\omega_{*})$. If no
such value exists, we define $\kappa^{*}(\omega_{*})=1$. We say that the
solution corresponding to $(\omega_{*},\kappa)$ is chaotic if
$\kappa>\kappa^{*}(\omega_{*})$. The value of $\kappa^{*}(\omega_{*})$ is
determined numerically from the simulations that produced the chaos measure in
Fig. 14. The methodology of computing $\kappa^{*}(\omega_{*})$ is depicted in
Fig. 15. In this figure, $\mathcal{C}$ is plotted as a function of discrete
amplitude factors ranging between $0\leq\kappa\leq 1$ for two fixed
frequencies. The horizontal and vertical lines described by $\mathcal{C}=0$
and $\kappa=\kappa^{*}(\omega_{*})$ are depicted for reference. It should be
noted that small and positive values of $\mathcal{C}$ do not necessarily imply
that the associated time series diverge. This is to be expected because the
flow variables need not oscillate sinusoidally and may favor more time near
the peak of the wave.

Figure 16: Monotonic, periodic, and chaotic interfaces of the outlet pressure
for the 50 (km) pipeline used in Fig. 12 with $\alpha_{1}=0.2$.

We conclude our analysis with a comparison of the interfaces in the regions of
boundary condition parameters $(\omega_{*},\kappa)$ that do or do not exhibit
monotonic, periodic, and chaotic responses to periodic boundary conditions for
a single pipeline. These interfaces, which we refer to as the monotonic
interface (MI), periodic interface (PI), and chaotic interface (CI), are
depicted in Fig. 16 for the single pipeline system used in Sections VII and
VIII. The MI for pressure in Fig. 16 is different from the MI for pressure in
Fig. 8 because of the different mean concentrations $\alpha_{1}$ used to
compute the two MIs. Moreover, the withdrawal rates of the three solutions
that are used here for the computation of the current MI are
$\bm{w}_{2}=40\pi(D/2)^{2}$, $\bm{w}_{2}=75\pi(D/2)^{2}$ and
$\bm{w}_{2}=110\pi(D/2)^{2}$, which are different from those used in Section
VII. There are a few key takeaways from Fig. 16. First, the region beneath the
MI is a subset of the region beneath both the PI and CI and is significantly
beneath these two interfaces for nonzero frequencies. This suggests that all
sinusoidal hydrogen blending conditions that show monotonicity properties will
not result in irregular or chaotic responses, as expected for the strongly
dissipative behavior of natural gas pipeline flows. This result is informative
for natural gas pipeline operations with significant hydrogen blending.
Second, the PI and CI are qualitatively similar to one another. Strictly
speaking, the PI is beneath the CI for all $\omega_{*}\leq 1.1$ and the CI is
beneath the PI for all $\omega_{*}>1.1$. We note that the interfaces depend on
the computational method and its parameter values, the discretization grid in
the space of boundary condition parameters, and the threshold chosen to
delimit the interfaces (i.e., $\mathcal{P}=0.3$ for the PI and
$\mathcal{C}=0.0$ for the CI). There may be other reasonable thresholds that
yield somewhat different quantitative results, such as one where the PI is
always beneath the CI or one where the PI is exactly equal to the CI at all of
the grid points. The qualitative similarity of the PI and CI suggests that
these regions may indeed coincide, and only one of these measures may be
sufficient to identify both non-periodic and chaotic behavior.

##  X Conclusions

We have developed a model for transporting heterogeneous mixtures of natural
gas and hydrogen through pipeline networks. The formulation may be applied to
real pipeline systems with time-varying operations of compressor and regulator
units, supply stations that inject gas into the network at defined pressure
and hydrogen fraction, and flow stations that withdraw the mixture from the
network. The nonlinear partial differential equation formulation is
discretized using a finite volume method to obtain a nonlinear input-to-state
system, for which we prove monotone ordering properties for injections with
constant hydrogen fraction, and prove that such monotonicity properties do not
hold in general for injections with time-varying hydrogen concentrations. This
result builds on previous work that assumed globally homogeneous gas [8], by
considering inhomogeneity in space (Proposition 1 on injections with constant
hydrogen fraction) and time (Proposition 2 on injections with time-varying
hydrogen fraction). The interface in the boundary condition parameter region
of concentration variation that partitions monotonic and non-monotonic
responses was analyzed numerically and the results were illustrated on a test
network. Operations outside of the monotone response region may create surges
with large pressure, energy, and concentration gradients, which do not occur
in flows of a homogeneous gas. The monotonic interface analysis indicates that
sufficiently slow variation in concentration about a constant profile will
likely maintain monotonicity of ordered solutions in overall system pressures,
and prevent large, rapid pressure transients. Such conditions are critical to
maintain a physical flow regime with behavior that is intuitive for pipeline
control room operators. This suggests that hydrogen may be blended into a
natural gas pipeline network as long as injection rates are changed only
gradually. The acceptable ramping rates depend significantly on the structure
of the network, and would have to be determined through numerous simulations.

This study demonstrates that heterogeneous location and time dependent
boundary conditions may result in non-periodic and chaotic flows when hydrogen
is blended into a natural gas pipeline. For a single pipeline, boundary
conditions with monotonic, periodic, and chaotic responses were analyzed
numerically and interfaces between regions where these properties do and do
not hold were estimated. The interface analysis of a single pipeline
demonstrates that sinusoidal boundary conditions for which monotonicity
properties are preserved lead to solutions that are not chaotic and will
eventually approach a periodic response. We have also demonstrated that the
interfaces that delimit periodic and chaotic responses to periodic boundary
conditions are qualitatively similar and may indeed coincide for appropriate
thresholds that are used to define these interfaces. This suggests that only
one of these interfaces may be sufficient to identify both non-periodic and
chaotic behavior. Characterizing the monotonic, periodic, and chaotic
interfaces with extensive simulations for specific network topologies may
enable a gas pipeline system designer to determine limitations on operating
their networks safely and predictably given blending of heterogeneous gases.

The results of this study enable new capabilities to plan for and operate the
transmission of gases through pipeline energy systems. The generalizable
modeling presented here can be used in addition to empirical studies to
quantify the effect of hydrogen blending on gas pipelines under transient
conditions [52]. The developed model can be applied in general natural gas
operations if two or more production sites supply notably different levels of
methane in their natural gas compositions. More generally, the derivation of
the gas mixture model may be extended from a mixture of two gases to any
finite number of distinct gases. We aim to present a self-contained analysis
of a variety of different physical properties of gas mixture flows in pipeline
networks. Much of our analysis was performed on a finite-dimensional
realization of the infinite-dimensional PDE system for a small network and a
single pipeline. Therefore, there are still a number of questions that remain
open in terms of monotonicity, periodicity, and chaoticity. Rigorous
definitions and theorems regarding the limits of monotonicity, periodicity,
and chaoticity of nonlinear PDEs on networks and their dependencies on graph
topology have yet to be examined. Such results may be used to establish a more
complete analysis of appropriate blending conditions with which monotone-
ordering and well-behaved flows are maintained. The onset of turbulence in
complex multi-dimensional fluid mixing flows is extensively studied [53, 54],
and here we show that complex, chaotic behavior can arise in a strongly
dissipative, essentially one-dimensional distributed system.

## Acknowledgements

The authors are grateful to Vitaliy Gyrya, Rodrigo Platte, Dieter Armbruster,
and Yan Brodskyi for numerous helpful discussions, and to E. Olga Skowronek
for drawing the networks in Figs. 1 and 2. This study was supported by the
U.S. Department of Energy’s Advanced Grid Modeling (AGM) project “Dynamical
Modeling, Estimation, and Optimal Control of Electrical Grid-Natural Gas
Transmission Systems”, as well as LANL Laboratory Directed R&D project
“Efficient Multi-scale Modeling of Clean Hydrogen Blending in Large Natural
Gas Pipelines to Reduce Carbon Emissions”. Research conducted at Los Alamos
National Laboratory is done under the auspices of the National Nuclear
Security Administration of the U.S. Department of Energy under Contract No.
89233218CNA000001.

## References

  * Nalley and LaRose [2022] Stephen Nalley and Angelina LaRose.  Annual energy outlook 2022 (AEO2022).  _U.S. Energy Information Agency_ , 2022.
  * Salvi and Subramanian [2015] B. L. Salvi and K. A. Subramanian.  Sustainable development of road transportation sector using hydrogen energy system.  _Renewable and Sustainable Energy Reviews_ , 51:1132–1155, 2015.
  * Götz et al. [2016] Manuel Götz, Jonathan Lefebvre, Friedemann Mörs, Amy McDaniel Koch, Frank Graf, Siegfried Bajohr, Rainer Reimert, and Thomas Kolb.  Renewable power-to-gas: A technological and economic review.  _Renewable energy_ , 85:1371–1390, 2016.
  * Ozturk and Dincer [2021] Merve Ozturk and Ibrahim Dincer.  A comprehensive review on power-to-gas with hydrogen options for cleaner applications.  _International Journal of Hydrogen Energy_ , 46(62):31511–31522, 2021.
  * van der Hoeven [2004] Tom van der Hoeven.  _Math in gas and the art of linearization_.  Energy Delta Institute Groningen, The Netherlands, 2004.
  * Melaina et al. [2013] Marc W. Melaina, Olga Antonia, and Michael Penev.  Blending hydrogen into natural gas pipeline networks: A review of key issues.  _National Renewable Energy Laboratory_ , (TP-5600-51995), 2013.
  * Witkowski et al. [2017] Andrzej Witkowski, Andrzej Rusin, Mirosław Majkut, and Katarzyna Stolecka.  Comprehensive analysis of hydrogen compression and pipeline transportation from thermodynamics and safety aspects.  _Energy_ , 141:2508–2518, 2017.  ISSN 0360-5442.
  * Misra et al. [2020] Sidhant Misra, Marc Vuffray, and Anatoly Zlotnik.  Monotonicity properties of physical network flows and application to robust optimal allocation.  _Proceedings of the IEEE_ , 108(9):1558–1579, 2020.
  * Uilhoorn [2009] F. E. Uilhoorn.  Dynamic behaviour of non-isothermal compressible natural gases mixed with hydrogen in pipelines.  _International journal of hydrogen energy_ , 34(16):6722–6729, 2009.
  * Chaczykowski et al. [2018] Maciej Chaczykowski, Filip Sund, Paweł Zarodkiewicz, and Sigmund Mongstad Hope.  Gas composition tracking in transient pipeline flow.  _Journal of Natural Gas Science and Engineering_ , 55:321–330, 2018.
  * Guandalini et al. [2017] Giulio Guandalini, Paolo Colbertaldo, and Stefano Campanari.  Dynamic modeling of natural gas quality within transport pipelines in presence of hydrogen injections.  _Applied energy_ , 185:1712–1723, 2017.
  * Elaoud et al. [2017] Sami Elaoud, Zahreddine Hafsi, and Lamjed Hadj-Taieb.  Numerical modelling of hydrogen-natural gas mixtures flows in looped networks.  _Journal of Petroleum Science and Engineering_ , 159:532–541, 2017.
  * Hafsi et al. [2019] Zahreddine Hafsi, Sami Elaoud, and Manoranjan Mishra.  A computational modelling of natural gas flow in looped network: Effect of upstream hydrogen injection on the structural integrity of gas pipelines.  _Journal of Natural Gas Science and Engineering_ , 64:107–117, 2019.
  * Fan et al. [2021] Di Fan, Jing Gong, Shengnan Zhang, Guoyun Shi, Qi Kang, Yaqi Xiao, and Changchun Wu.  A transient composition tracking method for natural gas pipe networks.  _Energy_ , 215:119131, 2021.
  * Agaie et al. [2017] Baba G. Agaie, Ilyas Khan, Ali Saleh Alshomrani, and Aisha M. Alqahtani.  Reduced-order modelling for high-pressure transient flow of hydrogen-natural gas mixture.  _The European Physical Journal Plus_ , 132(5):1–16, 2017.
  * Subani et al. [2017] Norazlina Subani, Norsarahaida Amin, and Baba Galadima Agaie.  Leak detection of non-isothermal transient flow of hydrogen-natural gas mixture.  _Journal of Loss Prevention in the Process Industries_ , 48:244–253, 2017.  ISSN 0950-4230.
  * Zlotnik et al. [2015] Anatoly Zlotnik, Michael Chertkov, and Scott Backhaus.  Optimal control of transient flow in natural gas networks.  In _54th IEEE Conference on Decision and Control (CDC)_ , pages 4563–4570. IEEE, 2015.
  * Grundel et al. [2013] Sara Grundel, Nils Hornung, Bernhard Klaassen, Peter Benner, and Tanja Clees.  Computing surrogates for gas network simulation using model order reduction.  In _Surrogate-Based Modeling and Optimization_ , pages 189–212. Springer, 2013.
  * Grundel et al. [2014] Sara Grundel, Lennart Jansen, Nils Hornung, Tanja Clees, Caren Tischendorf, and Peter Benner.  Model order reduction of differential algebraic equations arising from the simulation of gas transport networks.  In _Progress in Differential-Algebraic Equations: Deskriptor 2013_ , pages 183–205. Springer, 2014.
  * Wong and Larson [1968] P. Wong and R. Larson.  Optimization of natural-gas pipeline systems via dynamic programming.  _IEEE Transactions on Automatic Control_ , 13(5):475–481, 1968.
  * Percell and Ryan [1987] Peter B. Percell and Michael J. Ryan.  Steady state optimization of gas pipeline network operation.  In _PSIG annual meeting_. OnePetro, 1987.
  * Rachford and Carter [2000] Henry H. Rachford and Richard G. Carter.  Optimizing pipeline control in transient gas flow.  In _PSIG annual meeting_. OnePetro, 2000.
  * Zlotnik et al. [2019] Anatoly Zlotnik, Kaarthik Sundar, Aleksandr M Rudkevich, Aleksandr Beylin, and Xindi Li.  Optimal control for scheduling and pricing intra-day natural gas transport on pipeline networks.  In _58th Conference on Decision and Control_ , pages 4887–4884. IEEE, 2019.
  * Himpe et al. [2021] Christian Himpe, Sara Grundel, and Peter Benner.  Model order reduction for gas and energy networks.  _Journal of Mathematics in Industry_ , 11(1):1–46, 2021.
  * Zlotnik et al. [2016] Anatoly Zlotnik, Sidhant Misra, Marc Vuffray, and Michael Chertkov.  Monotonicity of actuated flows on dissipative transport networks.  In _2016 European Control Conference (ECC)_ , pages 831–836. IEEE, 2016.
  * Tziperman et al. [1994] Eli Tziperman, Lewi Stone, Mark A. Cane, and Hans Jarosh.  El Niño chaos: Overlapping of resonances between the seasonal cycle and the pacific ocean-atmosphere oscillator.  _Science_ , 264(5155):72–74, 1994.
  * Eisenman et al. [2005] Ian Eisenman, Lisan Yu, and Eli Tziperman.  Westerly wind bursts: Enso’s tail rather than the dog?  _Journal of Climate_ , 18(24):5224–5238, 2005\.
  * Pizza et al. [2008] Gianmarco Pizza, Christos E. Frouzakis, John Mantzaras, Ananias G. Tomboulides, and Konstantinos Boulouchos.  Dynamics of premixed hydrogen/air flames in microchannels.  _Combustion and Flame_ , 152(3):433–450, 2008\.
  * Alipoor and Mazaheri [2016] Alireza Alipoor and Kiumars Mazaheri.  Combustion characteristics and flame bifurcation in repetitive extinction-ignition dynamics for premixed hydrogen-air combustion in a heated micro channel.  _Energy_ , 109:650–663, 2016.
  * Gyrya and Zlotnik [2019] Vitaliy Gyrya and Anatoly Zlotnik.  An explicit staggered-grid method for numerical simulation of large-scale natural gas pipeline networks.  _Applied Mathematical Modelling_ , 65:34–51, 2019.
  * Osiadacz [1984] Andrzej Osiadacz.  Simulation of transient gas flows in networks.  _International Journal for Numerical Methods in Fluids_ , 4(1):13–24, 1984.
  * Herty et al. [2010] Michael Herty, Jan Mohring, and Veronika Sachers.  A new model for gas flow in pipe networks.  _Mathematical Methods in the Applied Sciences_ , 33(7):845–855, 2010.
  * Sundar and Zlotnik [2019] Kaarthik Sundar and Anatoly Zlotnik.  State and parameter estimation for natural gas pipeline networks using transient state data.  _IEEE Transactions on Control Systems Technology_ , 27(5):2110–2124, 2019.
  * Silberberg et al. [2006] Martin Stuart Silberberg, Patricia Amateis, Rashmi Venkateswaran, and Lydia Chen.  _Chemistry: The molecular nature of matter and change_ , volume 4.  McGraw-Hill New York, 2006.
  * Kunz and Wagner [2012] Oliver Kunz and Wolfgang Wagner.  The gerg-2008 wide-range equation of state for natural gases and other mixtures: an expansion of gerg-2004.  _Journal of chemical & engineering data_, 57(11):3032–3091, 2012.
  * Hante and Schmidt [2019] Falk M. Hante and Martin Schmidt.  Complementarity-based nonlinear programming techniques for optimal mixing in gas networks.  _EURO journal on computational optimization_ , 7(3):299–323, 2019.
  * Gugat et al. [2012] Martin Gugat, Michael Herty, Axel Klar, Günter Leugering, and Veronika Schleper.  Well-posedness of networked hyperbolic systems of balance laws.  In _Constrained optimization and optimal control for partial differential equations_ , volume 160, pages 123–146. Springer, 2012.
  * Angeli and Sontag [2003] David Angeli and Eduardo D. Sontag.  Monotone control systems.  _IEEE Transactions on automatic control_ , 48(10):1684–1698, 2003.
  * Hirsch and Smith [2006] Morris W. Hirsch and Hal Smith.  Monotone dynamical systems.  _Handbook of differential equations: ordinary differential equations_ , 2:239–357, 2006.
  * Baker et al. [2023] Luke S. Baker, Saif R. Kazi, Rodrigo B. Platte, and Anatoly Zlotnik.  Optimal control of transient flows in pipeline networks with heterogeneous mixtures of hydrogen and natural gas.  In _2023 American Control Conference (ACC)_ , pages 1221–1228. IEEE, 2023.
  * Khalil [2002] Hassan K. Khalil.  _Nonlinear Systems_.  Pearson Education. Prentice Hall, 2002.
  * Baker et al. [2021] Luke Baker, Dieter Armbruster, Anna Scaglione, and Rodrigo B. Platte.  Analysis of a model of a natural gas pipeline—a transfer function approach.  _Transactions of Mathematics and Its Applications_ , 5(1):tnab002, 2021.
  * Ascher and Greif [2011] Uri M. Ascher and Chen Greif.  _A first course on numerical methods_.  SIAM, 2011.
  * Oppenheim et al. [2001] Alan V. Oppenheim, John R. Buck, and Ronald W. Schafer.  _Discrete-time signal processing. Vol. 2_.  Upper Saddle River, NJ: Prentice Hall, 2001.
  * Zhang et al. [2009] Jie Zhang, Nan-Sheng Liu, and Xi-Yun Lu.  Route to a chaotic state in fluid flow past an inclined flat plate.  _Physical Review E_ , 79(4):045306(R), 2009.
  * Lorenz [1995] Edward N. Lorenz.  _The essence of chaos_.  University of Washington press, 1995.
  * Benettin et al. [1980] Giancarlo Benettin, Luigi Galgani, Antonio Giorgilli, and Jean-Marie Strelcyn.  Lyapunov characteristic exponents for smooth dynamical systems and for hamiltonian systems; a method for computing all of them. part 1: Theory.  _Meccanica_ , 15:9–20, 1980.
  * Wolf et al. [1985] Alan Wolf, Jack B. Swift, Harry L. Swinney, and John A. Vastano.  Determining lyapunov exponents from a time series.  _Physica D: nonlinear phenomena_ , 16(3):285–317, 1985.
  * Rosenstein et al. [1993] Michael T. Rosenstein, James J. Collins, and Carlo J. De Luca.  A practical method for calculating largest lyapunov exponents from small data sets.  _Physica D: Nonlinear Phenomena_ , 65(1-2):117–134, 1993.
  * Brown et al. [1991] Reggie Brown, Paul Bryant, and Henry D. I. Abarbanel.  Computing the lyapunov spectrum of a dynamical system from an observed time series.  _Physical review A_ , 43(6):2787, 1991.
  * Kennel et al. [1992] Matthew B. Kennel, Reggie Brown, and Henry D. I. Abarbanel.  Determining embedding dimension for phase-space reconstruction using a geometrical construction.  _Physical review A_ , 45(6):3403, 1992.
  * Stỳblo et al. [2022] Martin Stỳblo, Vít Meistr, and Marko Hauliš.  Effects of hydrogen admixtures on pipeline capacity.  In _PSIG Annual Meeting_ , pages PSIG–2217. PSIG, 2022.
  * Amon et al. [1996] Cristina H. Amon, Amador M. Guzmán, and Benoit Morel.  Lagrangian chaos, eulerian chaos, and mixing enhancement in converging–diverging channel flows.  _Physics of Fluids_ , 8(5):1192–1206, 1996.
  * Feigenbaum [1980] Mitchell J. Feigenbaum.  The transition to aperiodic behavior in turbulent systems.  _Communications in mathematical physics_ , 77(1):65–86, 1980.

##  Appendix A Model Comparison

The interested reader is referred to [30] for a model comparison of
homogeneous gas flow in the network shown in Fig. 2. Using the same initial
and boundary conditions as used in the previous study, we recover the same
solution, up to machine precision, with our mixed gas model. Our solution is
shown in Fig. 17.

Figure 17: Model comparison (see Fig. 11 in [30]).

##  Appendix B Chebyshev Spectral Differentiation

Consider a single pipeline of length $\ell$, diameter $D$, and friction factor
$\lambda$ with axial variable $x\in[0,\ell]$. Discretize the interval
$[0,\ell]$ with the $(N+1)$ discretization points
$x_{i}=\ell/2(1-\cos(i\pi/N))$ for $i=0,\dots,N$. Define the sampled variables
$\bm{\rho}^{(m)}_{i}(t)=\rho^{(m)}(t,x_{i})$ and
$\bm{\varphi}_{i}(t)=\varphi(t,x_{i})$. It follows from interpolating the
values of $\bm{\rho}^{(m)}_{i}(t)$ at the points $x_{i}$ using Lagrange
polynomials of order $N$ that (e.g., see [43])

$\partial_{x}\rho^{(m)}(t,x_{i})\approx\bm{D}\bm{\rho}^{(m)}_{i}(t),$ (34)

where

$\displaystyle\bm{D}_{ij}=\begin{cases}\sum\limits_{\begin{subarray}{c}n=0\\\
n\neq j\end{subarray}}^{n}\frac{1}{x_{j}-x_{n}},&i=j,\\\
\frac{1}{x_{j}-x_{i}}\prod\limits_{\begin{subarray}{c}n=0\\\ n\neq
i,j\end{subarray}}^{n}\frac{x_{i}-x_{n}}{x_{j}-x_{n}}&i\not=j.\end{cases}$
(35)

The discretized PDEs in Eqs. (3)-(4) become

$\displaystyle\dot{\bm{\rho}}^{(m)}+\bm{D}\left(\frac{\bm{\rho}^{(m)}}{\bm{\rho}^{(1)}+\bm{\rho}^{(2)}}\odot\bm{\varphi}\right)$
$\displaystyle=$ $\displaystyle 0,$ (36)
$\displaystyle\bm{D}\left(\sigma^{2}_{1}\bm{\rho}^{(1)}+\sigma^{2}_{2}\bm{\rho}^{(2)}\right)$
$\displaystyle=$
$\displaystyle-\frac{\lambda}{2D}\frac{\bm{\varphi}\odot|\bm{\varphi}|}{\bm{\rho}^{(1)}+\bm{\rho}^{(2)}}.\quad$
(37)

The boundary conditions are incorporated into the discretized equations by
replacing $\bm{\rho}^{(m)}_{0}(t)=\bm{s}_{0}^{(m)}(t)$ and
$\bm{\varphi}_{N}(t)=\bm{w}_{N}(t)/(0.25\pi D^{2})$.

##  Appendix C Nomenclature

Network Graph

$\displaystyle\mathcal{E}$ Directed Edge Set $\displaystyle\mathcal{V}$ Node
Set $\displaystyle\mathcal{V}_{s}\subset\mathcal{V}$ Slack Nodes
$\displaystyle\mathcal{V}_{d}\subset\mathcal{V}$ Non-slack Nodes
$\displaystyle\mathcal{V}_{w}\subset\mathcal{V}_{d}$ Non-slack Withdrawal
Nodes $\displaystyle\mathcal{V}_{q}\subset\mathcal{V}_{d}$ Non-slack Injection
Nodes $\displaystyle E$ Cardinality of $\mathcal{E}$ $\displaystyle V$
Cardinality of $\mathcal{V}$ $\displaystyle V_{s}$ Cardinality of
$\mathcal{V}_{s}$ $\displaystyle V_{d}$ Cardinality of $\mathcal{V}_{d}$
$\displaystyle V_{w}$ Cardinality of $\mathcal{V}_{w}$ $\displaystyle V_{q}$
Cardinality of $\mathcal{V}_{q}$ $\displaystyle k$ Edge Index $\displaystyle
i\;\&\;j$ Node Indices $\displaystyle k:i\mapsto j$ Edge $k$ Directs from Node
$i$ to Node $j$ ${}_{\mapsto}j$ Set of Edges Directed to Node $j$
$\displaystyle j_{\mapsto}$ Set of Edges Directed from Node $j$

Edge Variables

$\displaystyle\rho$ Total Density (kg m-3) $\displaystyle p$ Total Pressure
(MPa) $\displaystyle\varphi$ Total Mass Flux (kg m-2 s-1)
$\displaystyle\rho^{(m)}$ Partial Density (kg m-3) $\displaystyle p^{(m)}$
Partial Pressure (MPa) $\displaystyle\varphi^{(m)}$ Partial Mass Flux (kg m-2
s-1) $\displaystyle\eta^{(m)}$ Mass Fraction $\displaystyle\nu^{(m)}$
Volumetric Fraction $\displaystyle\underline{\psi}$ Evaluation of $\psi$ at
Edge Inlet $\displaystyle\overline{\psi}$ Evaluation of $\psi$ at Edge Outlet

Node Variables

$\displaystyle\bm{\rho}$ Total Density (kg m-3) $\displaystyle\bm{p}$ Total
Pressure (MPa) $\displaystyle\bm{E}$ Energy (GJ)
$\displaystyle\bm{\rho}^{(m)}$ Partial Density (kg m-3)
$\displaystyle\bm{\eta}^{(m)}$ Mass Fraction $\displaystyle\bm{\nu}^{(m)}$
Volumetric Fraction

Boundary Condition Variables

$\displaystyle\bm{s}^{(m)}$ Slack Node Partial Density (kg m-3)
$\displaystyle\bm{p}_{s}$ Slack Node Pressure (MPa) $\displaystyle\bm{\alpha}$
Slack Node Mass Fraction $\displaystyle\bm{\beta}$ Non-slack Injection Node
Mass Fraction $\displaystyle\bm{a}$ Slack Node Wave Speed
$\displaystyle\bm{b}$ Non-slack Injection Node Wave Speed
$\displaystyle\bm{w}$ Non-slack Withdrawal Node Mass Outflow (kg s-1)
$\displaystyle\bm{q}$ Non-slack Injection Node Mass Inflow (kg s-1)

Control Variables

$\displaystyle\underline{\mu}$ Compressor Ratio $\displaystyle\overline{\mu}$
Regulator Ratio

Network Graph Matrices

$\displaystyle M$ Weighted Incidence Matrix $\displaystyle M_{s}$ Weighted
Slack Node Incidence Submatrix $\displaystyle M_{d}$ Weighted Non-slack Node
Incidence Submatrix $\displaystyle Q$ Incidence Matrix $\displaystyle Q_{s}$
Slack Node Incidence Submatrix $\displaystyle Q_{d}$ Non-slack Node Incidence
Submatrix $\displaystyle\underline{A}$ Componentwise Negative Parts of $A$
$\displaystyle\overline{A}$ Componentwise Positive Parts of $A$

Parameters

$\displaystyle\sigma_{1}$ Natural Gas Wave Speed $\displaystyle\sigma_{2}$
Hydrogen Wave Speed $\displaystyle\sigma$ Local Mixture Wave Speed
$\displaystyle\ell$ Pipe Length (km) $\displaystyle L$ Diag($\ell_{k}$)
$\displaystyle D$ Pipe Diameter (m) $\displaystyle\chi$ Pipe Cross-sectional
Area (m2) $\displaystyle X$ Diag($\chi_{k}$) $\displaystyle\lambda$ Darcy-
Weisbach Friction Factor $\displaystyle\Lambda$
Diag($\sqrt{2D_{k}/\lambda_{k}\ell_{k}}$)