# Energy levels of graphene magnetic quantum dot in inhomogeneous gap

Fatima Belokda  Laboratory of Theoretical Physics, Faculty of Sciences,
Chouaïb Doukkali University, PO Box 20, 24000 El Jadida, Morocco  Ahmed Jellal
<EMAIL_ADDRESS>Laboratory of Theoretical Physics, Faculty of Sciences,
Chouaïb Doukkali University, PO Box 20, 24000 El Jadida, Morocco  Canadian
Quantum Research Center, 204-3002 32 Ave Vernon, BC V1T 2L7, Canada  El
Houssine Atmani  Laboratory of Nanostructures and Advanced Materials,
Mechanics and Thermofluids, Faculty of Sciences and Techniques, Hassan II
University, Mohammedia, Morocco

###### Abstract

We investigate the energy levels of charge carriers confined in a magnetic
quantum dot in graphene with an inhomogeneous gap through an electrical
potential. We solve the eigenvalue equation for two regions. We explicitly
determine the eigenspinors for both valleys $K$, $K^{\prime}$ and use the
boundary condition at the quantum dot interface to obtain the energy levels.
We show that the energy levels exhibit symmetric and asymmetric behavior under
appropriate conditions of the physical parameters. It has been found that
changing the energy levels by introducing an energy gap outside the quantum
dot changes the electrical properties.

###### pacs:

73.22.Pr, 72.80.Vp, 73.63.-b
Keywords: Graphene, dot dots, energy gap, energy levels.

##  I Introduction

Graphene is defined as a two-dimensional layer composed only of carbon atoms
arranged in a honeycomb pattern [1]. The structure and the nature of the bonds
between the atoms that compose it give graphene extraordinary and unique
electronic properties [2, 3]. The dispersion relation is linear such that the
valence band and the conduction band touch at two points, K and K’, called
Dirac points. At these points, the graphene electrons behave as massless Dirac
fermions [3, 4] and are described by the Dirac equation [3]. In the framework
of band theory, graphene is a zero-gap semiconductor. In general, the ultra-
relativistic nature of graphene’s charge carriers has led researchers to
question how they would respond to confinement [5]. It is precisely this
particular property that prohibits the use of fabrication techniques. The
question of confining Dirac fermions in graphene has led to many proposals.
There are several theoretical methods for confining Dirac fermions in single-
layer graphene [6, 7], such as inhomogeneous magnetic fields [9, 8], potential
cylindrical symmetry [7], spatial modulation of the Dirac gap [10], cutting
the flake into small nanostructures [11, 12], using the substrate to induce a
band gap [13] and so on.

Graphene quantum dots (QDs) are small disk-like pieces of graphene with a
radius of $r_{0}$ in which the electronic wave function is confined and
exhibits quantum confinement effects, regardless of size. Therefore, graphene
QDs have a non-zero band gap and are luminescent by excitation. This band gap
is tunable by changing the size of the graphene QDs. Since its discovery,
scientists have been trying to trap electrons in graphene-based QDs in
anticipation of the wide field of novel applications of QDs in electronic
devices [15], valves [8], photovoltaics [16], qubits [17], and gas detection
[18]. QDs made from nanostructures are very specific to the exact shape of the
edge, which is difficult to control [8].

Motivated by our previous work [19], we study the influence of two
inhomogenous gaps on the energy spectrum of the graphene QDs. For this, we
consider the confinement of charge carriers in a magnetic graphene quantum dot
with the gap $\Delta_{1}$ inside and the gap $\Delta_{2}$ outside the dot. To
obtain the eigenspinors, we solve the Dirac equation separately in each region
of the system. By using the boundary condition at the interface of the quantum
dot, we obtain an equation describing the energy levels as a function of the
physical parameters. We numerically study the energy levels as a function of
the radius of QDs, the magnetic field, the electrostatic potential, and two
gaps $\Delta_{1}$, $\Delta_{2}$.

This paper is organized as follows. In Sec. II, we set up the mathematical
tools needed to study the present system. In Sec. III, we determine the
eigenspinors that describe the fermions in our theoretical model. Using the
boundary condition, we derive a formula governing the energy levels as a
function of the physical parameters in Sec. IV. We numerically analyze the
energy levels in Sec. V. Finally, we conclude our results.

##  II Model and theory

Let us consider the magnetic confinement of Dirac fermions in graphene
subjected to an electrostatic potential and two inhomogeneous gaps. We create
a quantum dot with radius $r_{0}$ in the presence of two different gaps, one
inside the other

$\Delta(r)=\left\\{\begin{array}[]{ll}\Delta_{1},&r>r_{0}\\\
\Delta_{2},&r<r_{0}\end{array}\right.$ (1)

using the magnetic field along $z$-direction shown below

$\vec{B}=\left\\{\begin{array}[]{ll}B\vec{e}_{z},&r<r_{0}\\\
0,&r>r_{0}.\end{array}\right.$ (2)

Such a system can be illustrated, as shown in Fig. 1.

(a)  Figure 1: (color online) A circular quantum dot of radius $r_{0}$
surrounded by a graphene sheet in the presence of a perpendicular magnetic
field $B$ inside and inhomogeneous gaps is depicted schematically.

The following Hamiltonian can be used to describe the dynamics of carriers in
the honycomb lattice of covalently bound carbon atoms in a single graphene

$H=v_{F}(\vec{p}+e\vec{A})\cdot\vec{\sigma}+V(r)\mathbb{I}_{2}+\Delta{(r)}\sigma_{z}$
(3)

where the Fermi velocity is $v_{F}\approx 10^{6}$ m/s, $\vec{p}=(p_{x},p_{y})$
is the two-dimensional momentum operator, $\vec{\sigma}$ denotes the Pauli
matrices, $\mathbb{I}_{2}$ is the $2$ $\times$ $2$ identity matrix, and $V$ is
the applied potential. The vector potential can be calculated using (2) to get

$\vec{A}=\left\\{\begin{array}[]{ll}\dfrac{B}{2r}\left(r^{2}-R^{2}\right)\vec{e_{z}},&r<r_{0}\\\
0,&r>r_{0}.\end{array}\right.$ (4)

In polar coordinates $(r,\theta)$, the Hamiltonian (3) reduces to the form

$H=\begin{pmatrix}V_{+}&\pi^{+}\\\ \pi^{-}&V_{-}\end{pmatrix}$ (5)

and the momentum operators are introduced

$\pi^{\pm}=-i\hbar\nu_{F}e^{\pm i\tau\theta}\left(\partial_{r}\pm\dfrac{\tau
i}{r}\partial_{\theta}\right)$ (6)

where we have defined

$V_{\pm}=V\pm\Delta_{1}\pm\Delta_{2}.$ (7)

We can now calculate the energy spectrum using the eigenvalues equation,
$H\psi^{\tau}(r,\theta)=E\psi^{\tau}(r,\theta)$. Noting that the Hamiltonian
commutes with the total angular momentum
$J_{z}=L_{z}+\frac{\hbar}{2}\sigma_{z}$, the radial and angular components of
the eigenspinors $\psi^{\tau}(r,\theta)$ can be separated. As a result, we
have

$\psi^{\tau}{(r,\theta)}=e^{im\theta}\begin{pmatrix}\Phi_{A}^{\tau}(r)\\\
ie^{-i\tau\theta}\Phi_{B}^{\tau}(r)\end{pmatrix}$ (8)

where the angular-momentum quantum number is $m=0,\pm 1,\pm 2,\cdots$.
$\Phi_{A}^{\tau}$ and $\Phi_{B}^{\tau}$ are wave functions that represent
sublattices $A$ and $B$, respectively. The parameter $\tau$ has two values,
$\pm 1$, which distinguishes the two valleys $K$ and $K^{\prime}$.

##  III Eigenspinors

Using the mathematical tools established above, we can now get the solutions
to the energy spectrum. Indeed, the equations for the region 1 ($r>r_{0}$) are
obtained

$\displaystyle\left(\partial_{\rho}+\frac{m\tau}{\rho}\right)\Phi^{\tau}_{A}(\rho)=-\left((\varepsilon-v)+\delta_{1}\right)\Phi^{\tau}_{B}(\rho)$
(9)
$\displaystyle\left(\partial_{\rho}-\frac{m-\tau}{\rho}\right)\Phi^{\tau}_{B}(\rho)=\left((\varepsilon-v)-\delta_{1}\right)\Phi^{\tau}_{A}(\rho)$
(10)

where the variable change $\rho=\frac{r}{r_{0}}$ and dimensionless units
$\varepsilon=\frac{E}{E_{0}}$, $v=\frac{V}{E_{0}}$,
$\delta_{1}=\frac{\Delta_{1}}{E_{0}}$ are considered, such that
$E_{0}=\frac{\hbar\nu_{F}}{r_{0}}$. By injecting (9) into (10), we get a
second order differential equation for $\Phi^{\tau}_{A}(\rho)$

$\left(\rho^{2}\partial^{2}_{\rho}+\rho\partial_{\rho}-m^{2}+\alpha^{2}_{1}\rho^{2}\right)\Phi^{\tau}_{A}(\rho)=0$
(11)

showing the Bessel function of the first kind, regular at the origin, as a
solution

$\Phi_{A}(\rho)=C^{\tau}_{1}J_{m}(\alpha_{1}\rho)$ (12)

where the parameter $\alpha_{1}$ is given by

$\alpha_{1}=\sqrt{|(\varepsilon-v)^{2}-\delta_{1}^{2}|}$ (13)

and $C^{\tau}_{1}$ is a normalization constant. The second component of spinor
can be derived from (9) as

$\Phi_{B}(\rho)=-iC^{\tau}_{1}e^{-i\tau\theta}J_{m-\tau}(\alpha_{1}\rho).$
(14)

Finally in region 1, the eigenspinors take the form

$\psi^{\tau}_{1}{(\rho,\theta)}=C^{\tau}_{1}e^{im\theta}\begin{pmatrix}J_{m}(\alpha_{1}\rho)\\\
-ie^{-i\tau\theta}J_{m-\tau}(\alpha_{1}\rho)\end{pmatrix}.$ (15)

Remember that there is a magnetic field $B$ in region 2 ($r<r_{0}$), which
causes the momentum operators $\pi^{\pm}$ to take the following forms

$\pi^{\pm}=-i\hbar\nu_{F}e^{\pm i\tau\theta}\left(\partial_{r}\pm\dfrac{\tau
i}{r}\partial_{\theta}\mp\frac{\tau eBr}{2\hbar}\right).$ (16)

Using $H\psi^{\tau}_{2}=E\psi^{\tau}_{2}$ to get

$\displaystyle\left(\partial_{\rho}+\frac{m\tau}{\rho}+\tau\beta\rho\right)\Phi^{\tau}_{A}(\rho)=-\left((\varepsilon-v)+\delta_{2}\right)\Phi^{\tau}_{B}(\rho)$
(17)
$\displaystyle\left(\partial_{\rho}-\frac{m-\tau}{\rho}-\tau\beta\rho\right)\Phi^{\tau}_{B}(\rho)=\left((\varepsilon-v)-\delta_{2}\right)\Phi^{\tau}_{A}(\rho)$
(18)

where $\beta=\frac{eBr_{0}^{2}}{2\hbar}$ and
$\delta_{2}={\frac{\delta_{2}r_{0}}{\hbar\nu_{F}}}$ have been set. After
substitution of (17) into (18), one obtains

$\left(\rho^{2}\partial^{2}_{\rho}+\rho\partial_{\rho}-m^{2}-2\beta(m-\tau)\rho^{2}+\alpha^{2}_{2}\rho^{2}-\beta^{2}\rho^{4}\right)\Phi^{\tau}_{A}(\rho)=0$
(19)

and we have involved

$\displaystyle\alpha_{2}=\sqrt{|(\varepsilon-v)^{2}-\delta_{2}^{2}|}.$ (20)

To solve (19), let us take the first component of spinor as

$\Phi_{A}(\rho)=\rho^{|m|}e^{-\frac{\beta\rho^{2}}{2}}\varphi(\beta\rho^{2})$
(21)

and define a new variable $\eta=\beta\rho^{2}$ to end up with the confluent
hepergeometric ordinary differential equation

$\left(\eta{\partial^{2}_{\eta}}+(b-\eta){\partial_{\eta}}-a\right)\varphi(\eta)=0$
(22)

with the quantities

$\displaystyle
a=-\frac{(\varepsilon-v)^{2}-\delta^{2}_{2}}{4\beta}+\frac{m-\tau+|m|+1}{2}$
(23) $\displaystyle b=1+|m|.$ (24)

As a result, we obtain the solution

$\Phi_{A}(\rho)=\rho^{|m|}e^{-\frac{\beta\rho^{2}}{2}}C^{\tau}_{2}M\left(a,b,\beta\rho^{2}\right)$
(25)

where $M(a,b,\beta\rho^{2})$ is the confluent hypergeometric function and
$C^{\tau}_{2}$ is a normalization constant. Thus, the second component can be
extracted from (17) as

$\Phi_{B}(\rho)=i\frac{\rho^{|m|}e^{-\frac{\beta\rho^{2}}{2}}e^{-i\tau\theta}}{(\varepsilon-v)+\delta_{2}}\left[\left(\frac{\tau
m}{\rho}+\beta\tau\rho\right)M\left(a,b,\beta\rho^{2}\right)-aM\left(a+1,b+1,\beta\rho^{2}\right)\right].$
(26)

When we combine everything, we get the eigenspinors in region 2

$\psi^{\tau}_{2}(\rho,\theta)=C^{\tau}_{2}\rho^{|m|}e^{-\frac{\beta\rho^{2}}{2}}e^{im\theta}\begin{pmatrix}M\left(a,b,\beta\rho^{2}\right)\\\
\frac{ie^{-i\tau\theta}}{(\varepsilon-v)+\delta_{2}}\left[\left(\frac{\tau
m}{\rho}+\beta\tau\rho\right)M\left(a,b,\beta\rho^{2}\right)-aM\left(a+1,b+1,\beta\rho^{2}\right)\right]\end{pmatrix}.$
(27)

##  IV Energy levels

We will proceed by using the boundary condition at the interface, $r=r_{0}$,
which is equivalent to $\rho=1$, because obtaining the energy level explicitly
for the current system is difficult. In fact, the operation

$\displaystyle\psi^{\tau}_{1}(1,\theta)=\psi^{\tau}_{2}(1,\theta)$ (28)

produces the relations

$\displaystyle
C^{\tau}_{1}J_{m}(\alpha_{1})=C^{\tau}_{2}e^{-}\frac{\beta}{2}M(a,b,\beta)$
(29) $\displaystyle
C^{\tau}_{1}J_{m-\tau}(\alpha_{1})=-C^{\tau}_{2}\frac{e^{-\frac{\beta}{2}}}{(\varepsilon-v)+\delta_{2}}\left[({\tau
m}+\tau\beta)M(a,b,\beta)-aM(a+1,b+1,\beta)\right].$ (30)

It is convenient to write this set of equations in the matrix from

$M^{\tau}\begin{pmatrix}C^{\tau}_{1}\\\ C^{\tau}_{2}\end{pmatrix}=0$ (31)

where we have defined

$M^{\tau}=\begin{pmatrix}J_{m}(\alpha_{1})&e^{-\frac{\beta}{2}}M(a,b,\beta)\\\
J_{m-\tau}(\alpha_{1})&-\frac{e^{-\frac{\beta}{2}}}{(\varepsilon-v)+\delta_{2}}[({\tau
m}+\tau\beta)M(a,b,\beta)-aM(a+1,b+1,\beta)]\end{pmatrix}.$ (32)

By requiring a null determinant, $\det M_{\tau}=0$, we can now determine the
equation governing the energy level. This process yields

$-J_{m}(\alpha_{1})\frac{e^{-\frac{\beta}{2}}}{(\varepsilon-v)+\delta_{2}}[({\tau
m}+\tau\beta)M(a,b,\beta)-aM(a+1,b+1,\beta)]-J_{m-\tau}(\alpha_{1})e^{-\frac{\beta}{2}}M(a,b,\beta)=0.$
(33)

Since this involves several physical parameters, a numerical analysis is
required to underline the basic features of the energy levels. In fact, we
will discuss the results shown in various plots of energy under appropriate
conditions.

##  V Results and discussion

The energy levels for a quantum dot in graphene as a function of radius
$r_{0}$ for a magnetic field $B=25$ T and a gap $\Delta_{2}=100$ meV inside
the quantum dot are shown in Fig. 2. We choose three angular momentum $m$
values (a,d,g): $m=-1$, (b,e,h): $m=0$, (c,f,i): $m=1$, as well as three gap
values outside the quantum dot (a, b,c): $\Delta_{1}=0$ meV, $\Delta_{1}=50$
meV, and $\Delta_{1}=150$ meV. Here, the blue curves represent valley $K$
($\tau=1$) and the dashed black curves represent valley $K^{\prime}$
($\tau=-1$). We see that for very small values of $r_{0}$, the energy levels
for both valleys, $K$ and $K^{\prime}$, are degenerate. They demonstrate the
symmetry $E(m,\tau)=-E(m,-\tau)$ as well as the asymmetry $E(m,\tau)\neq
E(m,-\tau)$. When the value of $r_{0}$ approaches $40$ nm, the energy levels
become linear and exhibit the symmetry $E(m,\tau)=E(m,-\tau)$. Furthermore,
for $\Delta_{1}=50$ meV, Figs. 2(d,e,f) show a number of new levels between
the valence and conductance bands, particularly for higher $r_{0}$ values. In
Figs. 2(g,h,i), we see that the number of levels increases and occupies the
area of $-\Delta_{1}+V<E<\Delta_{1}+V$. When compared to the results obtained
in [20, 22], these show new features due to the presence of the gap
$\Delta_{1}$. Figs. 2(g,h,i) show that the energy levels move vertically with
respect to the potential $V$, which is consistent with [20] observations.

(a)

(b)

(c)

Figure 2: (color online) Energy levels as a function of the quantum radius
$r_{0}$ for $B=25$ T, $\Delta_{2}=100$ meV and $m=-1,0,1$. (a,b,c):
$\Delta_{1}=0$ and V= $0$ meV, (d,e,f): $\Delta_{1}=50$ meV and V= $0$ meV,
and (g,h,i): $\Delta_{1}=150$ and V= $100$ meV. The magenta curves: $\tau=1$
and the dashed black curves: $\tau=-1$.

Fig. 3 shows the energy levels as a function of magnetic field $B$ for
$r_{0}=70$ nm, $\Delta_{2}=100$ meV, and $m=-1,0,1$, where (a,b,c):
$\Delta_{1}=0$ meV and $V=0$ meV, (d,e,f): $\Delta_{1}=50$ meV and $V=0$ meV,
(g,h,i): $\Delta_{1}=150$ meV and $V=100$ meV. The blue curves represent
valley $K$ ($\tau=1$), while the dashed black lines represent valley
$K^{\prime}$ curves ($\tau=-1$). We find that the energy levels vary according
to three regimes, depending on the values of the two gaps, $\Delta_{1}$ and
$\Delta_{2}$. In the first case $\Delta_{2}<E<-\Delta_{2}$, when $B$ is small,
the energy levels become degenerate. By increasing $B$ we notice that the
degeneracy is removed by maintaining the same symmetry as obtained in [20].
For the second case $\Delta_{1}<E<\Delta_{2}$ or $-\Delta_{2}<E<-\Delta_{1}$,
we observe that the energy levels are linear and non-symmetric. In the third
case, $-\Delta_{1}<E<\Delta_{1}$, new energy levels appear for very small
values of $B$, breaking the symmetry for the angular momentum $m\neq$ $0$, as
shown in Figs. 3(g,i). We can see that increasing $\Delta_{1}$ increases the
number of energy levels and fulfill the symmetry $E(m,\tau)=E(-m,-\tau)$.
These levels disappear as $B$ increases, and for $m=0$, regardless of $B$
value, there is an energy difference of $2\Delta_{1}$ between the valence band
and the conductance band.

(a)

(b)

(c)

Figure 3: (color online) Energy levels as a function of the magnetic field $B$
with $r_{0}=70$ nm and $m=-1,0,1$. (a,b,c): $\Delta_{1}=0$ meV and $V=0$ meV,
(d,e,f): $\Delta_{1}=50$ meV and V = $0$ meV, (g,h,i): $\Delta_{1}=90$ meV and
V = $100$ meV. Blue curves: $\tau$= $1$ and dashed black curves: $\tau$=$-1$.

The energy levels as a function of potential $V$ are shown in Fig. 4 for
$r_{0}=70$ nm, $\Delta_{2}=100$ meV, $B=10$ T, and $m=-1,0,{1}$, with (a, b,
c): $\Delta_{1}=50$ meV and (d, e, f): $\Delta_{1}=150$ meV. Valley $K$
($\tau=1$) is represented by the magenta curves, and valley $K^{\prime}$
($\tau=-1$) is represented by the dashed black curves. We observe that the
energy levels are clearly linear and have the symmetry $E(m,\tau)=E(m,-\tau)$.
There is no energy level for $m=0$ in the energy zone
$-\Delta_{1}<E<\Delta_{1}$. However, when $E$ is close to zero, there are
equidistant levels for $m\neq 0$. By increasing $\Delta_{1}$, we see that the
gap between the two bands increases. We see a coincidence of energy levels for
$E=\pm\Delta_{1}$. We conclude that the positions of the energy levels are
strongly dependent on the gap $\Delta_{1}$ outside the quantum dot.

Figure 4: (color online) Energy levels as a function of the potential $V$ with
$r_{0}=70$ nm, $\Delta_{2}=100$ meV, $B=10$ T and $m=-1,0,1$, with (a,b,c):
$\Delta_{1}=50$ meV and (d,e,f): $\Delta_{1}=150$ meV. Magenta curves:
$\tau=1$ and dashed black curves: $\tau=-1$.

Figure 5: (color online) Energy levels as a function of the gap $\Delta_{1}$
for $B=10$ T, $r_{0}=70$ nm, $\Delta_{2}=50$ meV and $m=-1,0,1$, with (a,b,c):
$V=0$ meV and (d,e,f): $V=100$ meV.

The variation of the energy levels as a function of the gap $\Delta_{1}$ is
shown in Fig. 5 for $B=10$ T, $r_{0}=70$ nm, $\Delta_{2}=50$ meV, and
$m=-1,0,1$, with two potential values (a,b,c): $V=0$ meV and (d,e,f): $V=100$
meV. The blue curves represent valley $K$ ($\tau$=1) and the dashed black
curves are for valley $K^{\prime}$ ($\tau$=-1). For -150 meV$>E>150$ meV (the
maximum value of $\Delta_{1}$), we observe that the energy levels vary in
horizontal parabolic form and have the symmetry $E(m,\tau)=E(m,\tau)$. There
are new vertical parabolic levels for -150 meV $<E<150$ meV and 50 meV
$<\Delta_{1}<-50$ meV that differ from those obtained in [21].

Figure 6: (color online) Energy levels as a function of the gap $\Delta_{2}$
for $B=10$ T, $r_{0}=70$ nm, $V=0$ meV and $m=-1,0,1$, with (a,b,c):
$\Delta_{1}=50$ meV and (d,e,f): $\Delta_{1}=100$ meV. Blue curves: $\tau=1$
and dashed black curves: $\tau$=-$1$.

In Fig. 6, we show the energy levels as a function of the gap $\Delta_{2}$ for
$B=10$ T, $r_{0}=70$ nm, $V=0$ meV and $m=-1,0,1$. with (a,b,c):
$\Delta_{1}=50$ meV and (d,e,f): $\Delta_{1}$=100 meV. Valley $K$ ($\tau=1$)
is represented by the blue curves, and valley $K^{\prime}$ ($\tau=-1$) is
represented by the dashed black curves. We see that the energy levels vary
linearly and follow the symmetry $E(m,\tau)=E(m,-\tau)$. When
$E>\Delta_{1}=50$ meV and $E<-\Delta_{1}=-50$ meV, there is a $2\Delta_{1}$
gap between the valence and conductance bands for $m=0$, which decreases as
$m\neq 0$. As shown in Fig. 6(d,e,f), increasing $\Delta_{1}$ reduces the gap
between two bands for $m=0$, but for $m\neq 0$, when
$-\Delta_{1}<E<\Delta_{1}$, there are energy levels satisfying the symmetry
$E(m,\tau)=E(-m,-\tau)$.

##  VI Conclusion

In a graphene quantum dot with an inhomogeneous gap in the electrostatic
potential and a perpendicular magnetic field, we have investigated the
dynamics of charge carriers. We solved the Dirac equation to obtain the
eigenspinors inside and outside the quantum dot. We were able to derive an
equation describing the energy levels in accordance with the physical
parameters that define our system by using the continuity condition at the
quantum dot’s interface. The obtained energy levels are rich, which pushed us
to underline their basic features.

Our numerical results are shown by studying the energy spectrum. Indeed, the
study of the spectrum’s dependence on the QD’s radius $r_{0}$ revealed that
the energy spectrum maintained the degeneracy of the asymmetry $E(m,\tau)\neq
E(m,-\tau)$ when the limit $r_{0}\longrightarrow 0$ is satisfied. The energy
levels become linear as $r_{0}$ increases, representing the symmetry
$E(m,\tau)=E(m,-\tau)$. Furthermore, we see new and different energy levels
occupying the region $-\Delta_{1}<E{<}\Delta_{1}$, and the number of these
levels increases as the gap $\Delta_{1}$ outside the QD is increased. In terms
of the energy spectrum’s dependence on the magnetic field $B$, we have seen
that the spectra for $\pm\Delta_{1}<E<\pm\Delta_{2}$ are linear even when $B$
is close to zero. We observed the appearing energy levels as $\Delta_{1}$ was
increased, with the exception of the angular momentum $m\neq$ 0 and they
satisfy the symmetry $E(m,\tau)=E(-m,{-\tau)}$ This means that the transport
properties of charge carriers at the center of the graphene magnetic QD are
affected by $\Delta_{1}$. Furthermore, the electrostatic potential $V$
influenced the energy level to shift vertically from $V$. Finally, we have
seen that the introduction of $\Delta_{1}$ in each spectrum representation has
broken the symmetry and resulted in new energy levels of the value energy
$-\Delta_{1}<E<\Delta_{1}$.

## References

  * [1] K. S. Novoslov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A. Firsov, Science 306, 666 (2004).
  * [2] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, M. I. Katsnelson, I. V. Grigorieva, S. V. Dubonos, and A. A. Firsov, Nature 438, 197 (2005).
  * [3] A. H. C. Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, Rev. Mod. Phys 81, 109 (2009).
  * [4] A. K. Geim and K. S. Novoselov, Nat. Mater. 6 183 (2007).
  * [5] J. V. Gomes and N. M. R. Peres, J. Phys: Condens. Matter 20, 325221 (2008).
  * [6] A. V. Rozhkov, G . Giavaras, Yury P. Bliokha, Valentin Freilikher, and Franco Nori, Physics Reports 503, 77 (2011).
  * [7] H. Y. Chen, V. Apalkov, and T. Chakraborty, Phys. Rev. Lett. 98, 186803 (2007).
  * [8] T. Espinosa-Ortega, I. A. Luk’yanchuk, and Y. G. Rubo, Phys. Rev. B 87, 205434 (2013).
  * [9] A . De Martino, L . Dell’Anna, and R. Egger, Phys. Rev. Lett. 98, 066802 (2007).
  * [10] G . Giavaras and Franco Nori, Appl. Phys. Lett. 97, 243106 (2010); ibid Phys. Rev. B 83, 165427 (2011); ibid 85, 165446 (2012).
  * [11] D. P. Zebrowski, E. Wach, and B. Szafran, Phys. Rev. B 88, 165405 (2013).
  * [12] M. R. Thomsen and T. G. Pedersen, Phys. Rev. B 95, 235427 (2017).
  * [13] P. Recher, J. Nilsson, G. Burkard and B. Trauzettel, Phys. Rev. B 79, 085407 (2009).
  * [14] Y. Liu, X. Liu, Y. Zhang, Q. Xia, and J. He, Nanotechnology 28, 235303 (2017).
  * [15] Claire Berger, Zhimin Song, Xuebin Li, Xiaosong Wu, Nate Brown, Cécile Naud, Didier Mayou, Tianbo Li, Joanna Hass, Alexei N. Marchenkov, Edward H. Conrad, Phillip N. First, and Walt A. de Heer, Science 312, 1191 (2006).
  * [16] M. Bacon, S. J. Bradley, and T. Nann, Part. Syst. Charact. 31, 415 (2013).
  * [17] B. Trauzettel, D. V. Bulaev, D. Loss, and G. Burkard, Nat. Phys 3, 192 (2007).
  * [18] H. Sun , L. Wu, W. Wei, and X. Qu, Mater. Today 16, 433 (2013).
  * [19] F. Belokda, A. Jellal, and E. Atmani, Phys. Lett. A 448, 128325 (2022).
  * [20] A. Belouad, B. Lemaalem, A. Jellal, and H. Bahlouli, Mater. Res. Express 7, 015090 (2020).
  * [21] A. Farsi, A. Belouad, and A. Jellal, Eur. Phys. J. B 94, 9 (2021).
  * [22] M. Mirzakhani, M. Zarenia, S. A. Ketabi, D. R. da Costa, and F. M. Peeters, Phys. Rev. B 93, 165410 (2016).


# An Enterprise Control Assessment Case Study of the Energy-Water Nexus for
the ISO New England System

Steffi O. Muhanji  Thayer School of Engineering, Dartmouth College, Hanover,
NH 03755<EMAIL_ADDRESS>Clayton Barrows  National
Renewable Energy Lab (NREL), Golden CO, (United States)
<EMAIL_ADDRESS>Jordan Macknick<EMAIL_ADDRESS>Amro M.
Farid<EMAIL_ADDRESS>
###### Abstract

The electric power generation mix of ISO New England (ISO-NE) is fundamentally
changing. Nuclear, coal, and oil generation facilities are retiring while
natural gas, solar, and wind generation are being adopted to replace them.
Variable renewable energy (VREs) such as solar and wind present multiple
operational challenges that require new and innovative changes to how the
electricity grid is managed and controlled. Within the context of a New
England case study, this paper looks at ways in which the water supply systems
(water and wastewater treatment), and water dependent electricity generating
resources (hydro, and thermal power plants) can be operated flexibly to help
balance energy in an evolving grid. The study’s methodology employs the novel
but now well published Electric Power Enterprise Control System (EPECS)
simulator to study the electric power systems operation, and the System-Level
Generic Model (SGEM) to study the associated water consumption and
withdrawals. This work studies six potential 2040 scenarios for the energy-
water nexus of the ISO-NE system. It presents a holistic analysis that
quantifies power system imbalances, normal operating reserves, energy market
production costs, and water withdrawal and consumption figures. For scenarios
with a high penetration of VREs, the study shows great potential of water
resources to enhance grid flexibility through the provision of load-following,
ramping, and regulation reserves by water resources. The work also provides
significant insights on how to jointly control the water and energy supply
systems to aid in their synergistic integration.

###### keywords:

Renewable Energy Integration, Energy-Water-Nexus , ISO New England ,
Curtailment , Reserves ,

††journal: Renewable and Sustainable Energy Reviews

##  1 Introduction

The bulk electric power system of New England is fundamentally changing to
include more solar and wind generation resources. This evolving resource mix
has triggered changes to how the power grid is managed and controlled. The
bulk of these changes have been in capacity and transmission expansion.
However, with the growing uncertainty and variability introduced by variable
renewable energy, there is an even greater need for increased amounts of
operational flexibility [1, 2]. Water plays a fundamental role in the ISO New
England system. Conventional and run-of-river hydro make up over 9% of the
overall generation in the 6 New England states[3]. An additional 1% of
generation comes from the two main pumped-water storage facilities, Bearswamp
and Northfield[3]. In the meantime, over 83% of the current ISO-NE generation
fleet comes from thermal generation facilities which withdraw and consume
large quantities of water for cooling purposes[3]. In spite of the changing
resource mix, recent studies predict that thermal generation facilities will
still account for a significant percentage of future generation facilities in
2040[4, 5]. Fig. 1 illustrates the extent of the coupling between the water
and electricity generation resources in New England. From Fig. 1, it is clear
that most generating facilities are located near a water source and rely on
adequate water supply to perform their function. These factors not only
indicate significant coupling between the water and electricity supply systems
but they also emphasize the need for more coordination between the two
systems. Specifically, the potential synergies between the two systems cannot
be ignored especially as the electricity grid undergoes its sustainable energy
transition.

Figure 1: A map of New England’s electric power generation units and rivers.
Figure 2: Architecture of the Electric Power Enterprise Control System (EPECS)
simulator customized for ISO New England operations [6]

Concern about water security is growing especially with climate change
affecting hydrology patterns and the decline of freshwater resources[7, 8, 9].
At the same time, significant attention has gone into the integration of
variable renewable energy into the electricity grid as a means of
decarbonizing the electricity supply system. As discussed in the prequel to
this paper[10], the challenges of renewable energy integration and energy-
water-nexus are very much related. In addition to presenting low $CO_{2}$
emissions, VREs have very low life-cycle water intensities[11]. On the other
hand, water is easily stored and therefore, has the potential to serve as a
flexible energy-water resource on both the supply-side as well as the demand-
side[12]. As a result, the methodologies of energy-water-nexus and renewable
energy integration studies must converge in order to realize potential
synergies.

With the growing penetration of wind and solar generation in the ISO-NE
system, grid operators would benefit immensely from an increased number of
flexible resources that can be used to balance the grid in the real-time.
Similarly, water system operators could offer ancillary services, improve
their profits and also achieve a more robust operation of their systems.
Despite these benefits, renewable energy integration and EWN studies have not
yet converged to realize benefits. While some energy-water nexus studies have
quantified the withdrawals by thermal power plants, these studies have largely
been conducted in isolation of actual operation of the electricity generation
industry[7, 9, 13, 14, 15, 16]. Thus, the full impact on either infrastructure
is not assessed. Other EWN works have focused solely on optimizing the
operations of water systems such as in the optimal operation of water pumps
and optimal pump scheduling[17, 18, 19, 20, 21] in order to provide demand
response and other ancillary services while maximizing returns for water
system operations[22, 23, 24, 25]. Finally, a small subset have presented
mostly single-layer approaches to co-optimize the water and electricity
networks. Examples of such works include the optimal network flow in [26], the
economic dispatch in [27], and the unit commitment problem in [28] for a
combined water, power, and co-production facilities. Despite the large body of
work and research on the energy water nexus, there is still a lack of a
generic, case and geography-independent methodologies that encompass all flows
within, and between the water and energy systems.

On the other hand, renewable energy integration studies have often been case
and geography specific and have mostly utilized unit-commitment-economic-
dispatch (UCED) models of power system control to study the operation of
electricity markets with large penetrations of VREs[29, 30, 31, 32, 6][33, 34,
35]. A significant percentage of these studies have taken statistical
approaches to determine the impact of wind and solar forecast errors on
dispatch decisions. A majority of renewable integration studies have
recognized the vital role of reserves in the balancing performance of systems
with high VRE penetration and have thus focused on the acquisition of normal
operating reserves such as load-following, regulation, and ramping
reserves[29, 30, 31, 32, 6].

However, a recent review of renewable integration studies shows major
methodological limitations[36]. Firstly, while some studies focus on reserve
acquisition, the required quantity of reserves is usually based on the
experiences of grid operators which no longer applies to systems with high
penetrations of VREs[37, 38]. Secondly, most studies only consider either the
net load variability or the forecast error in determining the amount of
reserves despite evidence that shows that both of these variables contribute
towards normal operating reserve requirements[39, 37]. Lastly, although
studies have shown that VREs possess dynamics that span multiple timescales of
power system operation[40, 41, 42], most renewable energy integration studies
have largely neglected the effect of timescales on the various types of
operating reserve quantities[36]. Farid et al. [36] proposed a holistic
approach based on enterprise control to study the full impact of VREs on power
system balancing performance and reserve requirements while considering the
multi-timescale dynamics of VREs. _Enterprise control_ is an integrated and
holistic approach that allows operators to study and improve the technical
performance of the grid while realizing cost savings[36]. An application of
enterprise control in the form of the Electric Power Enterprise Control System
(EPECS) simulator has been proposed in literature[43, 44, 45, 46, 47, 48] and
tested on various case studies including the ISO New England system[6]. In
[6], the EPECS simulator is used to study the performance of the ISO-NE system
on 12 scenarios with varying penetrations of VREs. This study highlights the
key role of curtailment and normal operating reserves on the balancing
performance of the ISO-NE system. This paper extends the work in [6] and [10]
to quantify the flexibility afforded the ISO-NE system through flexible
operation of water resources.

Figure 3: A diagram of the physical flows between the physical infrastructures
(water supply system, wastewater management system, and electricity supply
system) and the natural surface environment. [49]

###  1.1 Contribution

This paper applies the methodology presented in the prequel [10] and [6] to
study the techno-economic performance of the energy-water-nexus for the ISO-NE
system focusing on six predefined scenarios in 2040. The study methodology
takes the blueyellow rectangle of Fig. 3 as its system boundary. Given this
specific choice of system boundary, this study is able to quantify the mass
and energy flows in and out of the defined system boundary regardless of the
test case or geography. The paper also provides insight into some of the
operational challenges presented by high penetrations of VREs and also
quantifies the amounts of normal operating reserves for the ISO-NE system for
each scenario. Given that the methodology presented in the prequel [10] is
generic and modular, the EPECS simulator is modified slightly to reflect the
ISO-NE methodology (fully presented in [6] and as shown in Fig. 2). In this
study, the following quantities are studied: 1) load-following, ramping, and
regulation reserves, 2) the demand response potential of water units, 3) the
fuel flows of thermal units and their carbon emissions, 4) water withdrawals
and consumption by thermal power plants, and 5) the effect of flexible
operation of water resources on the production cost operation of the New
England electricity grid.

###  1.2 Outline

The rest of the paper is structured as follows: Section 2 presents the
methodology for the ISO New England Energy-Water Nexus study. Section 3 gives
a detailed description of the case study data. Section 4 presents the results
of the study within the context of the key performance characteristics of the
power grid. Finally, the paper is concluded in Section 5.

##  2 Methodology

As shown in Figure 4, the methodology of the ISO New England Energy-Water
Nexus study is best viewed in two parts: planning and operations. Section 2.1
describes how the National Renewable Energy Laboratory’s (NREL) Regional
Energy Deployment System (ReEDS) was used to evolve the 2030 ISO New England
electric power generation capacity mixes to six distinct 2040 capacity mix
scenarios. From there, the remainder of the section describes the Electric
Power Enterprise Control System (EPECS) simulator as customized for ISO New
England’s operation[6, 10]. Typically, it includes simulation functionality
for two energy market layers: the Security Constrained Unit Commitment (SCUC)
and the Security Constrained Economic Dispatch (SCED), power system regulation
and a physical model of the power grid itself (i.e. power flow analysis). For
this study, the simulator has been customized for ISO-NE operations to include
the Real-Time Unit Commitment (RTUC) as shown in Fig. 2. Furthermore, the SGEM
model[50, 10] is used to capture the essential physics of cooling processes
for thermal power plants and in turn compute the water withdrawals and
consumption for each power plant.

Figure 4: Diagram of the Simulators Used in the ISO New England Energy-Water
Nexus Study.

###  2.1 Regional Energy Deployment System (ReEDS) for Capacity Planning

ReEDS is a capacity planning tool that was developed by the National Renewable
Energy Laboratory (NREL) starting in 2003. ReEDS is a tool that identifies the
long-term evolution of the electric power grid for various regions in the
United States[51, 52, 53]. At its core ReEDS is an optimization tool that
identifies the cost-optimal mix of generation technologies subject to
reliability, generation resource, and regulatory constraints[51, 52, 53]. The
optimization has a two-year time step for a total of 42 years ending in
2050[51, 52, 53]. The final output of the simulation is generation capacity by
technology, storage capacity, electricity costs among others[51, 52, 53]. This
optimization tool was used to determine the evolution of the ISO-NE system
from the 2030 scenarios to the 2040 scenarios. The model input assumptions
were selected from configurations defined by the 2018 Standard Scenarios[54]
(see Table 1) to align with the 2030 capacity mixes described in Section 3.1.
Details on added capacities for each scenario can be found in Section 3.

Table 1: ReEDS 2018 standard scenarios[54] used to evolve the SOARES 2030 scenarios into the 2040 scenarios. |  SOARES 2030 Scenarios |  ReEDS Scenarios
---|---|---
1 |  RPSs + Gas |  High RE Cost
2 |  ISO Queue |  Accelerated Nuclear Retirements
3 |  Renewables Plus |  Low RE Cost
4 |  No Retirements beyond
Forward Capacity Auctions (FCA) #10 |  Low Wind Cost
5 |  ACPs + Gas |  Extended Cost Recovery
6 |  Renewable Portfolio Standards (RPSs) +
Geodiverse Renewables |  Low Natural Gas Prices

###  2.2 The Security Constrained Unit Commitment (SCUC)

The power system balancing operation commences with the day-ahead resource
scheduling in form of the SCUC. It is performed the day before to determine
the best set of generators that can meet the hourly demand at a minimum cost.
The time step for the SCUC is 1-hour and it determines the optimal set of
generators for the next 24-hours. A simplified version of this program is
presented in [10] and the full version customized for ISO-NE operations is
presented here[6]. Note that the SCUC formulation used for this study extends
the methodology in [6] to also include ramping constraints for wind, solar,
and hydro resources[10].

###  2.3 Real-Time Unit Commitment (RTUC)

The same day resource scheduling is conducted every hour through the RTUC. It
uses an optimization program that is quite similar to that of SCUC but only
commits and de-commits _fast-start_ units. Fast-start units are defined by
their ability to go online and produce at full capacity within 15-30 minutes.
The RTUC runs every hour with a 15-minute time step and a 4-hour look-ahead.
The complete mathematics for the RTUC can be found in [6] with slight
modifications to include ramping constraints for wind, solar, and hydro
resources as presented in [10].

###  2.4 The Security Constrained Economic Dispatch (SCED)

The real-time balancing operation is implemented through the SCED which is run
every 10-minutes. The role of the SCED is to move available generator outputs
to new set points in a cost-effective way. The SCED does not bring online any
units but rather ramps up or down the available online units. The SCED
methodology is presented in [6, 10] and similar to SCUC and RTUC, it has been
extended to allow for the ramping of wind, solar, and hydro resources[10]. A
more comprehensive description of the EPECS methodology and mathematical
formulations for each control layer can be found in [6, 10]. This methodology
has been analyzed and validated by ISO-NE.

###  2.5 Regulation

A pseudo-steady-state approximation of the regulation service model that ties
directly to a power flow model of the physical power grid is also used in this
study. Normally, imbalances at the output of the regulation service would be
represented in the form of frequency changes[55]. However, for steady-state
simulations with 1-minute time step, the concept of frequency is not
applicable. Instead, a designated _virtual_ swing bus consumes the mismatches
between generation and load to make the steady state power flow equations
solvable[6].

###  2.6 Variable Renewable Energy

Variable renewable energy resources in the EPECS simulator are studied as
time-dependent, spatially distributed exogenous quantities that contribute
directly to the net load. Where the term _net load_ here is defined as the
difference between the aggregated system load and the total generation
produced by VREs, tieline profiles and any transmission losses[6].

As previously defined in [10], the EPECS simulator differentiates energy
resources into several classes:

######  Definition 2.6.1.

_Variable Renewable Energy Resources (VREs)_ : Generation resources with a
stochastic and intermittent power output. Wind, solar, run-of-river hydro, and
tie-lines are assumed to be VREs.

######  Definition 2.6.2.

_Semi-Dispatchable Resources_ : Energy resources that can be dispatched
downwards (i.e curtailed) from their uncurtailed power injection value. When
curtailment is allowed for VREs, they become semi-dispatchable. In this study,
wind, solar and tie-lines are treated as semi-dispatchable resources. Note
that for the purposes of this study, run-of-river and conventional hydro
resources can be curtailed and, therefore, are treated as semi-dispatchable in
the “flexible case” mentioned below. However, in the conventional case, run-
of-river and conventional hydro resources are _not_ semi-dispatchable.

######  Definition 2.6.3.

_Must-Run Resources_ : Generation resources that must run at their maximum
output at all times. In this study, nuclear generation units are assumed to be
_must run_ resources.

######  Definition 2.6.4.

_Dispatchable Resources_ : Energy resources that can be dispatched up and down
from their current value of power injection. In this study, all other
resources are assumed to be dispatchable.

The EPECS simulator employs the operating reserve concepts described in [56,
57] with only a few changes. This study focuses on the normal operating
reserves that are able to respond to real-time changes in wind and solar
generation. Specifically, how much of these reserve quantities comes from
water resources such as conventional hydro, run-of-river hydro, and water and
waste-water treatment facilities. Normal operating reserves are classified as
load following, ramping, and regulation reserves based on the mechanisms upon
which they are acquired and activated. For the purposes of this study, the
curtailment of VREs was assumed to provide both load-following and ramping
reserves in an upward direction to their forecasted value and in a downward
direction to their minimum operating capacity limit.

Figure 5:  The ISO-NE zonal network model represented as “pipes” and
“bubbles”[6]

These three types of operating reserves work together to respond to real-time
forecast errors and variability in the _net load_ during normal system
operation. Note that the actual quantities of these reserves are physical
properties of the power system and exist regardless of whether they are
monetized or not. The EPECS simulator provides as output quantities: system
imbalances, operating reserves (load-following, ramping and regulation),
generator set points, curtailed generation and line flows for every minute.

###  2.7 System-level Generic Model (SGEM)

The S-GEM was developed to study the water use of fossil
fuel,nuclear,geothermal and solar thermal power plants using either steam or
combined cycle technologies[58, 59, 60, 61, 62, 63]. This model is also
geography and case-independent; making it ideal for application to the ISO-NE
system. Three main cooling processes are applied in this paper: once-through
cooling, wet tower cooling and dry-air cooling. Majority of the older
generation power plants used once-through cooling technology while the newer
power plants were either recirculating or dry-cooling. The formulae for
computing water withdrawals and consumption are presented in [10].

With this information, the energy-water flows through the system boundary of
Fig. 3 can be easily quantified (as detailed in [10]) to determine, water
withdrawals and consumption by thermal power plants, as well as other aspects
such as fuel consumption and $CO_{2}$ emissions. As illustrated in Figure 3,
it is important to capture all the physical flows between the three physical
infrastructures(water supply system, wastewater management system, and
electricity supply system) and the natural surface environment. In this study,
however, each water resource fits within an electric power system load area
(or “‘bubble” as they commonly called within the New England Power Pool).
Therefore, full hydraulic modeling does not provide additional insight in the
provision of flexibility services to the electric power grid. The approach
presented here is sufficient to capture all the interfaces between the water
supply system and the electricity supply system and impose aggregate energy
constraints as necessary.

##  3 Case Study Scenarios and Data

Figure 6: Summary of available generation capacity as a percentage of total
available capacity by fuel type for all six 2040 scenarios.

###  3.1 Study Scenarios

The case study scenarios presented in this work are best understood in the
context of the twelve scenarios that were studied in the 2017 System
Operational Analysis and Renewable Energy Integration Study (SOARES) that was
commissioned by ISO-NE. These 12 scenarios distinguished between the amount
and diversity of dispatchable generation resources, load profiles, and the
penetration of VREs[6]. Of these scenarios, six were meant to describe 2025
while the other six were meant to describe 2030. In the study presented here,
the six 2030 (SOARES) scenarios were evolved forward ten years using the ReEDS
capacity expansion software[51, 52, 53, 54]. The final capacity mixes of the
six scenarios are summarized in Figure 6 and are described further below.

In order to assess the value of coordinated vs uncoordinated energy-water
nexus operation, each of these six scenarios were simulated twice; once with
energy-water resources as variable resources and another as semi-dispatchable
resources. These scenario variants are respectively referred to as the
“conventional” operating mode (as a control case) and the “flexible” operating
mode (as the experimental case).

####  3.1.1 Scenario 1: RPSs + Gas

In this scenario, the oldest oil and coal generation units are retired by 2030
and the retired units are replaced by natural gas combined-cycle (NGCC) units
at the same locations. Furthermore, the ReEDS model adds 50 MW of biomass, 233
MW of solar, 75MW of hydro and 6351 MW of natural gas (NG) to this scenario.
It also retires 870 MW of nuclear, 667 MW of NG and 1127 MW of oil generation.

####  3.1.2 Scenario 2: ISO Queue

The retired oil and coal units from Scenario 1 are replaced by renewable
energy resources instead of NGCC. The locations of the renewable energy
resources are determined according to the ISO-NE Interconnection Queue. The
ReEDS model resulted in the addition of 2498 MW of solar, 9.77 MW of hydro,
and 5831.75 MW of NG (mostly in New Hampshire). In addition, 2471 MW of
nuclear, 668 MW of natural gas and 25 MW of coal generation units were
retired.

####  3.1.3 Scenario 3: Renewables Plus

In this scenario, more renewable energy resources are used to replace the
retiring units. Additionally, battery energy systems, energy efficiency and
plug-in hybrid electric vehicles (PHEV) are added to the system. Moreover, two
new tie lines are added to increase the amounts of hydroelectricity imports.
The ReEDS model results in the following modifications to this scenario: 1)
addition of 2760 MW of solar, 9 MW of hydro, 2765 MW of NG, and 2) the
retirement of 378 MW of coal, 870 MW nuclear, 667 MW of NG and 1127 MW of oil.

####  3.1.4 Scenario 4: No Retirements beyond Forward Capacity Auctions (FCA)
#10

In contrast to other scenarios, no generation units are retired beyond the
known FCA resources. The FCA resources are replaced by NGCC located at the
Hub. This scenario is the _business-as-usual_ scenario. The ReEDS model
results in the following modifications to this scenario: 1) addition of 989 MW
of solar, 4.2 MW of hydro, and 3987 MW of NG, and 2) the retirement of 383 MW
of coal, 870 MW nuclear, 667 MW of NG and 1127 MW of oil.

####  3.1.5 Scenario 5: ACPs + Gas

In this scenario, the oldest oil and coal generation units are retired by 2030
and these units are replaced by new NGCC units to meet the net Installed
Capacity Requirement (NICR). The ReEDS model results in the following
modifications to this scenario: 1) addition of 3089 MW of solar, 11.1 MW of
hydro, and 2496 MW of NG, and 2) the retirement of 253 MW of coal, 870 MW
nuclear, 667 MW of NG and 1127 MW of oil.

####  3.1.6 Scenario 6: Renewable Portfolio Standards (RPSs) + Geodiverse
Renewables

This scenario is similar to Scenario 5 but instead of replacing the retired
units with NGCC units, additional renewable energy generation is used to meet
the RPSs and the NICR. However, the solar PV and offshore wind units are
located closer to the main load centers while the onshore wind is located in a
remote area in Maine. The ReEDS model results in the following modifications
to this scenario: 1) addition of 3011 MW of solar, 6.2 MW of hydro, and 2430
MW of NG, and 2) the retirement of 870 MW nuclear, 667 MW of NG and 1127 MW of
oil.

The system data is consolidated into the zonal network model shown in Figure
5. The zonal network captures the power flows between pre-defined load zones
(i.e. “bubbles”) along abstracted “pipes”; thus eliminating the need for
Critical Energy/Electric Infrastructure Information (CEII) clearance. The
EPECS simulator implements a lossless DC Power Flow Analysis to determine
these flows as described in [6, 10]. The high-level interface flow limits
between the various bubbles are indicative of the line congestion often
experienced in the ISO New England territory. In addition to the changes in
capacity mixes implemented in REEDS, interface limits were raised to reflect
the likely situation that New England would work to resolve line congestion
found in the 2025 and 2030 scenarios in the SOARES scenarios[6]. Finally, in
addition to the electric data, data on power consumption by water and waste-
water treatment facilities as well as the cooling mechanisms of thermal
generators were used to determine their share of the peak load. The cooling
data for thermal power plants was further enhanced by data from the Energy
Information Agency’s (EIA) databases[64, 65, 66].

Figure 7: The load and net load profiles from Scenario 2040-4 (top) and 2040-3
(bottom). Figure 8: A comparison of load and net load distributions for all
six 2040 scenarios.

###  3.2 Net Load Profiles

The net load profile comprised of the system load profile minus the wind,
solar, tie-line, run-of-river and pond-hydro generation profiles. Figure 7
contrasts the net load profile of Scenario 2040-4 as a “business-as-usual”
case to that of Scenario 2040-3 as a high VRE case. The latter exhibits
significant negative net load especially during low load periods such as the
Spring and Fall seasons. Figure 8 summarizes the statistics of the net load
profiles for all six scenarios. The system _peak load_ for Scenarios
2040-1/2/4/5/6 was $28594MW$ while that of Scenario 2040-3 was $22103MW$ due
to a higher penetration of energy efficiency measures. All scenarios had the
same profile for electricity demand by water and wastewater treatment
facilities. Run-of-river and pond-hydro generation profiles were curtailable
at a price of $\$4.5/MWh$ similar to the 2017 ISO-NE SOARES. In this study,
flexible water resources have load-shedding rather load-shifting capability
and are assumed to contribute to operating reserves. The $709GWh$ of available
pumped storage capacity is treated as dispatchable for all six scenarios
throughout the study. Table 2 summarizes the capacity data for these flexible
energy-water resources. Again, in order to assess the ”flexibility value” of
these energy-water resources, each of the six scenarios is simulated in a
conventional-uncoordinated mode of operation as well as a flexible-coordinated
mode.

Table 2: A summary of available flexible water resources in the system as
percentage of the peak load.

##  4 Results

Given the aforementioned scenarios, the value of flexible energy-water
resources is assessed from reliability, economic, and environmental
perspectives. From a reliability perspective, Section 4.1 presents the
relative improvements in the system’s balancing performance as quantified by
the available quantities of operating reserves (i.e. load-following, ramping,
and regulation reserves), curtailment, and the magnitude of system imbalances.
From an environmental perspective, Section 4.2 quantifies the improvements in
the quantities of water withdrawn and consumed as well as CO2 emitted. Here,
_water withdrawn_ refers to the volumetric flow rate of water withdrawn from
the natural surface environment and _water consumption_ refers to the amount
of water not returned to its original point of withdrawal (due to evaporative
losses). Finally, Section 4.3 quantifies the associated production costs in
the day-ahead and real-time energy markets.

###  4.1 Balancing Performance of Coordinated Energy-Water Operation

As mentioned above, this section presents the system balancing performance
improvements as result of coordinated energy-water operation in terms of: the
available quantities of operating reserves (i.e. load-following, ramping, and
regulation reserves), curtailment, and the magnitude of system imbalances.

####  4.1.1 Load-Following Reserves

In day-to-day operation, the upward and downward load-following reserves are
used in time to allow the system to respond to variability and uncertainty in
the net load. In the traditional operation of the electricity grid, having
sufficient load-following reserves is a primary concern especially in systems
with high penetrations of renewables. Both upward and downward load-following
reserves are equally important in ensuring system reliability. As upward load
following reserves are exhausted (approach zero), the ability of the system to
respond to fluctuation in the net load is constrained.

Figure 9: Distributions of the available upward and downward load following reserves for all six 2040 scenarios in both the conventional and flexible operating modes. Table 3: Change in downward and upward load-following reserves statistics (flexible minus conventional) for 2040 scenarios. $\Delta$ LFR (MW) |  2040-1 |  2040-2 |  2040-3 |  2040-4 |  2040-5 |  2040-6
---|---|---|---|---|---|---
Up Mean |  208.1
(5.77%) |  171.7
(2.86%) |  65.6
(1.83%) |  207.1
(5.78%) |  194.2
(5.08%) |  57.7
(1.24%)
Up STD |  8.4
(1.00%) |  -55.6
(-1.94%) |  -17.3
(-1.22%) |  -42.1
(-5.32%) |  -67.6
(-8.36%) |  -36.09
(-1.74%)
Up Max |  178.3
(3.07%) |  228.3
(1.56%) |  335.3
(2.32%) |  242.5
(4.37%) |  107.9
(1.92%) |  686.8
(3.94%)
Up Min |  211.9
(14.03%) |  311.1
(22.77%) |  -96.3
(-12.45%) |  221.2
(15.12%) |  212.6
(15.50%) |  422.6
(40.46%)
Up 95 percentile1 |  241.1
(10.51%) |  282.7
(11.59%) |  6.0
(0.31%) |  288.9
(12.35%) |  294.6
(11.83%) |  244.5
(9.15%)
Down Mean |  743.8
(8.48%) |  801.6
(7.41%) |  925.5
(12.66%) |  647.2
(7.83%) |  744.0
(8.77%) |  984.1
(9.68%)
Down STD |  8.75
(0.36%) |  16.29
(0.66%) |  36.01
(1.52%) |  2.98
(0.12%) |  9.50
(0.39%) |  67.97
(2.55%)
Down Max |  1177.0
(6.11%) |  932.5
(4.37%) |  1678.0
(10.27%) |  961.1
(5.22%) |  1086.0
(5.79%) |  1424.0
(6.77%)
Down Min |  540.3
(16.53%) |  267.9
(5.75%) |  1019.0
(82.96%) |  720.5
(21.91%) |  554.9
(17.30%) |  583.2
(18.97%)
Down 95 percentile |  749.0
(13.96%) |  790.6
(10.79%) |  1026.0
(28.55%) |  717.7
(14.73%) |  750.7
(14.99%) |  876.3
(14.43%)

Therefore, an enhanced balancing performance with respect to load following
reserves would show a significant trough around the zero LFR-axis in the
distributions of load following reserves shown in Figure 9. The larger the
trough is, the more the system is not using its load following reserves to
balance the system. Figure 9 shows that the flexible use of energy-water
resources (in black) widens the trough of load-following reserves around the
zero line relative to conventional operation (in red). These graphical results
are confirmed numerically in Table 3. Flexible operation enhances the mean
values of the upward and downward load following reserves (treated as separate
distributions) by 1.24%– 12.66% across all six scenarios. Furthermore, the
minimum upward and downward load following reserves are improved by flexible
operation by 5.75% – 82.96% across all but one of the six scenarios. The
minimum statistic is particularly important because it defines a type of worst
case “safety margin” that the system will always have available to ensure its
security. Similarly the 95 percentile statistic gives a measure of how much
this minimum level increases when 5% of the distribution is treated as
abnormal outlier behavior. The simulations show improvements in the 95
percentile statistic of 0.13–28.55% across all six scenarios; thus
demonstrating its robustness to not just the minimum worst-case point but also
the distribution tail that represents challenging periods of operation. The
maximum and standard deviation statistics are provided for completeness.

####  4.1.2 Ramping Reserves

Ramping reserves describe the total amount of power that the system can
respond up or down within a minute. Traditionally, only dispatchable resources
are assumed to contribute towards ramping reserves. In this study, renewable
energy resources are semi-dispatchable by virtue of curtailment. Consequently,
they are assumed to not just be able to ramp down or up to their minimum or
maximum values but also do so within five minutes given their power-
electronics based control. Five minutes, in this case, coincides with the
minimum time-step used in the real-time market. Similar to load-following
reserves, ramping reserves are key to ensuring that the system can respond in
time to fluctuations in the net load. Having sufficient amounts of both upward
and downward ramping reserves is equally important to ensuring reliable
performance. As the amount of ramping reserves approaches zero, the ability of
the system to respond to net load variability is significantly diminished.

Figure 10: Distributions of the available upward and downward ramping reserves for all six 2040 scenarios in both the flexible and conventional operating modes. Table 4: Change in downward and upward ramping reserves statistics (flexible minus conventional) for all six 2040 scenarios. $\Delta$ RampR Stats (MW/min) |  2040-1 |  2040-2 |  2040-3 |  2040-4 |  2040-5 |  2040-6
---|---|---|---|---|---|---
Up Mean |  334.9
(11.83%) |  259.4
(5.28%) |  291.3
(8.26%) |  308.7
(13.78%) |  325.3
(14.31%) |  287.7
(6.16%)
Up STD |  14.8
(2.86%) |  27.9
(5.42%) |  3.5
(0.31%) |  16.3
(3.40%) |  11.6
(2.55%) |  15.8
(1.48%)
Up Max |  430.7
(10.40%) |  354.7
(5.65%) |  271.0
(4.83%) |  361.5
(10.43%) |  372.9
(10.58%) |  331.1
(4.79%)
Up Min |  -59.3
(-3.89%) |  69.7
(3.07%) |  410.6
(47.32%) |  -4.4
(-0.40%) |  -5.6
(-0.49%) |  305.1
(15.21%)
Up 95 percentile |  310.6
(14.77%) |  195.5
(4.68%) |  314.9
(14.11%) |  300.0
(18.78%) |  318.0
(19.19%) |  42.5
(1.28%)
Down Mean |  339.7
(14.81%) |  261.8
(5.86%) |  292.3
(8.70%) |  317.3
(18.35%) |  325.8
(17.88%) |  288.9
(6.50%)
Down STD |  16.4
(3.69%) |  21.4
(4.81%) |  1.5
(0.13%) |  16.1
(3.67%) |  12.7
(2.94%) |  12.4
(1.20%)
Down Min |  294.2
(22.51%) |  22.1
(1.06%) |  199.7
(31.65%) |  -15.1
(-1.92%) |  -6.7
(-0.76%) |  293.9
(18.44%)
Down Max |  417.3
(15.37%) |  354.3
(7.06%) |  275.9
(5.64%) |  385.1
(17.38%) |  345.1
(14.42%) |  320.7
(5.40%)
Down 95 percentile |  344.3
(19.12%) |  208.5
(5.31%) |  308.0
(13.94%) |  328.3
(26.15%) |  337.4
(24.92%) |  42.1
(1.32%)

Similar to load-following reserves, both upward and downward ramping reserves
are enhanced through the flexible operation of energy-water resources. Figure
10 illustrates a widened trough in the flexible operating mode relative to the
conventional mode. This observation is supported by the statistics in Table 4.
The mean value for the upward ramping reserves is improved across all
scenarios by up to 14.31%. Likewise, the mean downward ramping reserves are
improved by up to 18.35%. Another key measure of sufficient ramping reserves
is the minimum level. As illustrated in Table 4, flexible operation enhances
the minimum downward ramping reserves by 31.65% and the minimum upward ramping
reserves by a maximum of 47.32%. However, in cases with a lower penetration of
VREs such as scenarios 2040-1/4/5, the minimum levels are slightly worse in
the flexible case than in the conventional case. Despite these anomalies,
flexible operation improved 95% percentile levels of upward and downward
ramping reserves in all cases (by 1.28%–26.15%). These results show that the
curtailment of VREs increases the flexibility to the system if they are used
to provide ramping reserves. A complete summary of ramping reserves statistics
for all six scenarios is found in Table 4.

####  4.1.3 Curtailment

By definition, flexible energy-water resources increase the amount of
generation available for curtailment. Recall that by Definition (2.6.2), run-
of-river and conventional hydro-pond resources are semi-dispatchable resources
that can be curtailed in a flexible operating mode. As illustrated in Figure
11, scenarios with a lower penetration of VREs such as scenario 2040-1/4/5
curtail infrequently and the amount of megawatt curtailed is generally zero.
For scenarios 2040-2/3/6, curtailment is used at least 40% of the time.
Although, the two case appear to have similar curtailment levels, a closer
look at Table 5 shows that the flexible case curtails for a smaller percentage
of the year (2.67% – 10.9%) less than the conventional case). Furthermore, the
two operating modes show nearly identical levels of total curtailed energy. In
the absence of sufficient load-following and ramping reserves, curtailment
serves a key role in ensuring system balance. This role is particularly
crucial for VREs located in remote areas (e.g. Maine) where it serves as the
only control option given topological constraints and distance from load
areas.

Figure 11: Curtailment duration curves for all six 2040 scenarios in both the flexible (above) and conventional (below) operating modes.  Table 5: Change in the curtailment statistics (flexible minus conventional) for all six 2040 scenarios. | 2040-1 | 2040-2 | 2040-3 | 2040-4 | 2040-5 | 2040-6
---|---|---|---|---|---|---
Tot. Semi-Disp. Res. (GWh) | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00
Tot. Curtailed Semi-Disp.
Energy (GWh) | 17.71 | -1.95 | 60.86 | 23.44 | 20.57 | -6.18
% Semi-Disp. Energy Curtailed | 0.03 | -0.00 | 0.07 | 0.05 | 0.04 | -0.01
% Time Curtailed | -10.42 | -2.67 | -5.97 | -10.90 | -10.74 | -3.08
Max Curtailment Level (MW) | 1.82 | 2.68 | 330.16 | -63.03 | -1.81 | 397.67
Figure 12: Regulation duration curves for all six 2040 scenarios in both the
flexible (above) and conventional (below) operating modes.

####  4.1.4 Regulation Service

The regulation service is used to correct system imbalances in real-time. This
control lever is used to meet any left-over imbalances after curtailment,
load-following and ramping reserves have been used up during market operation.
In both cases, all scenarios appear to use their regulation effectively as
shown in Figure 12. This is indicative of a system that has sufficient
regulation to mitigate real-time imbalances and maintain balancing
performance. A closer inspection of Table 6 illustrates that flexible
operation marginally increases the reliance on regulation (as shown by the
excess mileage) and exhausts its regulation (albeit for a small fraction of
the year 0.001) for all but scenarios 2040-3 and 2040-4.

Table 6: Change in regulation reserves statistics (flexible minus conventional) for all six 2040 scenarios. | 2040-1 | 2040-2 | 2040-3 | 2040-4 | 2040-5 | 2040-6
---|---|---|---|---|---|---
% Time Reg. Res
Exhausted | 0.001 | 0.001 | 0.000 | 0.000 | 0.001 | 0.001
Reg. Res.
Mileage (GWh) | 1.800 | 0.354 | 0.788 | 1.014 | 1.190 | 0.468
% Reg. Res.
Mileage | 1.349 | 0.251 | 0.638 | 0.777 | 0.909 | 0.326
Figure 13: Range (above) and standard deviation (below) statistics for all six
2050 scenarios in both the flexible (red) and conventional (blue) operation
modes.

####  4.1.5 System Imbalances

Balancing performance indicates the residual imbalances after the regulation
service has been deployed. Given that the regulation service was barely
saturated, the amount of imbalances are expected to be minimal. As shown in
Figure 13, flexible energy-water resources had a small impact on the range of
final imbalances of the system. Both systems appear to perform similarly with
all cases maintaining a standard deviation of less than 16MW across all six
scenarios. Table 7 illustrates that the flexible operating mode performs
slightly better than the conventional with up to a 6.48% improvement in
standard deviation. The minimum imbalances are lower in all cases except for
Scenarios 2040-1 and 2040-2. Similarly, the maximum imbalances are lower for
the flexible operating mode except for Scenarios 2040-2 and 2040-6 which
represent scenarios with high VREs.

Table 7: Change in range and standard deviations of imbalances (flexible minus conventional) for all six 2040 scenarios. Change in Imbalance | 2040-1 | 2040-2 | 2040-3 | 2040-4 | 2040-5 | 2040-6
---|---|---|---|---|---|---
Max (MW) | -0.384 | 0.597 | -1.767 | -0.682 | -2.911 | 1.902
% Max | -0.164 | 0.241 | -0.998 | -0.297 | -1.269 | 0.779
Min (MW) | 0.118 | 1.831 | -0.598 | -0.363 | -4.405 | -0.462
% Min | -0.050 | -0.733 | 0.335 | 0.156 | 1.887 | 0.189
Std. (MW) | -0.552 | -0.611 | -0.684 | -0.589 | -0.584 | -0.634
% Std. | -3.874 | -4.052 | -6.484 | -4.188 | -4.147 | -4.155

###  4.2 Environmental Performance of Coordinated Energy-Water Operation

As mentioned before, the environmental performance of coordinated energy-water
operation is assessed through overall reductions in water withdrawals,
consumption and CO2 emissions.

####  4.2.1 Water Withdrawals

Figure 14 shows the water withdrawal distributions for the flexible and
conventional operating modes. Flexible operation results in significantly
lower withdrawals compared to conventional operation because the flexible
energy-water resources are able to offset the use of thermo-electric power
plants in favor of VREs. This phenomena is seen in how the flexible withdrawal
distributions are shifted left towards zero. The associate water withdrawal
statistics are summarized in Table 8 indicating improvements in mean
withdrawals of up to 25.58%. These improvements are most pronounced in
Scenarios 2040-2/3/6 with high penetrations of VREs. Indeed, the integration
of several percent (on capacity basis) of flexible energy-water resources as
shown in Table 2, serves to reduce water withdrawals by many multiples of that
percentage. Such a phenomena can potentially appear in any scenario where VRE
curtailment serves as a major lever of balancing control. Nevertheless, the
flexible operation of energy-water resources reduces water withdrawals across
all six scenarios.

Figure 14: Distributions of water withdrawals for all six 2040 scenarios in both the flexible and conventional operating modes. Table 8: Change in water withdrawals statistics (conventional minus flexible) for all six 2040 scenarios. $\Delta H_{2}0$ Withdrawals |  2040-1 |  2040-2 |  2040-3 |  2040-4 |  2040-5 |  2040-6
---|---|---|---|---|---|---
Mean ($m^{3}/min$) |  905.0
(0.70%) |  21370.0
(17.29%) |  24050.0
(20.59%) |  965.5
(0.74%) |  837.6
(0.65%) |  32460.0
(25.58%)
STD ($m^{3}/min$) |  106.7
(0.20%) |  714.1
(1.35%) |  -9537.0
(-19.92%) |  161.8
(0.31%) |  85.3
(0.16%) |  -12790.0
(-24.40%)
Max ($m^{3}/min$) |  1251.0
(0.45%) |  924.6
(0.34%) |  976.6
(0.39%) |  1290.0
(0.47%) |  1534.0
(0.56%) |  1289.0
(0.47%)
Min ($m^{3}/min$) |  40.1
(0.11%) |  27260.0
(88.22%) |  25630.0
(75.82%) |  431.1
(1.17%) |  575.7
(1.54%) |  26830.0
(75.99%)
Total ($m^{3}/min\times 10^{6}$) |  475.7 |  11230.0 |  12640.0 |  507.5 |  440.2 |  18090.0
Percent change ($\%$) |  0.70 |  17.29 |  20.59 |  0.74 |  0.65 |  25.58

####  4.2.2 Water Consumption

Electric power system water consumption occurs through the evaporative losses
from cooling towers in recirculating cooling systems. Figure 15 shows the
water consumption distribution for both the conventional and flexible
operating modes. While the effect is not large, the flexible mode of operation
shifts the distribution slightly towards the zero mark. Specifically, flexible
operation consumes 1.07–4.51% less water than the conventional operation
across all six scenarios, as shown in Table 9. This relatively small
percentage nevertheless accounts for $258\times 10^{3}m^{3}$ of water saved
every year. Scenarios 2040-3 and 2040-6 have the least savings. Due to high
penetrations of VREs, these scenarios require faster ramping generation which
mostly comes from fast-ramping natural gas units with recirculating cooling
systems. In short, the water saving effect of integrating VREs is a diminished
to a certain extent by the need for operating reserves from water-consuming
but flexible NGCC plants. If demand side resources (from water loads or
otherwise) played a large balancing role, then the water saving role of
integrating VREs would be more pronounced.

Figure 15: Distributions of water consumption for all six 2040 scenarios in both the flexible and conventional operating modes. Table 9: Change in evaporative loss statistics (conventional minus flexible) for all six 2040 scenarios. $\Delta$ Evap Losses |  2040-1 |  2040-2 |  2040-3 |  2040-4 |  2040-5 |  2040-6
---|---|---|---|---|---|---
Mean ($m^{3}/min$) |  2.67
(3.96%) |  1.63
(3.11%) |  0.30
(1.44%) |  3.37
(5.03%) |  1.51
(2.84%) |  0.31
(1.03%)
STD ($m^{3}/min$) |  1.10
(2.77%) |  1.05
(2.97%) |  0.74
(5.58%) |  1.23
(3.33%) |  0.61
(2.61%) |  0.68
(3.05%)
Max ($m^{3}/min$) |  5.71
(2.45%) |  3.42
(1.44%) |  6.40
(6.02%) |  -0.00
(-0.00%) |  1.80
(1.11%) |  0.07
(0.04%)
Min ($m^{3}/min$) |  -0.62
(-3.50%) |  -0.00
(-0.00%) |  -0.13
(-1.65%) |  0.47
(2.56%) |  -0.12
(-0.83%) |  -0.06
(-0.52%)
Total ($m^{3}\times 10^{3}$) |  1402 |  859 |  158 |  1769 |  794 |  165
Percent change ($\%$) |  4.12 |  3.21 |  1.46 |  5.30 |  2.92 |  1.03

####  4.2.3 $CO_{2}$ Emissions

Finally, as shown in Figure 16, the overall $CO_{2}$ emissions are
significantly reduced through flexible operation. It reduces the overall
$CO_{2}$ emissions by 2.10%–3.46%, as shown in Table 10. The mean, max, and
standard deviation of emissions are all improved. This CO2 emissions reduction
occurs because flexible energy-water resources 1.) eliminate the need for some
generation through reduced electricity consumption, 2.) enable greater VRE
generation through a reduction in curtailment and 3.) displace fossil-fueled
conventional generation.

Figure 16: Distributions of $CO_{2}$ emissions for all six 2040 scenarios in both the flexible and conventional operating modes. Table 10: Change in $CO_{2}$ emissions statistics (flexible minus conventional) for all six 2040 scenarios. $\Delta CO_{2}$ Emissions |  2040-1 |  2040-2 |  2040-3 |  2040-4 |  2040-5 |  2040-6
---|---|---|---|---|---|---
Mean ($kg$) |  82280
(3.46%) |  60330
(3.28%) |  21900
(3.17%) |  82390
(3.11%) |  71840
(2.90%) |  23120
(2.10%)
STD ($kg$) |  31460.0
(2.44%) |  32230.0
(2.66%) |  36350.0
(5.75%) |  30660.0
(2.69%) |  29540.0
(2.71%) |  28830
(2.96%)
Max ($kg$) |  51500
(0.71%) |  176000
(2.38%) |  222500
(5.54%) |  90040
(1.26%) |  121800
(1.72%) |  103100
(1.59%)
Min ($kg$) |  8189.00
(2.07%) |  -3313.00
(-1.08%) |  -2383.00
(-1.35%) |  -5755.00
(-1.14%) |  1179.00
(0.31%) |  92.23
(0.03%)
Total ($kg\times 10^{6}$) |  43240 |  31710 |  11510 |  43300 |  37760 |  12150
Percent change ($\%$) |  3.46 |  3.28 |  3.17 |  3.11 |  2.90 |  2.10

###  4.3 Economic Performance of Coordinated Energy-Water Operation

The economic performance of coordinated energy-water operation is assessed in
terms of the day-ahead and real-time production costs.

####  4.3.1 Day-Ahead Energy Market Production Costs

Figure 17 shows flexible operation reduced the total production cost in the
day-ahead energy market for all 2040 scenarios. Table 11 summarizes the
associated statistics. Flexible operation reduced total production costs by
29.3–68.09M$ or between 1.22–1.76%. As illustrated in Figure 17, Scenarios
2040-2/3/6 have much lower day-ahead production costs due to a high
penetration of VREs. In contrast, scenarios 2040-1/4/5 have significantly
higher costs as they are forced to commit expensive thermal power plants. In
short, the day-ahead energy market production costs are lower because the
flexible mode of operation represents an optimization program that is less
constrained than the program associated with the conventional mode of
operation.

Figure 17: Total production cost in the day-ahead energy market for all 2040 scenarios in both the flexible and conventional operating modes.  Table 11: Change in day-ahead energy market production cost statistics (flexible minus conventional) for all six 2040 scenarios. $\Delta$ Day-Ahead Costs |  2040-1 |  2040-2 |  2040-3 |  2040-4 |  2040-5 |  2040-6
---|---|---|---|---|---|---
Mean ($/hr) |  6115.1
(1.22%) |  5909.4
(1.49%) |  3345.2
(1.76%) |  7712.7
(1.41%) |  7773.1
(1.49%) |  4388.1
(1.64%)
STD ($/hr) |  4859.0
(2.09%) |  4355.7
(1.89%) |  5336.3
(3.89%) |  5327.3
(2.62%) |  6160.9
(3.05%) |  6095.2
(3.02%)
Max ($/hr) |  -16071.5
(-0.95%) |  38820.1
(2.65%) |  66093.4
(5.44%) |  -76701.8
(-4.56%) |  15683.0
(0.83%) |  476535.0
(23.20%)
Min ($/hr) |  19290.1
(18.95%) |  -2738.0
(-3.14%) |  15922.7
(19.18%) |  -706.4
(-0.45%) |  -419.0
(-0.36%) |  -10860.0
(-12.17%)
Total (million $) |  53.57 |  51.77 |  29.30 |  67.56 |  68.09 |  38.44
% Reduction |  1.22 |  1.49 |  1.76 |  1.41 |  1.49 |  1.64

####  4.3.2 Real-Time Energy Market Production Costs

Figure 18 illustrates the total real-time energy market production cost for
all six scenarios. Similar to the day-ahead energy market, Scenarios
2040-1/4/5 have significantly higher production costs as they are forced to
dispatch more expensive thermal power plants. Meanwhile, Scenarios 2040-2/3/6
have lower real-time energy market production costs due to a greater
utilization of renewable energy. As detailed in Table 12, flexible operation
reduces the average real-time market production costs by 2.46%–3.70% (or
19.58-70.83M$) across all six scenarios.

Figure 18: A comparison of the real-time production costs for flexible and conventional operation. Table 12: A summary of the real-time production cost statistics (flexible minus conventional). $\Delta$ Real-Time Cost |  2040-1 |  2040-2 |  2040-3 |  2040-4 |  2040-5 |  2040-6
---|---|---|---|---|---|---
Mean ($/min) |  1347.5
(3.70%) |  1013.5
(3.65%) |  372.5
(3.59%) |  1304.9
(3.12%) |  1173.1
(2.96%) |  412.5
(2.46%)
STD ($/min) |  493.5
(2.31%) |  533.2
(2.62%) |  553.8
(5.21%) |  497.8
(2.58%) |  545.8
(2.90%) |  536.9
(3.30%)
Max ($/min) |  895.8
(0.58%) |  3976.9
(2.69%) |  385.2
(0.36%) |  3163.4
(2.02%) |  -5845.8
(-3.41%) |  40662.3
(23.52%)
Min ($/min) |  88.4
(2.76%) |  75.5
(3.45%) |  -0.0
(-0.00%) |  65.3
(0.98%) |  -0.0
(-0.00%) |  157.3
(3.78%)
Total (million $) |  70.83 |  53.27 |  19.58 |  68.58 |  61.66 |  21.7
% Reduction |  3.70 |  3.65 |  3.59 |  3.12 |  2.96 |  2.46

##  5 Conclusion

This work has used a novel enterprise control assessment methodology to study
the energy-water nexus for the ISO New England System. Six scenarios were
studied representing plausible electric power capacity mixes in 2040. The
study specifically sought to understand the impact of flexible coordinated
operation of energy-water resources on the holistic behavior of these six
scenarios. In short, the flexible operation energy-water resources
demonstrated truly “sustainable synergies” with respect to balancing,
environmental, and economic performance. Table 13 summarizes the most
important results of the study in a balanced sustainability scorecard and
highlights the synergistic improvements caused by flexible coordinated
operation of the energy-water nexus. These results show that as VRE resources
become an ever-important part of the electric power system landscape, so too
must the electric power system evolve to engage energy-water resources as
control levers. In some cases, such resources – like hydro-power plants – are
mainstays of traditional operation. In other cases, particularly water utility
electric loads, these resources will have to evolve their operation to become
true electric power grid participants.

Table 13: Balanced Sustainability Scorecard: The range of _improvements_ caused by coordinated flexible operation of the energy-water nexus. Balancing Performance |
---|---
Average Load Following Reserves |  1.24–12.66%
Average Ramping Reserves |  5.28–18.35%
Percent Time Curtailed |  2.67–10.90%
Percent Time Exhausted Regulation Reserves |  0%
Std. Dev. of Imbalances |  3.874–6.484%
Environmental Performance |
Total Water Withdrawals |  0.65–25.58%
Total Water Consumption |  1.03–5.30%
Total CO2 Emissions |  2.10–3.46%
Economic Performance |
Total Day-Ahead Energy Market Production Cost |  29.30–68.09M$
Total Real-Time Energy Market Production Cost |  19.58–70.83M$

## References

  * [1] A. M. Farid, “The liines commitment to open-information,” The Labortory for Intelligent Integrated Networks of Engineering Systems, Tech. Rep., 2016. [Online]. Available: http://engineering.dartmouth.edu/liines/wpblog/2016/02/07/the-liines-commitment-to-open-information/
  * [2] EPRI, “Electric power system flexibility: Challenges and opportunities,” Electric Power Research Institute, Tech. Rep. 3002007374, February 2016.
  * [3] ISO-NE. (2019) Resource Mix-ISO New England. [Online]. Available: https://www.iso-ne.com/about/key-stats/resource-mix
  * [4] K. T. Sanders, M. F. Blackhurst, C. W. King, and M. E. Webber, “The Impact of Water Use Fees on Dispatching and Water Requirements for Water-Cooled Power Plants in Texas,” _Environ. Sci. Technol._ , vol. 48, no. 12, p. 140602120931006, Jun 2014.
  * [5] N. R. Armstrong, R. C. Shallcross, K. Ogden, S. Snyder, A. Achilli, and E. L. Armstrong, “Challenges and opportunities at the nexus of energy, water, and food: A perspective from the southwest united states,” _MRS Energy & Sustainability_, vol. 5, 2018.
  * [6] A. Muzhikyan, S. Muhanji, G. Moynihan, D. Thompson, Z. Berzolla, and A. M. Farid, “The 2017 ISO New England System Operational Analysis and Renewable Energy Integration Study,” _Energy Reports_ , vol. 5, pp. 747–792, July 2019.
  * [7] J. Rogers, K. Averyt, S. Clemmer, M. Davis, F. Flores-Lopez, D. Kenney, J. Macknick, N. Madden, J. Meldrum, S. Sattler, and E. Spanger-Siegfried, “Water-Smart Power: Strengthening the U.S. Electricity System in a Warming World,” Union for Concerned Scientists, Cambridge, MA, Tech. Rep., 2013.
  * [8] T. Kanyerere, S. Tramberend, A. D. Levine, P. Mokoena, P. Mensah, W. Chingombe, J. Goldin, S. Fatima, and M. Prakash, “Water futures and solutions: Options to enhance water security in sub-saharan africa,” in _Systems Analysis Approach for Complex Global Challenges_. Springer, 2018, pp. 93–111.
  * [9] K. Averyt, J. Macknick, J. Rogers, N. Madden, J. Fisher, J. Meldrum, and R. Newmark, “Water use for electricity in the united states: an analysis of reported and calculated water use information for 2008,” _Environmental Research Letters_ , vol. 8, no. 1, p. 015001, Jan 2013.
  * [10] S. O. Muhanji and A. M. Farid, “An Enterprise Control Methodology for the Techno-Economic Assessment of the Energy Water Nexus,” _(under revision)_ , vol. 1, no. 1, p. 31, [SPG-JR05] 2019. [Online]. Available: https://arxiv.org/abs/1908.10469
  * [11] M. Al-Nory and M. El-Beltagy, “An energy management approach for renewable energy integration with power generation and water desalination,” _Renewable Energy_ , vol. 72, pp. 377–385, 2014.
  * [12] W. N. Lubega and A. M. Farid, “A reference system architecture for the energy–water nexus,” _IEEE Systems Journal_ , vol. 10, no. 1, pp. 106–116, 2016.
  * [13] J. Meldrum, S. Nettles-Anderson, G. Heath, and J. Macknick, “Life cycle water use for electricity generation: a review and harmonization of literature estimates,” _Environmental Research Letters_ , vol. 8, no. 1, p. 015031, Mar 2013.
  * [14] J. Macknick, R. Newmark, G. Heath, and K. C. Hallett, “Operational water consumption and withdrawal factors for electricity generating technologies: a review of existing literature,” _Environmental Research Letters_ , vol. 7, no. 4, p. 045802, Dec 2012.
  * [15] K. Averyt, J. Fisher, A. Huber-Lee, A. Lewis, J. Macknick, N. Madden, J. Rogers, and S. Tellinghuisen, “Freshwater use by us power plants: Electricity’s thirst for a precious resource,” Union of Concerned Scientists, Cambridge, MA, USA, Tech. Rep., 2011.
  * [16] J. Macknick, S. Sattler, K. Averyt, S. Clemmer, and J. Rogers, “The water implications of generating electricity: water use across the united states based on different electricity pathways through 2050,” _Environmental Research Letters_ , vol. 7, no. 4, p. 045803, Dec 2012.
  * [17] S. A. Bagloee, M. Asadi, and M. Patriksson, “Minimization of water pumps’ electricity usage: a hybrid approach of regression models with optimization,” _Expert Systems with Applications_ , 2018.
  * [18] A. M. Bagirov, A. Barton, H. Mala-Jetmarova, A. Al Nuaimat, S. Ahmed, N. Sultanova, and J. Yearwood, “An algorithm for minimization of pumping costs in water distribution systems using a novel approach to pump scheduling,” _Mathematical and Computer Modelling_ , vol. 57, no. 3-4, pp. 873–886, 2013.
  * [19] B. Ulanicki, J. Kahler, and H. See, “Dynamic optimization approach for solving an optimal scheduling problem in water distribution systems,” _Journal of Water Resources Planning and Management_ , vol. 133, no. 1, pp. 23–32, 2007\.
  * [20] M. López-Ibáñez, T. D. Prasad, and B. Paechter, “Ant colony optimization for optimal control of pumps in water distribution networks,” _Journal of Water Resources Planning and Management_ , vol. 134, no. 4, pp. 337–346, 2008.
  * [21] Z. Ghelichi, J. Tajik, and M. S. Pishvaee, “A novel robust optimization approach for an integrated municipal water distribution system design under uncertainty: A case study of mashhad,” _Computers & Chemical Engineering_, vol. 110, pp. 13–34, 2018.
  * [22] C. Diaz, F. Ruiz, and D. Patino, “Modeling and control of water booster pressure systems as flexible loads for demand response,” _Applied Energy_ , vol. 204, pp. 106–116, 2017.
  * [23] S. Takahashi, H. Koibuchi, and S. Adachi, “Water supply operation and scheduling system with electric power demand response function,” _Procedia Engineering_ , vol. 186, pp. 327–332, 2017.
  * [24] R. Menke, E. Abraham, P. Parpas, and I. Stoianov, “Extending the envelope of demand response provision though variable speed pumps,” _Procedia Engineering_ , vol. 186, pp. 584–591, 2017.
  * [25] ——, “Demonstrating demand response from water distribution system through pump scheduling,” _Applied Energy_ , vol. 170, pp. 377–387, 2016.
  * [26] A. Santhosh, A. M. Farid, and K. Youcef-Toumi, “Optimal Network Flow for the Supply Side of the Energy-Water Nexus,” in _2013 IEEE International Workshop on Intelligent Energy Systems_ , Vienna, Austria, 2013, pp. 1–6. [Online]. Available: http://dx.doi.org.libproxy.mit.edu/10.1109/IWIES.2013.6698578
  * [27] A. Santhosh, A. M. Farid, A. Adegbege, and K. Youcef-Toumi, “Simultaneous Co-optimization for the Economic Dispatch of Power and Water Networks,” in _The 9th IET International Conference on Advances in Power System Control, Operation and Management_ , Hong Kong, China, 2012, pp. 1–6. [Online]. Available: http://dx.doi.org/10.1049/cp.2012.2148
  * [28] W. Hickman, A. Muzhikyan, and A. M. Farid, “The Synergistic Role of Renewable Energy Integration into the Unit Commitment of the Energy Water Nexus,” _Renewable Energy_ , vol. 108, no. 1, pp. 220–229, 2017. [Online]. Available: https://dx.doi.org/10.1016/j.renene.2017.02.063
  * [29] E. Ela, M. Milligan, B. Parsons, D. Lew, and D. Corbus, “The evolution of wind power integration studies: past, present, and future,” in _Power & Energy Society General Meeting, 2009. PES’09. IEEE_. IEEE, 2009, pp. 1–8.
  * [30] A. S. Brouwer, M. van den Broek, A. Seebregts, and A. Faaij, “Impacts of large-scale Intermittent Renewable Energy Sources on electricity systems , and how these can be modeled,” _Renewable and Sustainable Energy Reviews_ , vol. 33, pp. 443–466, 2014.
  * [31] H. Holttinen, M. O. Malley, J. Dillon, and D. Flynn, “Recommendations for wind integration studies – IEA task 25,” International Energy Agency, Helsinki, Tech. Rep., 2012.
  * [32] H. Holttinen, A. Orths, H. Abilgaard, F. van Hulle, J. Kiviluoma, B. Lange, M. OMalley, D. Flynn, A. Keane, J. Dillon, E. M. Carlini, J. O. Tande, A. Estanquiro, E. G. Lazaro, L. Soder, M. Milligan, C. Smith, and C. Clark, “Iea wind export group report on recommended practices wind integration studies,” International Energy Agency, Paris, France, Tech. Rep., 2013.
  * [33] GE-Energy, “New england wind integration study,” GE Energy and ISO New England, Schenectady, New York, Tech. Rep. May, 2010.
  * [34] E. Shlatz, L. Frantzis, T. McClive, G. Karlson, D. Acharya, S. Lu, P. Etingov, R. Diao, J. Ma, N. Samaan, V. Chadliev, M. Smart, R. Salgo, R. Sorensen, B. Allen, B. Idelchik, A. Ellis, J. Stein, C. Hanson, Y. V. Makarov, X. Guo, R. P. Hafen, C. Jin, and H. Kirkham, “Large-scale pv integration study,” Navigant Consulting, Las Vegas, NV, USA, Tech. Rep., 2011.
  * [35] R. Piwko, X. Bai, K. Clark, G. Jordan, N. Miller, and J. Zimberlin, “The Effects of Integrating Wind Power on Transmission System Planning, Reliability and Operations,” GE Energy, Schnectady, New York, Tech. Rep., 2005\.
  * [36] A. M. Farid, B. Jiang, A. Muzhikyan, and K. Youcef-Toumi, “The Need for Holistic Enterprise Control Assessment Methods for the Future Electricity Grid,” _Renewable & Sustainable Energy Reviews_, vol. 56, no. 1, pp. 669–685, 2015. [Online]. Available: http://dx.doi.org/10.1016/j.rser.2015.11.007
  * [37] H. Holttinen, M. Milligan, B. Kirby, T. Acker, V. Neimane, and T. Molinski, “Using standard deviation as a measure of increased operational reserve requirement for wind power by wind engineering using standard deviation as a measure of increased operational reserve requirement for wind power,” vol. 44, no. 0, 2008.
  * [38] A. Robitaille, I. Kamwa, A. H. Oussedik, M. de Montigny, N. Menemenlis, M. Huneault, A. Forcione, R. Mailhot, J. Bourret, and L. Bernier, “Preliminary impacts of wind power integration in the hydro-quebec system,” _Wind Engineering_ , vol. 36, no. 1, pp. 35–52, Feb 2012.
  * [39] B. C. Ummels, M. Gibescu, E. Pelgrum, W. L. Kling, and A. J. Brand, “Impacts of Wind Power on Thermal Generation Unit Commitment and Dispatch,” _IIEEE Transactions on Energy Conversion_ , vol. 22, no. 1, pp. 44–51, 2007\. [Online]. Available: http://dx.doi.org/10.1109/TEC.2006.889616http://ieeexplore.ieee.org/lpdocs/epic03/wrapper.htm?arnumber=4106021
  * [40] A. E. Curtright and J. Apt, “The character of power output from utility-scale photovoltaic systems,” _Progress in Photovoltaics: Research and Applications_ , vol. 16, no. 3, pp. 241–247, 2008.
  * [41] J. Apt and A. Curtright, “The spectrum of power from utility-scale wind farms and solar photovoltaic arrays,” Carnegie Mellon Electricity Industry Center Working Paper, Pittsburgh, PA, United states, Tech. Rep., 2007.
  * [42] F. Milano, _Power system modelling and scripting_ , 1st ed. New York: Springer, 2010. [Online]. Available: http://www.uclm.es/area/gsee/web/Federico/psat.htm
  * [43] A. Muzhikyan, T. Mezher, and A. M. Farid, “Power System Enterprise Control with Inertial Response Procurement,” _IEEE Transactions on Power Systems (in press)_ , vol. 33, no. 4, pp. 3735 – 3744, 2018. [Online]. Available: http://dx.doi.org/10.1109/TPWRS.2017.2782085
  * [44] A. Muzhikyan, A. M. Farid, and K. Youcef-Toumi, “Relative Merits of Load Following Reserves and Energy Storage Market Integration Towards Power System Imbalances,” _International Journal of Electrical Power and Energy Systems_ , vol. 74, no. 1, pp. 222–229, 2016. [Online]. Available: http://dx.doi.org/10.1016/j.ijepes.2015.07.013
  * [45] ——, “An Enterprise Control Assessment Method for Variable Energy Resource Induced Power System Imbalances Part 2: Results,” _IEEE Transactions on Industrial Electronics_ , vol. 62, no. 4, pp. 2459 – 2467, 2015. [Online]. Available: http://dx.doi.org/10.1109/TIE.2015.2395380
  * [46] ——, “An Enterprise Control Assessment Method for Variable Energy Resource Induced Power System Imbalances Part 1: Methodology,” _IEEE Transactions on Industrial Electronics_ , vol. 62, no. 4, pp. 2448–2458, 2015\. [Online]. Available: http://dx.doi.org/10.1109/TIE.2015.2395391
  * [47] ——, “A Power Grid Enterprise Control Method for Energy Storage System Integration,” in _IEEE Innovative Smart Grid Technologies Conference Europe_ , Istanbul, Turkey, 2014, pp. 1–6. [Online]. Available: http://dx.doi.org/10.1109/ISGTEurope.2014.7028898
  * [48] ——, “Enterprise Control Mitigation of Variable Energy Resource Induced Power System Imbalances,” in _9th Carnegie Mellon University Electricity Conference_ , Pittsburgh, PA, United states, 2014, pp. 1–41.
  * [49] H. Abdulla and A. M. Farid, “Extending the energy-water nexus reference architecture to the sustainable development of agriculture, industry & commerce,” in _First IEEE International Smart Cities Conference_ , Guadalajara, Mexico, 2015, pp. 1–7. [Online]. Available: http://dx.doi.org/10.1109/ISC2.2015.7366166
  * [50] M. J. Rutberg, A. Delgado, H. J. Herzog, and A. F. Ghoniem, “A system-level generic model of water use at power plants and its application to regional water use estimation,” in _ASME 2011 International Mechanical Engineering Congress and Exposition_. American Society of Mechanical Engineers, 2011, pp. 513–523.
  * [51] S. Cohen, J. Becker, D. Bielen, M. Brown, W. Cole, K. Eurek, W. Frazier, B. Frew, P. Gagnon, J. Ho, P. Jadun, T. Mai, M. Mowers, C. Murphy, A. Reimers, J. Richards, N. Ryan, E. Spyrou, D. Steinberg, Y. Sun, N. Vincent, and M. Zwerling, “Regional Energy Deployment System (ReEDS) Model Documentation: Version 2018,” National Renewable Energy Lab. (NREL), Golden, CO (United States), Tech. Rep. NREL/TP-6A20- 72023, April 2019.
  * [52] W. Short, P. Sullivan, T. Mai, M. Mowers, C. Uriarte, N. Blair, D. Heimiller, and A. Martinez, “Regional Energy Deployment System (ReEDS),” National Renewable Energy Lab.(NREL), Golden, CO (United States), Tech. Rep. NREL/TP-6A20-46534, December 2011.
  * [53] K. Eurek, W. Cole, D. Bielen, N. Blair, S. Cohen, B. Frew, J. Ho, V. Krishnan, T. Mai, B. Sigrin, _et al._ , “Regional Energy Deployment System (ReEDS) Model Documentation: Version 2016,” National Renewable Energy Lab.(NREL), Golden, CO (United States), Tech. Rep. NREL/TP-6A20-67067, November 2016.
  * [54] W. Cole, W. Frazier, P. Donohoo-Vallett, T. Mai, and P. Das, “2018 Standard Scenarios Report: A U.S. Electricity Sector Outlook,” National Renewable Energy Lab (NREL), Tech. Rep. NREL/TP-6A20-71913, November 2018.
  * [55] A. Gómez-Expósito, A. de la Villa Jaén, C. Gómez-Quiles, P. Rousseaux, and T. Van Cutsem, “A taxonomy of multi-area state estimation methods,” _Electric Power Systems Research_ , vol. 81, no. 4, pp. 1060–1069, Apr 2011.
  * [56] H. Holttinen, M. Milligan, E. Ela, N. Menemenlis, J. Dobschinski, B. Rawn, R. J. Bessa, D. Flynn, E. Gomez-Lazaro, and N. K. Detlefsen, “Methodologies to determine operating reserves due to increased wind power,” _IEEE Trans. Sustain. Energy_ , vol. 3, no. 4, pp. 713–723, Oct 2012.
  * [57] E. Ela, B. Kirby, E. Lannoye, M. Milligan, D. Flynn, B. Zavadil, and M. O’Malley, “Evolution of operating reserve determination in wind power integration studies,” in _Power and Energy Society General Meeting, 2010 IEEE_. IEEE, 2010, pp. 1–8.
  * [58] M. J. Rutberg, “Modeling water use at thermoelectric power plants,” Ph.D. dissertation, Massachussets Institute of Technology, 2012.
  * [59] Integrated Pollution Prevention and Control (IPPC) and E. Commission, “Reference document on the application of best available techniques to industrial cooling systems december 2001,” European Commission, Tech. Rep. December, 2001.
  * [60] ECOFYS, “– geo-localised inventory of water use in cooling processes, assessment of vulnerability and of water use management measures,” European Commission, End Report ENV.D1/SER/2013/0004, July 2014.
  * [61] Eurelectric, “Bat for cooling systems, working group ”environmental protection”,” Eurelectric, Tech. Rep., 1999.
  * [62] J. L. Tsou, J. Maulbetsch, J. Shi, and EPRI, “Power plant cooling system overview for researchers and technology developers,” EPRI, Tech. Rep., 2013.
  * [63] S.-Y. Pan, S. W. Snyder, A. I. Packman, Y. J. Lin, and P.-C. Chiang, “Cooling water use in thermoelectric power generation and its associated challenges for addressing water-energy nexus,” _Water-Energy Nexus_ , vol. 1, no. 1, pp. 26–41, 2018.
  * [64] E. I. A. (EIA). Form eia-923 detailed data with previous form data (eia-906/920). [Online]. Available: https://www.eia.gov/electricity/data/eia923/
  * [65] Energy Information Agency (EIA). Form eia-860 detailed data with previous form data (eia-860a/860b). [Online]. Available: https://www.eia.gov/electricity/data/eia860/
  * [66] (2019) Form eia-767 historical data files. [Online]. Available: https://www.eia.gov/electricity/data/eia767/