# Stochastic Finite Volume Method for Uncertainty Management in
Gas Pipeline Network Flows

Saif R. Kazi⋆, Sidhant Misra⋆, Svetlana Tokareva⋆, Kaarthik Sundar†, and
Anatoly Zlotnik⋆ This project was supported by the LDRD program project
“Stochastic Finite Volume Method for Robust Optimization of Nonlinear Flows”
at Los Alamos National Laboratory. Research conducted at Los Alamos National
Laboratory is done under the auspices of the National Nuclear Security
Administration of the U.S. Department of Energy under Contract No.
89233218CNA000001. Report No. LA-
<EMAIL_ADDRESS>Applied Mathematics &
Plasma Physics, Los Alamos National Laboratory, Los Alamos, NM
<EMAIL_ADDRESS>Information Systems & Modeling, Los Alamos National
Laboratory, Los Alamos, NM 87545

###### Abstract

Natural gas consumption by users of pipeline networks is subject to increasing
uncertainty that originates from the intermittent nature of electric power
loads serviced by gas-fired generators. To enable computationally efficient
optimization of gas network flows subject to uncertainty, we develop a finite
volume representation of stochastic solutions of hyperbolic partial
differential equation (PDE) systems on graph-connected domains with nodal
coupling and boundary conditions. The representation is used to express the
physical constraints in stochastic optimization problems for gas flow
allocation subject to uncertain parameters. The method is based on the
stochastic finite volume approach that was recently developed for uncertainty
quantification in transient flows represented by hyperbolic PDEs on graphs. In
this study, we develop optimization formulations for steady-state gas flow
over actuated transport networks subject to probabilistic constraints. In
addition to the distributions for the physical solutions, we examine the dual
variables that are produced by way of the optimization, and interpret them as
price distributions that quantify the financial volatility that arises through
demand uncertainty modeled in an optimization-driven gas market mechanism. We
demonstrate the computation and distributional analysis using a single-pipe
example and a small test network.

##  I Introduction

The increasing use of wind and solar electricity production concurrently with
growing reliance on gas-fired power generation for both base load and to
compensate for intermittent renewables results in more volatile and uncertain
demands on natural gas delivered by pipelines [1]. This issue has been
addressed in part by developing probabilistic scheduling methods based on
power grid physics and operations for calibrated planning of gas-fired
electricity generation and its fuel consumption [2]. For instance, stochastic
optimization methods were developed based on sampling demand and price
scenarios associated with wind uncertainty [3]. These methods take advantage
of the understanding that power grid physics are approximated well for
operations scheduling using direct current (DC) power flow modeling. In
contrast, it has been well known for decades that probabilistically
constrained optimization for gas pipeline network flow allocation is highly
challenging because of nonlinearity [4].

The notion of robust optimization involves the formulation of mathematical
programs to determine optimal decisions that are feasible under a range of
parameter values [5, 6]. Formulations based on scenario sampling have been
developed for managing demand uncertainty in operations of steady and
transient gas network flows [7, 8], including recent approaches that account
for composition uncertainty [9]. Recent formulations have proposed integrated
optimization formulations to jointly manage day-ahead uncertainty in a power
grid and interconnected gas pipeline, including compensation for interval [10]
and intertemporal [11] variability. Despite these advances, the more
fundamental steady-state optimal gas flow under demand uncertainty remains
challenging to analyze and scale to large networks.

The interconnections between physical and financial management of the power
grid and gas pipelines in the presence of fuel price volatility and renewable
generation uncertainty raise complex economic issues [12, 13]. Whereas
economics are connected to operations in the power grid through well-
established optimization-driven wholesale electricity markets [14], such
mechanisms for gas pipeline management have been explored much more recently
[15]. Deterministic optimization-driven market mechanisms for steady-state and
dynamic hourly scheduling and pricing of natural gas delivery have been
developed [16, 17], and such formulations were shown to scale well to
realistic systems [18]. However, the extension of market mechanisms for gas
pipelines to stochastic formulations have not been examined to date. It is
well known that optimization of gas pipeline flows under uncertainty is highly
complex, because of the substantial nonlinearity and computational complexity
involved [4, 7]. The challenge of modeling the effect of intertemporal
uncertainty in time-dependent simulation of gas pipelines has to an extent
been met by a new stochastic finite volume (SFV) method for uncertainty
quantification (UQ) in gas pipeline flows [19], which was shown to model the
propagation of uncertainty in associated initial boundary value problems.
Moreover, finite volume methods have been used in financial analysis of
markets that do not involve underlying physics, such as options pricing [20].

In this paper, we adapt the SFV representation for solutions of the nonlinear
gas flow equations with uncertain boundary conditions to optimize the steady-
state flow allocation in a pipeline system subject to uncertain demands. We
propose this type of uncertainty modeling because of its computationally
efficient representation of the probability densities of state variables
throughout a pipeline system, and in particular its amenability to higher-
order reconstructions in representing probabilistic inequality constraints.
The primary contribution of our study is a chance-constrained optimization
formulation that extends the problem of steady-state flow allocation for a gas
pipeline system [16], resulting in a single solution for compressor
configuration and optimized flows that extremizes an operational and/or
economic objective in expectation. By utilizing the proven method of chance-
constraints, our pipeline demand uncertainty management approach provides
probabilistic robustness guarantees that can be calibrated based on risk
tolerance. Furthermore, when the objective function is interpreted as an
expected economic value created for users of the pipeline system, we examine
the dual variables of the optimization solution and use them to reconstruct
functions that we interpret as the probability distributions of locational
values of gas.

The paper is structured as follows. We describe gas pipeline flow modeling and
the steady-state optimal gas flow optimization formulation in Section II. This
is extended to the chance-constrained stochastic steady-state optimal gas flow
formulation that accounts for uncertain loads using the SFV method in Section
III. In Section IV, we describe our computational implementation and two small
interpretable case studies, then briefly propose an economic interpretation of
the results, and finally conclude in Section V.

##  II Steady-State Optimal Gas Flow Problem

We formulate a general robust optimization problem for steady-state gas flow,
in which the objective minimizes average gas compressor power (for maximal
operating efficiency) and maximizes the expected economic value for users of
the system in the form of revenue collected by the pipeline operator. The
formulation admits situations in which gas consumptions at demand nodes are
specified, as well as problems where the gas delivery to one or more nodes is
an optimized variable. In the latter case, price-quantity bids for withdrawal
(by consumers) and injection (by suppliers) are provided, and we suppose that
a market authority solves the steady-state flow allocation problem subject to
uncertain flows that are unknown but within a specified range.

###  II-A Gas Flow Modeling

We consider a minimal turbulent flow approximation relating pressure to mass
flow for compressible gas in a pipe. The flow of gas through a pipe in the
turbulent regime can be described by a simplification of the Euler equations
[21],

$\displaystyle\partial_{t}\rho+\partial_{x}\varphi=0,\,\,\,\partial_{t}\varphi+\partial_{x}p=-\frac{\lambda}{2D}\frac{\varphi|\varphi|}{\rho},$
(1)

where $\varphi$ is per-area mass flux (kg/m2/s) of the gas, $p$ is pressure
(Pa), and $\rho$ is gas density (kg/m3), which are defined on a domain
$x\in[0,L]$ at each time $t\in[0,T]$. The parameters that determine flow
capacity are the pipe-specific non-dimensional friction factor $\lambda$, pipe
diameter $D$, and pipe length $L$, and the wave propagation speed $a$ that is
a property of the gas that depends on composition and temperature. Recall that
the per-area mass flux is $\varphi=\rho\\!\cdot\\!u$, and is related to the
total flow $\phi$ through a pipe by $\varphi=\phi/A$ where $A=\pi D^{2}/4$ is
the cross-sectional area of the pipe. We assume that the wave speed $a$ is
uniform system-wide, and that gas pressure $p$ and density $\rho$ are related
by the ideal gas law $p=a^{2}\rho$. We simplify this way to focus on
probabilistic modeling, and the results can be extended to more complex
settings such as non-ideal gases. In steady-state, mass flux $\varphi$ is
constant, so that equations (1) will have the form

$\displaystyle\partial_{x}\varphi=0,\qquad\partial_{x}p=-\frac{\lambda}{2D}\frac{\varphi|\varphi|}{\rho}.$
(2)

The second equation above defines the change of gas density along the pipe,
and integrating it along space yields

$(\rho(L))^{2}-(\rho(0))^{2}=-\beta\varphi|\varphi|,\qquad\beta=\lambda
L/(a^{2}D).$ (3)

###  II-B Pipeline Network Modeling

A pipeline network can be modeled as a set of edges $\mathcal{E}$ connected at
junctions in a set $\mathcal{V}$ where gas is withdrawn from or injected into
the network. Gas flow is actuated by compressors in a set
$\mathcal{C}\subset\mathcal{E}$, which are modeled as multiplicative pressure
boosters that preserve through-flow. Expressing equation (3) in terms of
squared pressures $\Pi_{i}$ and $\Pi_{j}$ at nodes $i,j\in\mathcal{V}$ and
flow $\phi_{ij}$ on pipe $(i,j)\in\mathcal{E}$ yields a flow equation for each
edge in the set of pipes $\mathcal{P}\subset\mathcal{E}$,

$\Pi_{j}-\Pi_{i}=\kappa_{ij}\phi_{ij}|\phi_{ij}|,\quad\forall(i,j)\in\mathcal{P},$
(4)

where $\kappa_{ij}=a^{2}\lambda_{ij}L_{ij}/(A_{ij}^{2}D_{ij})$. The pressure
boost achieved by compressor stations is modeled according to

$\Pi_{j}=\alpha_{ij}\Pi_{i},\quad\forall(i,j)\in\mathcal{C},$ (5)

where $\alpha_{ij}$ is a factor that relates the squares of the compressor
discharge and suction pressures. The network nodes $\mathcal{V}$ are
categorized either as slack (pressure) nodes in the set $\mathcal{V}_{s}$ or
nonslack (flow) nodes in the set $\mathcal{V}_{q}$. The slack nodes are
characterized by a given squared pressure $\Pi_{j}$, and nonslack nodes are
characterized by a given withdrawal flow $q_{j}$ (negative if an injection),
according to

$\Pi_{j}=p_{j}^{2},\,\,j\in\mathcal{V}_{s},\quad\text{and}\quad
q_{j}=d_{j}-s_{j},\,\,j\in\mathcal{V}_{q},$ (6)

where the values of $p_{j}$ and withdrawals $d_{j}\geq 0$ or injections
$s_{j}\geq 0$ must be specified. We suppose that each node $j\in\mathcal{V}$
is either consumer or supplier, so only one of $d_{j}$ or $s_{j}$ can be
positive. Conservation of flow is enforced at each node in the network as

$\sum_{i\in\partial_{+}j}\phi_{ij}-\sum_{k\in\partial_{-}j}\phi_{jk}=q_{j},\,\,\forall
j\in\mathcal{V},$ (7)

with the notation
$\partial_{+}j=\left\\{i\in\mathcal{V}\mid(i,j)\in\mathcal{E}\right\\}$ and
$\partial_{-}j=\left\\{k\in\mathcal{V}\mid(j,k)\in\mathcal{E}\right\\}$ used
to denote the sets of nodes connected to $j$ by incoming and outgoing edges,
respectively. To ensure mass flow balance for the system, pressure is fixed at
a slack node into which flow is free, which represents a large supply source
with an injection (negative withdrawal). The withdrawal $d_{j}$ or supply
$s_{j}$ at a flow node $j\in\mathcal{V}_{q}$ specified in (6) may be a
decision variable. If the withdrawal $d_{j}\geq 0$ is optimized, we write that
$j\in\mathcal{O}_{d}\subset\mathcal{V}_{q}$, and if supply $s_{j}\geq 0$ is
optimized, we write that $j\in\mathcal{O}_{s}\subset\mathcal{V}_{q}$.

###  II-C Deterministic Optimal Gas Flow

The deterministic steady-state optimal gas flow problem is comprised of an
objective function solved subject to the physical flow constraints and
boundary conditions defined in equations (4)-(7) as well as inequality
constraints, which we define below. These limitations are imposed due to
engineering and operating requirements. We suppose that the gas pressure at
each location $j\in\mathcal{V}$ in the network is bounded according to

$\Pi_{j}^{min}\leq\Pi_{j}\leq\Pi_{j}^{max},\,\,\,\forall j\in\mathcal{V},$ (8)

and that the compressor ratios are bounded according to

$1\leq\alpha_{ij}\leq\alpha_{ij}^{max},\,\,\,\forall(i,j)\in\mathcal{C}.$ (9)

The objective function balances two quantities, the first of which
approximates the cost of operating the pipeline by using energy to operate
compressors. Following our previous studies [22, 17], we approximate the
compressor power by the energy needed for adiabatic compression, which takes
the form

$W_{c}=\sum_{(i,j)\in\mathcal{C}}\eta_{ij}\phi_{ij}(\alpha_{ij}^{m}-1),$ (10)

where $0<m=(\gamma_{g}-1)/\gamma_{g}/2<1$ where $\gamma_{g}$ is the heat
capacity ratio of the gas [23], and $\eta_{ij}$ is a constant coefficient. The
second component of the objective function quantifies the economic value
produced by the pipeline for its users, which is defined as

$W_{e}=\sum_{j\in\mathcal{V}}(c^{d}_{j}d_{j}-c^{s}_{j}s_{j}),$ (11)

where $c_{j}^{d}$ is the bid price of a consumer and $c_{j}^{s}$ is the offer
price of a supplier. In practice, we suppose that the economic value $W_{e}$
is at least an order of magnitude greater than $W_{c}$ for an pipeline system
design. Finally, in a market mechanism for gas flow scheduling, we may
constrain any optimized demand flows by

$\displaystyle 0\leq d_{j}\leq d_{j}^{\max},\quad j\in\mathcal{O}_{d},$ (12a)
$\displaystyle 0\leq s_{j}\leq s_{j}^{\max},\quad j\in\mathcal{O}_{s}.$ (12b)

Synthesizing the objective function subject to equality and inequality
constraints yields the deterministic steady-state optimal gas flow problem:

$\displaystyle\min_{\alpha,d,s}$ $\displaystyle\quad W_{c}-W_{e}\,\,\text{as
in eq. \eqref{eq:obj_comp} and \eqref{eq:obj_econ}}$ (13a) s.t. pipe momentum
conservation (4) (13b) compressor actions (5) (13c) node flow balance (7)
(13d) pressure constraints (8) (13e) compressor ratio limits (9) (13f)
nomination limits (12) (13g)

The form of this optimization problem can be modified by collecting variables
and parameters into vectors. Let us define the vectors $\boldsymbol{c}_{d}$
and $\boldsymbol{c}_{s}$ of nodal demand and supply prices for optimized flows
$\boldsymbol{d}$ and $\boldsymbol{s}$ at nodes in the sets $\mathcal{O}_{d}$
and $\mathcal{O}_{s}$, respectively. Define also the vectors $\boldsymbol{q}$
of nodal gas withdrawals, $\boldsymbol{\kappa}$ of pipe resistances
$\kappa_{ij}$, $\boldsymbol{\phi}$ of flows through all pipes,
$\boldsymbol{\Pi}$ of all nodal pressures, and $\boldsymbol{\alpha}$ of all
compressor ratios. We then define the incidence matrix $\boldsymbol{A}$ of the
graph

$A_{ik}=\left\\{\begin{array}[]{ll}1&\text{edge $k=(j,i)$ enters node $i$},\\\
-1&\text{edge $k=(i,j)$ leaves node $i$},\\\ 0&\text{else}.\end{array}\right.$
(14)

The problem (13) can then be stated as

$\displaystyle\min_{\alpha,d,s}$ $\displaystyle\quad
W_{c}-(\boldsymbol{c}_{d}^{T}\boldsymbol{d}-\boldsymbol{c}_{s}^{T}\boldsymbol{s})$
(15a) s.t.
$\displaystyle\quad\Pi_{j}-\Pi_{i}=\kappa_{ij}\phi_{ij}|\phi_{ij}|,\qquad\forall(i,j)\in\mathcal{P}$
(15b)
$\displaystyle\quad\Pi_{j}=\alpha_{ij}\Pi_{i},\qquad\qquad\qquad\forall(i,j)\in\mathcal{C}$
(15c) $\displaystyle\quad\boldsymbol{A}\boldsymbol{\phi}=\boldsymbol{q}$ (15d)
$\displaystyle\quad\boldsymbol{\Pi}_{min}\leq\boldsymbol{\Pi}\leq\boldsymbol{\Pi}_{max}$
(15e)
$\displaystyle\quad\boldsymbol{1}\leq\boldsymbol{\alpha}\leq\boldsymbol{\alpha}_{max}$
(15f)
$\displaystyle\quad\boldsymbol{0}\leq\boldsymbol{d}\leq\boldsymbol{d}_{max},\,\,\boldsymbol{0}\leq\boldsymbol{s}\leq\boldsymbol{s}_{max}$
(15g)

Formulations similar to the optimization problem (15) have been examined in
previous studies [24, 25], and extended to transient flows [17] and gas
mixtures [22], and economic interpretations were also developed [16, 26].
Steady-state formulations such as problem (13) can be viewed as providing a
day-ahead flow allocation for an operating day, assuming steady ratable
offtakes by customers. Here we extend the steady-state optimal gas flow
problem to account for uncertain demands using the SFV approach.

##  III Chance-Constrained Optimal Gas Flow

In the steady-state flow setting, we aim to develop a chance-constrained
formulation to provide probabilistic guarantees given a known distribution in
one or more gas flow withdrawals from a pipeline network. This type of
formulation can be more flexibly calibrated than robust optimization
approaches that compensate for arbitrary interval uncertainty [24, 27]. The
key concern for pipeline systems is to assure delivery of gas to all customers
given uncertainty in some loads while maintaining adequate pipeline pressure.
The lower bound in constraint (8) (i.e., (15e)) is somewhat flexible in
practice, as long as it is not violated significantly and pressures can be
returned to adequate levels for the subsequent operating day. We enforce this
limit using chance constraints, allowing for a small probability $\epsilon$ of
violation.

Suppose that a subset of gas pipeline nodes $\mathcal{S}\subset\mathcal{V}$
has stochastic gas consumptions of the form

$d_{j}(\omega)=d_{j}+r_{j}(\omega),\quad j\in\mathcal{S},$ (16)

where we suppose that $(\Omega_{j},\mathcal{B}_{j},\mu_{j})$ is a probability
space where $r_{j}:\Omega_{j}\to R_{j}$ is a random variable taking values on
a compact interval
$R_{j}=[\underline{r}_{j},\overline{r}_{j}]\subset\mathbb{R}$ and where
$\mathcal{B}_{j}$ is the Borel $\sigma$-algebra. The fixed parameters $d_{j}$
denote the nominal baseline load. Letting
$\Omega=\prod_{j=1}^{|\mathcal{S}|}\Omega_{j}$, we consider samples
$\omega\in\Omega$ as instances of stochastic loads across the pipeline
network. We suppose that values of random parameters $r_{j}(\omega)$
corresponding to various values of $\boldsymbol{q}$ in problem (15) would
result in different optimal solutions of that problem. Thus in the stochastic
extension below we consider the physical solutions $\boldsymbol{\Pi}$ and
$\boldsymbol{\phi}$ to also depend on the sample $\omega$, and we seek to find
a single decision for the value of $\boldsymbol{\alpha}$ and any optimized
flows $d_{j}$ or $s_{j}$ for $j\in\mathcal{O}\subset\mathcal{V}$ such that the
inequality constraints corresponding to (15e) are satisfied in a probabilistic
sense. The flow equations (4) and compressor conditions (5) must be evaluated
for all samples $\omega$ as

$\displaystyle\\!\\!\\!\\!\\!\\!\Pi_{j}(\omega)\\!-\\!\Pi_{i}(\omega)\\!=\\!\kappa_{ij}\phi_{ij}(\omega)|\phi_{ij}(\omega)|,\,\forall(i,j)\\!\in\\!\mathcal{P},\forall\omega\\!\in\\!\Omega,\\!\\!\\!\\!$
(17)
$\displaystyle\Pi_{j}(\omega)=\alpha_{ij}\Pi_{i}(\omega),\quad\forall(i,j)\in\mathcal{C},\,\,\forall\omega\in\Omega.$
(18)

Conservation of flow (7) is enforced for all samples as

$\sum_{i\in\partial_{+}j}\phi_{ij}(\omega)\\!-\\!\\!\\!\sum_{k\in\partial_{-}j}\phi_{jk}(\omega)\\!=\\!q_{j}(\omega),\,\,\forall
j\\!\in\\!\mathcal{V},\,\,\forall\omega\\!\in\\!\Omega,$ (19)

where $q_{j}(\omega)=d_{j}+r_{j}(\omega)-s_{j}$. The chance constraint for
minimum pipeline pressure is expressed using a quadratic penalty function

$\Gamma(z)=\begin{cases}\gamma z^{2},\quad&\mbox{if }x\geq 0,\\\
0,\quad&\mbox{otherwise},\end{cases}$ (20a) and a penalized minimum pressure
violation variable
$v^{min}_{j}(\omega)=\Gamma\left(\Pi_{j}^{\min}-\Pi_{j}(\omega)\right).$ (20b)

The quadratic penalty is chosen in order to reflect the practical
consideration that larger violations of the lower pressure limits are more
problematic. Additionally, because it is twice differentiable, it facilitates
a well-behaved nonlinear program. The chance constraint is then

$\mathbb{E_{\omega}}\\!\left[v^{min}_{j}(t,\omega)\right]\leq\epsilon_{j},\qquad\forall
t\in T,$ (21)

where the acceptable violation probability $\epsilon_{j}$ may depend on the
network node $j\in\mathcal{V}$. We approximate the constraint (21) in a
deterministic manner (i.e., without requiring Monte Carlo simulation) by
applying the SFV method as follows. We discretize the stochastic space
$\Omega$, which has one-to-one correspondence with the interval $R_{j}$ on
which the random consumption $r_{j}(\omega)$ appears, into $M$ cells delimited
by $M+1$ uniformly spaced boundary points, each of which corresponds to a
value of $r_{j}{\omega_{m}}$. We construct the penalty variable using a third
order spline expansion on $\Omega$ of form

$v^{min}_{j}(\omega)=\Gamma\left(\Pi_{j}^{\min}-\Pi_{j}(\omega)\right)=\sum_{m\in
M}a_{jm}b_{jm}(\omega),$ (22a) where $b_{jm}(\omega)$ is the $m$-th spline
function on the uniform stochastic space grid that is completely known. The
constraint (21) is expressed as an expectation over $\omega$ as
$\mathbb{E}_{\omega}[v^{min}_{j}(\omega)]=\sum_{m\in
M}a_{jm}\int_{\omega}b_{jm}(\omega)d\mu_{j}\leq\epsilon_{j},$ (22b)

where the coefficients $a_{jm}$ serve as decision variables that can be
optimized. Aside from minimum pressure limits, all other constraints including
the maximum operating pressure,

$\Pi_{j}(\omega)\leq\Pi_{j}^{max},\,\,\,\forall j\in\mathcal{V},$ (23)

are strictly enforced for all values of $\omega$. The objective function
components must also be evaluated in expectation,

$\displaystyle\mathbb{E}_{\omega}[W_{c}]$
$\displaystyle=\sum_{(i,j)\in\mathcal{C}}\eta_{ij}(\alpha_{ij}^{m}-1)\cdot\mathbb{E}_{\omega}[\phi_{ij}],$
(24) $\displaystyle\mathbb{E}_{\omega}[W_{e}]$
$\displaystyle=\sum_{j\in\mathcal{V}}(c^{d}_{j}(d_{j}+\mathbb{E}_{\omega}[r_{j}])-c^{s}_{j}s_{j}).$
(25)

Assembling the objective function and probabilistic constraints in the
stochastic setting results in the chance-constrained optimal gas flow problem
given by

$\displaystyle\min_{\alpha,d,s}$
$\displaystyle\quad\mathbb{E}_{\omega}[W_{c}]-\mathbb{E}_{\omega}[W_{c}]\,\,\text{as
in eq. \eqref{eq:obj_comp_exp} and \eqref{eq:obj_econ_exp}}$ (26a) s.t. pipe
momentum conservation (17) (26b) node flow balance (19) (26c) compressor
actions (18) (26d) pressure constraints (22a)-(22b) and (23) (26e) compressor
ratio limits (9) (26f) nomination limits (12) (26g)

Invoking the vector notations used to express formulation (13) as (15) above,
and using $\boldsymbol{r}$ to state the vector of random withdrawals provided
at prices $\boldsymbol{c}_{r}$ for $j\in\mathcal{S}$, we re-write problem (26)
as

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\min_{\alpha,d,s}$
$\displaystyle\,\,\mathbb{E}_{\omega}[W_{c}]-(\boldsymbol{c}_{d}^{T}\boldsymbol{d}-\boldsymbol{c}_{s}^{T}\boldsymbol{s})-\boldsymbol{c}_{r}^{T}\mathbb{E}_{\omega}[\boldsymbol{r}]$
(27a) s.t.
$\displaystyle\,\,\Pi_{j}(\omega)\\!-\\!\Pi_{i}(\omega)\\!=\\!\kappa_{ij}\phi_{ij}(\omega)|\phi_{ij}(\omega)|,$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\,\,\forall(i,j)\in\mathcal{P},\,\,\forall\omega\in\Omega,$
(27b)
$\displaystyle\,\,\Pi_{j}(\omega)=\alpha_{ij}\Pi_{i}(\omega),\quad\quad\forall(i,j)\in\mathcal{C},\,\,\forall\omega\in\Omega,$
(27c)
$\displaystyle\,\,\boldsymbol{A}\boldsymbol{\phi}(\omega)=\boldsymbol{q}(\omega),\quad\forall\omega\in\Omega$
(27d)
$\displaystyle\,\,\boldsymbol{\Pi}(\omega)\leq\boldsymbol{\Pi}_{max},\quad\forall\omega\in\Omega$
(27e)
$\displaystyle\,\,\boldsymbol{1}\leq\boldsymbol{\alpha}\leq\boldsymbol{\alpha}_{max}$
(27f)
$\displaystyle\,\,\Gamma\left(\Pi_{j}^{\min}-\Pi_{j}(\omega)\right)=\sum_{m\in
M}a_{jm}b_{jm}(\omega),$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\,\,\forall
j\in\mathcal{V},\,\forall\omega\in\Omega,\,\,$ (27g) $\displaystyle\sum_{m\in
M}\\!\\!a_{jm}\int_{\omega}b_{jm}(\omega)d\mu_{j}\leq\epsilon_{j},\,\,\forall
j\\!\in\\!\mathcal{V},\,\forall\omega\\!\in\\!\Omega,\,\,\,\,\,\,$ (27h)
$\displaystyle\,\,\boldsymbol{0}\leq\boldsymbol{d}\leq\boldsymbol{d}_{max},\,\,\boldsymbol{0}\leq\boldsymbol{s}\leq\boldsymbol{s}_{max}$
(27i)

In the rest of the paper, we examine the computational solution of problem
(27) for to small test networks by estimating the distributions of physical
solutions. Because the optimization objective (27a) is formulated as an
expected economic value, we also examine the possibility of reconstructing the
distributions of dual variables and interpreting them as stochastic locational
values of natural gas.

##  IV Computational Studies

We explore solutions to problem (27) obtained using our SFV-based uncertainty
management method using two example systems. The first consists of a single
pipe with a compressor at the sending end, and the second is an 8-node test
network that consists of 5 pipes and 3 compressor stations, and which was
described in our previous studies [28, 19, 22]. We examine the distributions
of physical flows and pressures that arise given uniform or truncated normal
probability distributions in withdrawal gas flows. The problem (27) is
implemented as a nonlinear constrained optimization model in JuMP [29], an
algebraic modeling language in Julia, and solving using an open source
interior point method based optimization solver IPOPT [30]. For computational
well-conditioning, we re-scale the pressure, flow, distance, and velocity
variables according to the non-dimensionalization in a previous study [31]. We
suppose for simplicity in this study that the value of the exponent $m$ in
equation (10) is $m=1$, so that the specific heat capacity ratio of the gas is
0.5.

Figure 1: Single pipe test system and component parameters.

###  IV-A Single Pipe with Compressor

The first test network consists of a single pipe with a compressor connected
to the inlet node, as shown in Figure 1. The injection node (N1) acts as a
slack node with fixed pressure set at $p^{slack}_{1}=4.3367$ MPa. The pressure
bounds at the withdrawal node (N3) are specified as $P_{min}=4.0$ MPa and
$P_{max}=6.0$ MPa. The compressor ratio is bounded by $\alpha\in[1,1.4]$. The
uncertain withdrawal flow is parameterized by an uncertain parameter
($\omega$) as $q_{3}(\omega)=d_{3}+r_{3}(\omega)$, where $d_{3}$ represents
the nominal value that is fixed at 250 kg/s. We consider a scenario where the
uncertain flow $r_{3}(\omega)$ follows a uniform distribution
$r_{3}\sim\mathcal{U}[-50,50]$, and a scenario where $r_{3}(\omega)$ follows a
truncated normal distribution with zero mean and standard deviation
$\sigma=50/3$, i.e., $r_{3}\sim\mathcal{N}(0,50/3)$ with support on
$[-50,50]$. The distributions of the withdrawal $q_{3}(\omega)$ are shown in
Figure 2.

Figure 2: Single pipe example probability distribution functions
$q_{3}\sim\mathcal{U}[200,300]$ and $q_{3}\sim N(200,50/3)$ with truncated
tails.

We discretize the sample space $\Omega_{3}$ into $K=100$ stochastic cells that
uniformly partition the interval $[200,300]$ and solve the reformulated chance
constrained optimization problem (27) for several values of the acceptable
violation probability parameter $\epsilon_{3}$, namely,
$\epsilon_{3}\in\\{0.01,0.05,0.1\\}$. The distribution of the pressure at the
withdrawal node N3 is approximated by the kernel density estimate method
(`ksdensity` in MATLAB) using 10000 empirical samples of pressure distribution
obtained from the optimal solution. The pressure probability density function
(PDF) estimate at the withdrawal node is plotted for these constraint
violation probabilities and shown in Figure 3.

Figure 3: Pressure Probability Distribution at Withdrawal Node $J3$. Top:
Truncated Normal Uncertainty $q_{3}\sim N(200,50/3)$; Bottom: Uniform
Uncertainty $q_{3}\sim\mathcal{U}[200,300]$

Observe that the withdrawal pressure PDF retains a bell curve type shape for
normally distributed uncertain withdrawal flow but the pressure PDF for
uniform uncertain withdrawal has a non-uniform distribution over the pressure
values. The withdrawal pressure PDFs are consistent as they show higher
violation for higher values of $\epsilon$ and vice-versa. We also note that
the estimated probability of violation for the pressure constraint ($P\leq
P_{min}$) is higher in the uniform uncertainty case than the normal
distribution for same chance constraint violation probability parameter
($\epsilon_{3}$). This happens because the pressure drop inside the pipe
increases with increasing withdrawal flow resulting in lower pressure at the
withdrawal node, and the probability of withdrawal flow near the upper end of
the uncertainty interval $[200,300]$ is higher given uniform uncertainty than
when considering the normally distributed uncertainty scenario. A uniform
distribution for the uncertain flow also results in higher compressor ratio
solution $\alpha=\\{1.1985,1.192,1.187\\}$ compared to normal uncertainty
($\alpha=\\{1.183,1.171,1.164\\}$), where these enumerations correspond to the
values $\epsilon_{3}=\\{0.01,0.05,0.1\\}$ for the acceptable violation
probability.

Figure 4: 8-node Gas Pipeline Network System

###  IV-B 8-node Pipeline Test Network

In our second case study, we apply the SFV-based optimization approach to an
uncertainty management scenario for an 8-node gas pipeline test network with 5
pipes, 3 compressors, 2 withdrawal (at nodes J3 and J5), and 1 slack injection
node (J1) with constant pressure $P^{slack}_{1}=5.0$ MPa (shown in Fig.4).
This test network has been examined for steady-state and transient flow
algorithms in several previous studies [28, 19, 26, 22]. The lower bounds for
pressure at withdrawal nodes is specified as $P^{min}_{5}=4.0$ MPa at node J5,
and $P^{min}_{3}=3.0$ MPa at the other nodes. We suppose that the withdrawal
flow at node J5 is uncertain and parameterized using a uniformly distributed
uncertain parameter, whereas the withdrawal flow $q_{3}$ at node J3 is a
decision variable in the optimization problem. The uncertain withdrawal flow
at node J5 is modeled as $q_{5}(\omega)=d_{5}+r_{5}(\omega)$, where $d_{5}=64$
kg/s and $r_{j}=\mathcal{U}[0,32]$. In our implementation, $\Omega_{5}=[0,32]$
is discretized using $K=50$ uniform stochastic cells. We solve three cases
with different values of the withdrawal nomination constraint bound
$q_{3}^{\max}\equiv d_{3}^{\max}$, namely, $q_{3}^{\max}\equiv
d^{max}_{3}\in\\{200,300,\infty\\}$ kg/s, using the same value of violation
probability parameter $\epsilon=0.1$ in each case.

Estimates of the pressure PDFs at two withdrawal nodes J3 and J5 are shown in
Figure 5. The pressure profile shows no violation at withdrawal node J3 for
any instance. In contrast, the pressure at withdrawal node J5 exhibits bound
violation for all instances and the pressure profile for the instance with no
upper bound ($q_{3}^{max}=\infty$) has values less than $P^{min}_{5}=4.0$ MPa
for each uncertain scenario. The pressure PDFs also show a higher standard
deviation in withdrawal pressure at node J3 for $q^{max}_{3}=\infty$ and in
withdrawal pressure at node J5 for $q_{3}^{max}\in\\{200,300\\}$ kg/s.

Figure 5: Pressure probability distributions for 8-node case study. Top: Node
3; and Bottom: Node 5.

The optimal withdrawal flow $q_{3}$ at J3 is shown in Fig 6. We observe that
the withdrawal is at upper bound limits for $q_{3}^{\max}\equiv
d_{3}^{\max}\in\\{200,300\\}$ kg/s, but has a continuous PDF estimate when
there is no upper bound for the withdrawal flow variable.

Figure 6: Withdrawal Flow Probability Distribution at Node 3

###  IV-C Economic Interpretation of Lagrange Multipliers

We now explore the possibility of economic analysis of solutions to problem
(27) by use of the associated Lagrange multipliers to reconstruct stochastic
locational prices. In particular, by solving problem (27) we obtain the
Lagrange multipliers $\boldsymbol{\lambda}_{q}(\omega)$ that are associated
with the nodal mass flow balance constraint (27d), as well as Lagrange
multipliers $\boldsymbol{\lambda}_{d}(\omega)$ associated with upper bound
constraint (27i) for the optimized withdrawal flows the withdrawal bound
constraint (in particular, $d_{3}(\omega)\leq d^{max}_{3}$ the 8-node
example). These multipliers at the optimal solution, which we interpret as
marginal prices, follow the Karush-Kuhn Tucker (KKT) conditions, and can be
derived by assembling terms in the Lagrangian as shown in equation (28),
differentiating with respect to $\boldsymbol{d}$ to obtain the first order
condition in equation (29), and solving for the extremal value. For the 8-node
example, the relevant parts of the Lagrangian take the form

$\displaystyle
L(.,\boldsymbol{d})=-c_{3}^{d}(d_{3}+\mathbb{E}_{\omega}[r_{3}])+\sum_{\omega\in\Omega}\boldsymbol{\lambda}_{q}^{T}(\boldsymbol{A}\boldsymbol{\phi}(\omega)-\boldsymbol{q}(\omega))$
$\displaystyle\qquad\qquad+...+\sum_{\omega\in\Omega}\lambda_{d,3}(d_{3}-d^{max}_{3}),$
(28)

and taking the variation with respect to the optimized withdrawal $d_{3}$
yeilds

$\displaystyle\quad\nabla
L_{q_{3}}(.)=-\frac{c_{3}}{K}+\lambda_{q,3}+\lambda_{d,3}=0.$ (29)

Recall that the parameter $K$ denotes the number of stochastic cells used to
discretize the interval $\Omega$, which in this case is
$\Omega_{3}\equiv[200,300]$.

Figure 7: Marginal price probability distribution associated with the flow
balance constraint at Node J5

The probability distribution functions of $\lambda_{f,3}$ and $\lambda_{q,3}$
are shown in Figure 8, and in this case they are discrete. Observe that the
dual variables satisfy the KKT condition where $c_{3}=20$ and $K=50$ and the
sum of the dual variables is equal to $c_{3}/K=0.4$. We also observe that for
instances with $q^{max}_{3}\in\\{200,\infty\\}$ kg/s, the dual variable
$\lambda_{q,3}$ is constant at 0.02 and 0.39 respectively for each scenario,
whereas the dual variable $\lambda_{d,3}\approx$ 0 and 0.38 respectively. For
the instance with $q^{max}_{3}=300$ kg/s, the dual variables
($\lambda_{q,3},\lambda_{d,3}$) are dependent on the scenario, and their sum
is equal to 0.4 as expected based on (29). Based on these observations, we see
that the SFV-based computational method satisfies the KKT conditions for this
case study.

Figure 8: Marginal price probability distributions at Node 3. Top: associated
with the withdrawal bound constraint; Bottom: associated with the flow balance
constraint

Let us now examine the PDF estimate for locational marginal price
$\lambda_{q,5}$ at withdrawal node J5, which is shown in Figure 7. The
marginal price denotes the sensitivity of the optimal objective function value
$J(\alpha,d,s)$ to variation in the uncertain withdrawal flow $q_{5}(\omega)$
as $\lambda_{q,5}(\omega)=\partial J(\alpha,d,s)/\partial q_{5}(\omega)$.
Because $q_{5}(\omega)$ is a random variable, then we the multiplier
$\lambda_{q,5}(\omega)$ is a random variable as well, and its distribution can
be reconstructed from the dual variables corresponding to flow balance
constraints (27d) in the optimal solution to problem (27). We observe that
distribution of the marginal price $\lambda_{q,5}(\omega)$ of gas at node J5
is always zero for $q_{3}^{\max}=200$ and sometimes 3.18 when
$q_{3}^{\max}=300$. If there is no specified upper bound on optimized
withdrawal to node J3, i.e., $q_{3}^{\max}=\infty$, then the optimization
solution is driven by the probabilistic lower pressure consraint at node J5
and as a result the marginal price distribution at node J3 is continuous.

##  V Conclusion

We have developed a stochastic finite volume representation for solutions of
the steady-state gas flow equations on a pipeline network with compressors
subject to uncertain boundary conditions. This representation is used to
formulate a stochastic optimization problem for steady-state flow allocation
in a pipeline system subject to uncertain demands. This type of uncertainty
modeling enables the generalizable representation of arbitrary probability
distributions in problem parameters, and facilitates higher-order
reconstructions in computationally efficient calibration of probabilistic
inequality constraints. The primary contribution of our study is a chance-
constrained optimization formulation for steady-state flow allocation for a
gas pipeline system, and a method for computing solutions for compressor
configuration and nomination flows that optimize an operational and/or
economic objective in expectation. By extending the proven method of chance-
constraints using higher-order nonlinear reconstructions, our pipeline demand
uncertainty management approach provides probabilistic robustness guarantees
that can be calibrated based on risk tolerance. Furthermore, when the
objective function is interpreted as an expected economic value created for
users of the pipeline system, we show that the dual variables of the
optimization solution can be used to reconstruct distributions that we
interpret as probability densities of locational gas prices. Such solutions
could be used to quantify the financial volatility of energy in an
optimization-driven gas market mechanism that accounts for demand, supply,
transportation capacity, and price volatility. Compelling future work would
include analysis of convergence in modeling accuracy and probability measure,
as well as an investigation of the trade-offs between computational cost and
stochastic degrees of freedom examined. We note that, because consumption
uncertainty throughout a network may be correlated, it is possible to
represent wide-area uncertainty using only a few stochastic dimensions.

## References

  * [1] Jingwei Yang, Ning Zhang, Chongqing Kang, and Qing Xia.  Effect of natural gas flow dynamics in robust generation scheduling under wind uncertainty.  IEEE Transactions on Power Systems, 33(2):2087–2097, 2017.
  * [2] Line Roald and Göran Andersson.  Chance-constrained AC optimal power flow: Reformulations and efficient algorithms.  IEEE Transactions on Power Systems, 33(3):2906–2918, 2017.
  * [3] Bining Zhao, Antonio J Conejo, and Ramteen Sioshansi.  Unit commitment under gas-supply uncertainty and gas-price variability.  IEEE Transactions on Power Systems, 32(3):2394–2405, 2016.
  * [4] Reuven R. Levary and Burton V. Dean.  A natural gas flow model under uncertainty in demand.  Operations Research, 28(6):1360–1374, 1980.
  * [5] Nikolaos V. Sahinidis.  Optimization under uncertainty: state-of-the-art and opportunities.  Computers & chemical engineering, 28(6-7):971–983, 2004.
  * [6] Dimitris Bertsimas, David B. Brown, and Constantine Caramanis.  Theory and applications of robust optimization.  SIAM review, 53(3):464–501, 2011.
  * [7] Hesam Ahmadian Behrooz.  Managing demand uncertainty in natural gas transmission networks.  Journal of Natural Gas Science and Engineering, 34:100–111, 2016\.
  * [8] Victor M. Zavala.  Stochastic optimal control model for natural gas networks.  Computers & Chemical Engineering, 64:103–113, 2014.
  * [9] Kai Liu, Lorenz T. Biegler, Bingjian Zhang, and Qinglin Chen.  Dynamic optimization of natural gas pipeline networks with demand and composition uncertainty.  Chemical Engineering Science, 215:115449, 2020.
  * [10] Line A. Roald, Kaarthik Sundar, Anatoly Zlotnik, Sidhant Misra, and Göran Andersson.  An uncertainty management framework for integrated gas-electric energy systems.  Proceedings of the IEEE, 108(9):1518–1540, 2020.
  * [11] Saif R. Kazi, Kaarthik Sundar, Sidhant Misra, Svetlana Tokareva, and Anatoly Zlotnik.  Intertemporal uncertainty management in gas-electric energy systems using stochastic finite volumes.  In 2024 Power Systems Computation Conference (PSCC). IEEE, 2024\.
  * [12] Jiang-Bo Geng, Fu-Rui Chen, Qiang Ji, and Bing-Yue Liu.  Network connectedness between natural gas markets, uncertainty and stock markets.  Energy Economics, 95:105001, 2021.
  * [13] Jiqian Wang, Feng Ma, Elie Bouri, and Juandan Zhong.  Volatility of clean energy and natural gas, uncertainty indices, and global economic conditions.  Energy Economics, 108:105904, 2022.
  * [14] Fred C. Schweppe, Michael C. Caramanis, Richard D. Tabors, and Roger E. Bohn.  Spot pricing of electricity.  Springer Science & Business Media, 2013.
  * [15] W. Pepper, B. J. Ring, E. G. Read, and S. R. Starkey.  Implementation of a scheduling and pricing model for natural gas.  Handbook of Networks in Power Systems II, pages 3–35, 2012.
  * [16] Aleksandr M. Rudkevich and Anatoly Zlotnik.  Locational marginal pricing of natural gas subject to engineering constraints.  In Proc. of the 50th Hawaii International Conference on System Sciences, pages 3092–3101, 2017.
  * [17] Anatoly Zlotnik, Kaarthik Sundar, Aleksandr M Rudkevich, Aleksandr Beylin, and Xindi Li.  Optimal control for scheduling and pricing intra-day natural gas transport on pipeline networks.  In 2019 IEEE 58th Conference on Decision and Control (CDC), pages 4887–4884. IEEE, 2019.
  * [18] Anatoly Zlotnik, Aleksandr M. Rudkevich, Evgeniy Goldis, Pablo A. Ruiz, Michael Caramanis, Richard Carter, Scott Backhaus, Richard Tabors, Richard Hornby, and Daniel Baldwin.  Economic optimization of intra-day gas pipeline flow schedules using transient flow models.  In PSIG Annual Meeting, pages PSIG–1715. PSIG, 2017.
  * [19] Svetlana Tokareva, Anatoly Zlotnik, and Vitaliy Gyrya.  Stochastic finite volume method for uncertainty quantification of transient flow in gas pipeline networks.  Applied Mathematical Modelling, 125:66–84, 2024.
  * [20] C.-S. Huang, C.-H. Hung, and Song Wang.  A fitted finite volume method for the valuation of options on assets with stochastic volatilities.  Computing, 77:297–320, 2006.
  * [21] Michael Herty, Jan Mohring, and Veronika Sachers.  A new model for gas flow in pipe networks.  Mathematical Methods in the Applied Sciences, 33(7):845–855, 2010\.
  * [22] Saif R. Kazi, Kaarthik Sundar, Shriram Srinivasan, and Anatoly Zlotnik.  Modeling and optimization of steady flow of natural gas and hydrogen mixtures in pipeline networks.  International Journal of Hydrogen Energy, 54:14–24, 2024.
  * [23] E. S. Menon.  Gas pipeline hydraulics.  CRC Press, 2005.
  * [24] Marc Vuffray, Sidhant Misra, and Michael Chertkov.  Monotonicity of dissipative flow networks renders robust maximum profit problem tractable: General analysis and application to natural gas flows.  In 2015 54th IEEE Conference on Decision and Control (CDC), pages 4571–4578. IEEE, 2015.
  * [25] Fei Wu, Harsha Nagarajan, Anatoly Zlotnik, Ramteen Sioshansi, and Aleksandr M Rudkevich.  Adaptive convex relaxations for gas pipeline network optimization.  In 2017 American Control Conference (ACC), pages 4710–4716. IEEE, 2017.
  * [26] Mo Sodwatana, Saif R. Kazi, Kaarthik Sundar, and Anatoly Zlotnik.  Optimization of hydrogen blending in natural gas networks for carbon emissions reduction.  In 2023 American Control Conference (ACC), pages 1229–1236. IEEE, 2023.
  * [27] Sidhant Misra, Marc Vuffray, and Anatoly Zlotnik.  Monotonicity properties of physical network flows and application to robust optimal allocation.  Proceedings of the IEEE, 108(9):1558–1579, 2020.
  * [28] Vitaliy Gyrya and Anatoly Zlotnik.  An explicit staggered-grid method for numerical simulation of large-scale natural gas pipeline networks.  Applied Mathematical Modelling, 65:34–51, 2019.
  * [29] Iain Dunning, Joey Huchette, and Miles Lubin.  Jump: A modeling language for mathematical optimization.  SIAM Review, 59(2):295–320, 2017.
  * [30] Andreas Wächter and Lorenz T. Biegler.  On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming.  Mathematical Programming, 106(1):25–57, 2006.
  * [31] Shriram Srinivasan, Kaarthik Sundar, Vitaliy Gyrya, and Anatoly Zlotnik.  Numerical solution of the steady-state network flow equations for a non-ideal gas.  IEEE Transactions on Control of Network Systems, 2022.


# A Longitudinal Analysis of Online Ad-Blocking Blacklists

Saad Sajid Hashmi  Macquarie University
<EMAIL_ADDRESS>Muhammad Ikram  Macquarie University
University of Michigan
<EMAIL_ADDRESS>Mohamed Ali Kaafar  Macquarie University
Data61, CSIRO
<EMAIL_ADDRESS>
###### Abstract

Websites employ third-party ads and tracking services leveraging cookies and
JavaScript code, to deliver ads and track users’ behavior, causing privacy
concerns. To limit online tracking and block advertisements, several ad-
blocking (black) lists have been curated consisting of URLs and domains of
well-known ads and tracking services. Using Internet Archive’s Wayback Machine
in this paper, we collect a retrospective view of the Web to analyze the
evolution of ads and tracking services and evaluate the effectiveness of ad-
blocking blacklists. We propose metrics to capture the efficacy of ad-blocking
blacklists to investigate whether these blacklists have been reactive or
proactive in tackling the online ad and tracking services. We introduce a
stability metric to measure the temporal changes in ads and tracking domains
blocked by ad-blocking blacklists and a diversity metric to measure the ratio
of new ads and tracking domains detected. We observe that ads and tracking
domains in websites change over time, and among the ad-blocking blacklists
that we investigated, our analysis reveals that some blacklists were more
informed with the existence of ads and tracking domains, but their rate of
change was slower than other blacklists. Our analysis also shows that Alexa
top 5K websites in the US, Canada, and the UK have the most number of ads and
tracking domains per website, and have the highest  proactive scores. This
suggests that ad-blocking blacklists are updated by prioritizing ads and
tracking domains reported in the popular websites from these countries.

##  I Introduction

For almost as long as the commercial world wide web has existed, third-party
advertisements (in short ads) and web tracking practices have been used to
support free services [36]. To this end, first-party websites leverage various
techniques such as third-party iframes and JavaScript to show ads and track
users’ activities on websites. While this helps to support the business model
of most content providers, its use in showing ads and profiling users’
activities raises serious privacy concerns. Besides, many third-party services
such as ads are being used as malvertisements [34, 8].

A number of studies have characterized and measured the ads and tracking
ecosystem to protect user privacy and limit intrusive ads, resulting in
numerous ad-blocking tools such as Ghostery [4], Adblock [1], Adblock Plus
[2], Disconnect [7], and Privacy Badger [6] for the web [32] and mobile
platforms [33]. Most of these tools are fueled by community-driven public
blacklists (such as EasyPrivacy [16], EasyList [14], FanboyList [17], and
hpHosts [18]) and are evaluated for their (in)effectiveness [44, 32]. Although
the research work in this domain has been on-going for some time and ad-
blocking tools have been developed, we believe these studies typically consist
of short-term measurements of specific tracking techniques and the ad-blocking
tools have not been comprehensively evaluated over time. Given that the online
ads, web tracking, and prevention landscape is continuously evolving, studies
performed at a snapshot or longitudinally starting in the present may not
comprehensively illuminate online ads/tracking to improve ad-blocking.

To fill the gap, we leverage the historical data extracted from the Internet
Archive’s Wayback Machine [5] to analyze online ads (resp. web tracking) and
ad-blocking blacklists (listed in Table I). From online ads and web tracking
perspective, we look at the prevalence and the evolution of ads and tracking
domains since 2009. From an ad-blocking and tracking prevention perspective,
we measure the longitudinal performance of various ad-blocking blacklists with
respect to their rate of change and update speed. This paper makes the
following main contributions:

Prominence of a&t domains. We study the evolution of the number of ad and
tracking (a&t for short) domains embedded in websites over time. Globally,
since 2010 the presence of a&t domains has been dynamic on an average of 55%
of the websites crawled, with 45% exhibiting an annual increase and 10% of the
websites exhibiting an annual decrease in the number of a&t domains embedded.
We observe a bias in the presence of a&t domains across different countries.
with the Alexa top-5K websites in the US, Canada, and the UK exhibiting the
highest number of a&t domains per website, whereas China and Iran had the
least. We also examine the number of a&t services detected by each of the
analyzed ad-blocking blacklist (see Table I). We observe that the blacklists
such as Mahakala and hpHosts report the highest number of a&t services while
Cybercrime and Sa-blacklist detect the least.

Rate of change of tracking prevention lists. To measure the churn of third-
party a&t domains within a blacklist, we introduce two new metrics: stability
and diversity. Although each year, the ad-blocking blacklists mostly detect
a&t domains that have already been reported in the preceding years, they
continue to detect new third-party a&t services in the Alexa top 5K websites.
We observe a gradual drop of new a&t services detected by ad-blocking
blacklists: on average from 19% in 2011 to 10% in 2017.

Update speed of ad-blocking blacklists. We introduce two new metrics that
measure the speed of updates of ad-blocking blacklists by (de)listing third-
party a&t domains. Most of the ad-blocking blacklists have proactively
detected a&t services since 2009. Among the analyzed ad-blocking blacklists,
we found that AdZHosts has the highest average reactive score of 0.6 and
Mahakala has the highest average proactive score of 0.49.

##  II Ads and Tracking on Web

We briefly discuss the mechanism of third-party online a&t services and
overview ad-blocking blacklists.

Online ads and web tracking: Many online businesses and their revenue depend
on online advertisements delivered by third-parties. With the advent of new
technologies and platforms, third-party online a&t services (resp. websites’
publishers) have adopted new ways to deliver ads and track users’ interaction.
In Web platform, webpages isolate advertisements from content by placing ads
in an iframe, where ad and tracking content hosted in the iframe is isolated
from the hosting webpage, and browsers allow only specific cross frame
interactions [27]. The website publisher includes the ad and tracking
components (i.e., JavaScript codes and cookies belonging to a third-party ad
and tracking service) in their websites and allocate iframes to display ads.
The website publisher is paid by the ad provider based on the number of ad
clicks or impressions, or both. Most often, these techniques track users
across websites for targeted ads—delivering ads based on their browsing/search
or purchase history—revealing sensitive information about website users.

Ads and tracking services operate by storing the user’s personal identifiable
information (PII) in the form of browser cookies, HTML5 local storage and
flash cookies which is then shared with (other) third-parties either
implicitly via HTTP referrer or explicitly through tracker-provided JavaScript
programs. If the user is tracked within the website, only by the first party
domain (website directly visited by the user), it is termed as first-party
tracking. An application of first-party tracking is to count the repeated
visits of a user to the website for analytics. If the user is tracked on the
website by another domain, this is termed as third-party tracking. For
example, if the user visits bbc.com, and the advertisement on bbc.com is
fetched from advertiser.com, then the user’s browser sends a request to
advertiser.com without the user being aware of it. In our study, we focus on
the JavaScript programs containing third-party domains which are embedded in
the websites that are popular according to Alexa websites ranking. Therefore,
in this paper, the terms ‘ads and trackers’ and ‘third-party ads and trackers’
are used interchangeably to mean third-party ads and tracking domains.

Ad-blocking blacklists: To counter the privacy threat posed by third-party a&t
services, several ad-blocking tools have been developed. They include pop-up
blockers, privacy preserving web proxies (like Privoxy [23], Proxomitron [24],
and Pi-hole [22]) to filter and block a&t traffic at network-layer and ad-
blocking browser extensions (e.g., Adblock Plus, Ghostery, and Disconnect)
that prevent web resources from rendering in the browser [32, 33]. The most
popular of these tools are the ad-blocking browser extensions which rely on
filter-lists (also referred to as tracking prevention lists or ad-blocking
lists), comprising of sets of URLs and domains of a&t services. These
blacklists are either crowd-sourced with feedback from web users or maintained
by ad-blocking tools developers. Table I overviews several well-known a&t
blacklists used by ad-blocking tools.

##  III Dataset and Methodology

In this section, we provide an overview of our dataset and our measurement
methodology.

###  III-A Dataset

We aim to study the prevalence of a&t services in popular websites worldwide.
We choose two sets of popular websites: (i) Alexa top 5K global websites (5K);
and (ii) Alexa top 5K websites in the following fourteen countries: Australia
(AU), Brazil (BR), Canada (CA), China (CN), Germany (DE), Israel (IL), India
(IN), Iran (IR), Russia (RU), Saudi Arabia (SA), Singapore (SG), Ukraine (UA),
United Kingdom (UK), and United States (US).

To conduct our longitudinal study, we use web data from the Internet Archive’s
Wayback Machine. Since 1996, the Wayback Machine has archived full websites,
including JavaScript codes, style sheets, and any multimedia resources that it
can identify statically from the site’s content. We refer to a single capture
of a webpage (resp. a&t domains in an ad-blocking blacklist) as a snapshot.
Wayback Machine mirrors past snapshots of these websites on its own servers.
We scrape the HTML DOM (Document Object Model, where HTML elements are defined
as objects) of Alexa top 5K global and country-wise websites. The purpose of
scraping the HTML DOM is to extract all JavaScript codes (along with their
sources) from the websites, which is then used for our analysis. Wayback
Machine has a number of archived snapshots, with varying intervals, for each
website. We utilize Memento API [20] to capture snapshots at intervals of
three months to detect most—if not all—of the ads and tracking domains that
appeared on a website for a given year, during the period of 2009 to 2017.
Memento API provides the nearest time-stamp for the archived snapshot of a
website (resp. an ad-blocking blacklist) from the date provided.

We empirically observed that 85% of the snapshots are within an interval of 6
months. From these snapshots, we then obtain third-party sources from the
embedded JavaScript codes, and by comparing it with the second level domains
of ad-blocking blacklists, we obtain the a&t domains in the crawled websites.
For example, if http://ad.doubleclick.net/dot.gif appears in the script tag of
sportsbet.com, then the second level domain (doubleclick.net) is extracted
from the URL, and compared against the domains in the ad-blocking blacklists.

Similarly, by using Internet Archive’s Wayback Machine, we obtain the
snapshots of ad-blocking blacklists. Table I lists the number of distinct
domains and the snapshots period of the analyzed ad-blocking blacklists.

TABLE I: Overview of the analyzed ad-blocking blacklists. |  | Snapshots | # Distinct
---|---|---|---
Blacklist | Description | Period | Domains
AdZHosts [11] | Ad/tracker/malicious hosts | 2017-17 | 81,294
AdGuard [10] | Ad/tracker/phishing hosts | 2017-17 | 12,881
Cameleon [12] | Ad/tracker hosts | 2009-16 | 47,390
CyberCrime [13] | Tracker/malware/phishing hosts | 2013-16 | 6,557
EasyList [14] | Ad/trackers/analytics | 2010-17 | 13,882
FanboySocial [17] | Social buttons/widgets hosts | 2013-17 | 282
EasyPrivacy [16] | Trackers/analytics hosts | 2011-17 | 7,050
QuidsUp [25] | Tracker/malicious hosts | 2017-17 | 865
EasyList_China [15] | Chinese EasyList websites | 2014-17 | 18,639
hpHosts [18] | Ad/tracker/malicious hosts | 2009-17 | 856,498
MVPs [21] | Ad/tracker/malicious hosts | 2011-17 | 414,461
AdAway [9] | Mobile ad/tracker hosts | 2013-16 | 137
Disconnect [7] | Ad/tracker hosts | 2015-16 | 15,238
Mahakala [19] | Ad/tracker hosts | 2015-17 | 1,846,198
Sa-blacklist [26] | Ad/tracker hosts | 2009-13 | 1,253,243

###  III-B Methodology

The objective of this work is three fold. We study the: (i) prominence of a&t
domains on the Alexa top 5K worldwide as well as top 5K websites in different
countries, (ii) rate of change of ad-blocking blacklists, and (iii) update
speed of ad-blocking blacklists.

Prominence of a&t domains: For our first objective, we compare the second-
level domains of JavaScript sources (observed in retrospectively crawled top
5K websites from the years 2009 to 2017) with those of known a&t domains from
the ad-blocking blacklists. This comparison reveals the number of a&t domains
blocked in different ad-blocking blacklists’ yearly snapshots and the number
of websites that embed them.

By measuring the number of a&t domains on the Alexa top 5K websites of
different countries, we also measure the prevalence of a&t domains in a given
country’s top websites. Besides detecting existing a&t domains in the wild, an
effective ad-blocking blacklist also has to evolve and detect new a&t domains
over time to counter the ever-changing tracking landscape. Therefore, we also
investigate the number of new a&t domains detected by the ad-blocking
blacklists from 2009 to 2017.

Rate of change of ad-blocking blacklists: For our second objective, we compare
the a&t domains detected in two consecutive years. This comparison reveals
whether there has been any change in the a&t domains blocked compared with the
preceding year, and indicates the evolution of a&t domains in the ad-blocking
blacklist. To this end, we propose the stability metric:

$Stability=\frac{\mid T_{i}\cap T_{i-1}\mid}{\mid T_{i}\cup T_{i-1}\mid},$ (1)

where ${\bf T}_{i}$ and ${\bf T}_{i-1}$ are the sets of distinct a&t domains
detected in the current year, $i$, and the preceding year, $i-1$. Stability
score of 1 implies that the list of a&t domains blocked in a given year is
identical to the list from the preceding year, and stability score of 0
implies that the list of a&t domains blocked in a given year is completely
different than the list from the preceding year. By default, the stability of
any ad-blocking blacklist for the first year is always zero.

From Eq. 1, the stability score of an ad-blocking blacklist may be low,
however, it can detect new a&t domains, because the stability score considers
the given year and the preceding year only. Therefore, if a&t domains that
were blocked in years before the preceding year are detected again, they have
no impact on the measurement of stability. To overcome this limitation, we
introduce the diversity metric, which measures the ratio of new a&t domains
detected by the ad-blocking blacklist for a given year. The diversity of an
ad-blocking blacklist for a given year is measured as:

$Diversity=1-\frac{\mid T_{i}\cap T_{a}\mid}{\mid T_{i}\mid},$ (2)

where ${\bf T}_{a}$ is the set of distinct a&t domains detected in all
previous years. So a diversity of $1$ implies that all the a&t domains blocked
in a given year are new, i.e., they have not been detected by that ad-blocking
blacklist in the previous years, and a diversity of 0 implies that none of the
a&t domains blocked in a given year are new. By default, the diversity score
for any ad-blocking blacklist is always $1$ for the first year.

Update speed of ad-blocking blacklists: For our third objective, we study the
speed by which each blacklist updated itself to encounter new threats. The
time difference is taken to be the number of days between the first occurrence
of an a&t domain in a website in the given year, and its addition to the
blacklist. The time difference is taken to be positive if an a&t domain is
seen on the web before it gets added to the blacklist. It is negative
otherwise. We consider the blacklist is reactive if the time difference is
positive, and proactive otherwise. We aim to introduce the metrics for which a
high score implies high performance. The performance of a blacklist is
considered good if it detects more a&t domains in the least amount of time or
it contains a&t domains prior to them being observed on the web. Therefore, to
measure the reactive score, we divide the number of a&t domains per website
detected reactively by the average time difference, and compute proactive
score by taking the product of the number of a&t domains per website detected
proactively and average time difference. The scores are then normalized
between $0$ and $1$. The reactive and proactive scores of an ad-blocking
blacklist for a given year are measured as:

$Reactive=\;\frac{(\frac{\sum^{domains}\mid
T_{react}\mid}{t_{days}})}{Reactive_{max}}$ (3)
$Proactive=\;\frac{(\sum^{domains}\mid T_{proact}\mid\;)\;*\;\mid
t_{days}\mid}{Proactive_{max}}$ (4)

where $T_{react}$ and $T_{proact}$ are the set of distinct a&t domains per
website detected reactively and proactively, respectively, in the given year
and $t_{days}$ is the average time difference in days to add those a&t domains
in blacklists. Since for proactive a&t domains blocked the average time
difference is negative, we only consider the magnitude and ignore the negative
sign.

A low reactive (resp. proactive) score implies that the ad-blocking blacklist
is either blocking low number of a&t domains or it is taking a large amount of
time to add a&t domains in the blacklist (or in case of proactive score, adds
a&t domains just recently before being observed in the top 5K websites) or
both. A high reactive (resp. proactive) score implies that the ad-blocking
blacklist is either blocking large quantities of a&t domains or it is taking
less amount of time to add those a&t domains in the list (or in case of
proactive score, adds a&t domains prior to when they are observed in the top
5K websites) or both.

Note that if an ad-blocking blacklist started its service (or we capture its
first snapshot from Wayback Machine) after the first occurrence of a given a&t
domain on the web, the reactive score is computed by taking the time
difference from when the a&t domains are added to a given ad-blocking
blacklist. If an a&t domain appears in the first version (or snapshot) of a
given blacklist then we consider the time difference as $1$ instead of $0$ (to
avoid division by $0$).

Similarly, for a given ad-blocking blacklist, we compute the reactive and
proactive scores per country per year. We then average the scores over the
countries and years. The normalized reactive and proactive scores are then
computed by dividing the score by the maximum observed score at the lowest
aggregation level i.e., total a&t domains per year per country.

In Section IV-B, we use the above four metrics for different countries. All
these metrics take the set of a&t domains as input. For the ad-blocking
blacklists, that set of a&t domains comprises the distinct a&t domains blocked
by the given ad-blocking blacklist. For countries, the set of a&t domains
consists of the distinct a&t domains that occurred in ad-blocking blacklists
on the Alexa top 5K websites of the given country. For example, the stability
of a country measures the similarity in the a&t domains detected (by the union
of all the ad-blocking blacklists) on the Alexa top 5K websites of that
country compared with the preceding year. Similarly, we use the other three
metrics for the countries.

##  IV Analysis and Results

Next, we present results of our analysis. We measure the prominence of a&t
domains to establish whether some countries’ top 5K websites consistently had
high number of a&t domains over the period of time, and study whether websites
added new a&t domains each year to avoid ad-blocking tools. In order to
understand the variation in effectiveness of ad-blocking blacklists in recent
years, we investigate how the ad-blocking blacklists have performed in recent
years in terms of number of a&t domains detected on the Alexa top 5K websites,
and their rate of change measured by stability and d͡iversity metrics. We also
investigate the update speed of ad-blocking blacklists in terms of reactive
and proactive scores.

###  IV-A Prominence of Ads and Tracking Domains

Annual Measurement of a&t domains: From the crawled data, we measure the
number of a&t services observed in each year. Figure 1 shows the number of a&t
domains blocked by ad-blocking blacklists. We observe that a&t domains peak in
2011, before gradually declining to the mean number 317,357 of a&t domains for
all the years. One possible explanation for the decline after 2011 is that
some websites embed new a&t domains over time either because the a&t domain is
not functioning anymore or to evade ad-blocking tools [45]. Also, all the new
a&t domains may not have been detected yet by the analyzed ad-blocking
blacklists, and thus not reflected in the annual a&t domains measurement. For
example, we note that in 2011, wikia.com embed google-analytics.com,
quantserve.com, and nocookie.net in its HTML whereas in 2016 wikia.com had
only optimizely.com as a third-party a&t domain.

Country-wise measurement: Next, we map the a&t domains observed on the Alexa
top 5K websites of different countries. We acknowledge that there are some
limitations in capturing snapshots with Wayback Machine, and many websites
from top 5K had to be discarded due to the challenges mentioned in §VI.
Therefore, the number of websites do not remain equal for all the countries,
and instead of measuring the total number of a&t domains in all the websites,
we measure the average number of a&t domains per website for each country.

(a)

(b)

(c)

(d)

Figure 1: (a) Annual sum of distinct a&t domains per website (b) Distinct a&t
domains per website on the Alexa top 5K websites of each country (c) Annual
Percentage of websites with an increase or decrease in number of a&t domains
from preceding year (d) Distribution of change in number of a&t domains from
the websites in Figure 1.

Figure 1 shows the union of distinct a&t domains per website from 2009 to
2017. We observe that the US, Canada, and the UK have the most number of a&t
domains per website. As the ad-blocking blacklists are mostly updated through
crowd-sourcing [43], the prevalence of a&t domains across countries also
suggests that the crowd-sourcing users for the analyzed ad-blocking blacklists
are more active in the countries reporting the most number of a&t domains,
compared to countries like Iran and China, that have detected the lowest
number of a&t domains per website.

Figure 4 shows the annual number of a&t domains blocked on the Alexa top 5K
websites of the top seven countries in Figure 1. Note that for ease of
presentation, we select countries with most number of a&t domains detected on
the Alexa top 5K websites. The a&t domains country-wise trend is similar to
the overall trend in Figure 1, suggesting that the new a&t domains may not yet
have been detected by the analyzed ad-blocking blacklists.

Domain-wise measurement of a&t domains: Figure 1 shows the percentage of
websites with a change in the number of a&t domains from the preceding year.
We observe that, from 2010 to 2017, the presence of a&t domains has been
dynamic on an average of 55% of the analyzed websites, with 45% exhibiting an
annual increase and 10% of the websites exhibiting an annual decrease.
Moreover, nearly 80% of those changes are a difference of one or two a&t
domains (cf. Figure 1). For example, in tomshardware.com there are four a&t
domains in 2009, ten in 2010, six in 2011, one in 2012, five in 2013, two in
2014, none in 2015, two in 2016, and one in 2017.

Measurement of a&t domains in blacklists: Table II lists the number of a&t
domains detected by the ad-blocking blacklists, and the mean time (in days)
taken to add them in blacklist. The table also lists the number of a&t domains
detected reactively and proactively for each ad-blocking blacklist. The number
of a&t domains detected is the annual sum of distinct a&t domains per website
for the years 2009 to 2017. We observe (cf. Table II) that the number 2,460K
of a&t domains blocked by Mahakala is the highest. Not surprisingly, Mahakala
also contains the most number 1,846K of distinct a&t domains We also observe
that EasyList and EasyPrivacy contain 13,882 and 7,050 a&t domains,
respectively, and these two blacklists are able to detect 393K instances of
those a&t domains within the analyzed top 5K websites. Surprisingly, Sa-
blacklist that contain 1,253K a&t domains is only able to detect 65,698
instances of those a&t domains. This reveals that EasyList and EasyPrivacy
though fairly lean contain mostly useful rules to block a&t domains in the top
5K websites. Conversely Sa-blacklist contains mostly stale or otherwise
useless rules to block a&t domains on the Alexa top 5K websites. CyberCrime
blocks the least number 3,349 of a&t domains, and one reason why this number
is low is that a significant fraction of domains in CyberCrime belongs to
malvertisements [3].

Table II lists ad-blocking blacklists ranked by the number of a&t domains
blocked on the Alexa top 5K websites. Annual number of a&t domains blocked by
the top seven ad-blocking blacklists in figure 3 indicates that the prevalence
of a&t domains peaks in 2011 for all the lists and then gradually declines and
remains almost unaltered over the remaining period. This trend is similar to
the number of a&t domains observed in Figure 1.

(a)

(b)

Figure 2: (a) Average stability and diversity score of ad-blocking blacklists
(b) Distribution of time taken (in days) to add a&t domains in each blacklist.
The blacklists on left side are proactive whereas those on right side are
reactive.

(a)

(b)

(c)

Figure 3: (a) Annual a&t domains detected by blacklists (b) Stability and (c) Diversity of blacklists. With the gradual rise in stability and gradual decline in diversity, blacklists are blocking more a&t domains similar to the preceding years, and less ratio of new a&t domains. TABLE II: Number of a&t domains detected and time taken (in days) to add them in ad-blocking blacklists. The number of a&t domains detected is the annual sum of distinct a&t domains per website for the years 2009 to 2017. Negative values for days indicate that the a&t domains were added in blacklists before being observed on the Alexa top 5K websites. Black- | Total a&t | Mean | # Reactive | Mean | # Proactive | Mean
---|---|---|---|---|---|---
List | Domains | Days | a&t Domains | Days | a&t Domains | Days
Mahakala | 2,460,303 | -70 | 1,915,941 | 35 | 544,362 | -438
hpHosts | 662,287 | -289 | 240,112 | 1,220 | 422,175 | -1,148
MVPs | 642,023 | 294 | 399,207 | 1,093 | 242,816 | -1,020
AdZHosts | 566,228 | -5 | 509,887 | 1 | 56,341 | -46
EasyList_China | 423,514 | -78 | 290,945 | 135 | 132,569 | -545
EasyList | 393,990 | -381 | 150,384 | 815 | 243,606 | -1,120
EasyPrivacy | 393,914 | -510 | 143,316 | 321 | 250,598 | -985
AdGuard | 386,220 | -3 | 370,247 | 1 | 15,973 | -80
Cameleon | 365,851 | -1,159 | 28,337 | 489 | 337,514 | -1,297
Disconnect | 199,773 | -93 | 155,389 | 7 | 44,384 | -441
FanboySocial | 135,163 | -283 | 76,398 | 53 | 58,765 | -719
QuidsUp | 116,297 | -8 | 101,852 | 1 | 14,445 | -68
AdAway | 79,861 | -290 | 47,370 | 25 | 32,491 | -749
Sa-blacklist | 65,598 | -810 | 14,117 | 546 | 51,481 | -1,182
CyberCrime | 3,349 | 10 | 3,214 | 22 | 135 | -257

###  IV-B Rate of Change

Stability and diversity analysis of blacklists: Figure 2 shows the average
stability and diversity scores of the ad-blocking blacklists in detecting a&t
domains in the top 5K global as well as the top 5K websites in different
countries from 2009 to 2017. We observe that the ad-blocking blacklists having
high stability scores tend to have low diversity scores and vice-versa. For
instance, AdAway and Fanboy Social have highest average stability and lowest
average diversity, whereas CyberCrime and AdGuard have lowest average
stability and highest average diversity. The average stability scores of top
seven ad-blocking blacklists (in Table II) range between 0.5 and 0.6, and the
average diversity scores are in the range of 0.2 to 0.3.

The annual average stability and diversity scores of the ad-blocking
blacklists in Figure 3 and 3 show the gradual rise in stability and the
gradual decline in diversity of ad-blocking blacklists, respectively. The
stability is $0$ and diversity is $1$ in 2009. In 2010, we observe that the
set of distinct a&t domains blocked by the ad-blocking blacklists is
significantly different from the one in 2009. This is evident from the
stability and diversity scores, which are at their lowest and highest peaks,
respectively. Since the diversity score is less than 0.5, it implies that the
set of distinct a&t domains blocked in 2010 had less than half the a&t domains
that were newly detected in 2010. We also observe a sharp increase in
stability and steep decline in diversity in 2011, and no significant variation
in the later years. This suggests that although the annual set of a&t domains
blocked by each ad-blocking blacklist is significantly similar (at least 65%
in 2017) to the annual set of a&t domains from the preceding years, new a&t
domains are still appearing on the Alexa top 5K websites.

Stability and diversity analysis of countries: Figure 5 shows the average
stability and diversity scores of different countries based on the a&t domains
blocked on the Alexa top 5K websites. The score for each country is computed
by measuring the stability and diversity for each year and then taking the
average over the years. We observe that Singapore, UK, and US have the lowest
diversity and conversely the highest stability scores, whereas Ukraine,
Israel, and Saudi Arabia have the highest diversity and conversely the lowest
stability scores.

Figures 4 and 4 show the annual stability and diversity scores, respectively,
of the top seven countries (in Figure 1). Although the stability scores are
fluctuating sharply for some countries, we notice that the diversity scores
remain stable from 2011 onwards. This suggests that the percentage of new a&t
domains blocked on the Alexa top 5K websites in these countries drops slightly
from 20% in 2011 to 15% in 2017 at the most.

(a)

(b)

(c)

Figure 4: (a) Annual a&t domains per website blocked on the Alexa top 5K of
countries (b) Average stability and (c) Average diversity scores of top seven
countries.

(a)

(b)

Figure 5: (a) Average stability and diversity (b) Average reactive and
proactive scores of countries. Countries with high proactive scores also tend
to have a high number of a&t domains on their Alexa top 5K websites (cf.
Figure 1).

###  IV-C Update Speed

Reactive and proactive analysis of blacklists: Figure 2 shows the distribution
of time taken in days to add a&t domains observed on the Alexa top 5K
websites, in the ad-blocking blacklists. In the distribution, we have omitted
the a&t domains that were observed on the Alexa top 5K websites before the ad-
blocking blacklists’ first snapshot. We notice that since we started capturing
the a&t domains in the ad-blocking blacklists (or the ad-blocking blacklists
started their service), barring MVPs and CyberCrime, the ad-blocking
blacklists have reported most of the a&t domains proactively.

Next, we measure the update speed of ad-blocking blacklists with reactive and
proactive metrics. Figure 6 shows the average reactive and proactive scores
for the ad-blocking blacklists. Mahakala has the highest average proactive
score whereas AdZHosts has the highest average reactive score. We note that
among the top seven blacklists in Table II, six are in the top seven
blacklists by proactive score.

Figure 6 shows the annual reactive scores for the ad-blocking blacklists. All
the a&t domains blocked by AdZHosts from 2009 to 2016 were included in the
first captured snapshot of AdZHosts. Therefore, the time difference was taken
as $1$ day which resulted in a very high reactive score. The rise in the
reactive score of AdZHosts from 2009 to 2011 is due to the rise in the number
of a&t domains blocked by AdZHosts in that period (see Figure 3). Mahakala’s
reactive score peaks in 2012, and then declines gradually from 2013 to 2015
because of an increase in the average time difference. Furthermore, in 2016,
the reactive score of Mahakala decreases sharply because of the low number of
reactive a&t domains identified.

In Figure 6 we notice a gradual increase in proactive score from the year in
which we first captured the snapshot of the ad-blocking blacklist, because of
an increase in the average time difference with each subsequent year. The
proactive score does not rise significantly for AdZHosts in 2017 because the
average time difference is least compared to the other ad-blocking blacklists.
The sharp rise in Mahakala’s proactive score in 2016 and 2017 is due to
Mahakala blocking the most number of a&t domains in the analyzed top 5K
websites. This suggests that since ad-blocking blacklists started their
service, their proactive scores have been increasing gradually at nearly the
same rate, apart from Mahakala which has risen sharply since 2015. Therefore,
based on the number of a&t domains blocked and the reactive and proactive
score, we conclude that in the last two years, Mahakala has been the most
effective ad-blocking blacklist. The number of black-listed domains in
Mahakala are also the highest, and if we are applying ad-blocking blacklist on
a performance constrained device then a leaner ad-blocking blacklist is most
effective. In that case, EasyList, EasyList_China and EasyPrivacy are the
recommended ad-blocking blacklists.

Reactive and proactive analysis of countries: Figure 5 illustrates the average
reactive and proactive scores of different countries from 2009 to 2017. US has
the highest reactive and proactive scores, whereas Iran has the lowest
reactive and proactive scores. Out of the top seven countries ranked by the
number of a&t domains per website (in Figure 1), six have the highest
proactive scores. Furthermore, among the bottom five countries ranked by a&t
domains per website, four are the bottom most countries by proactive scores.

(a)

(b)

(c)

Figure 6: (a) Average reactive and proactive scores of blacklists (b) Reactive
and (c) Proactive scores of top seven blacklists.

(a)

(b)

Figure 7: (a) Reactive and (b) Proactive scores of countries. The number of
a&t domains detected reactively has decreased in recent years, whereas those
identified proactively has risen.

Figures 7 and 7 show the annual reactive and proactive scores, respectively,
of the top seven countries by the number of a&t domains per website. We notice
a gradual increase in the reactive score of countries from 2009 to 2015. The
reactive score depends on the number of distinct a&t domains per website
detected and the average time difference in days (cf Eq. 3). We notice a steep
increase in reactive scores in 2016 because of less average time difference.
Then in 2017, we observe a steep decline in reactive scores. Although the
average time difference is less in 2017, the number of distinct a&t domains
per website detected is also less. In Figure 7 we observe a gradual rise in
proactive score from 2009 to 2015 due to an increase in average time
difference (cf. Eq. 4). Though the time difference increases steadily in 2016
and 2017, the number of a&t domains blocked proactively in 2016 and 2017 is
significantly greater than the previous years. Therefore, we observe a steep
increase in proactive scores from 2016 onwards. This suggests that since 2016,
the prevalence of a&t domains detected proactively has significantly
increased.

##  V Related Work

Several research works focus on characterization and prevention of online
tracking [30], [32], [36], [35], [38], [39], [40], [42], [41]. We classify the
Web tracking literature in two classes and discuss them below.

Prevalence of online ads and tracking: Krishnamurthy and Wills [35], use
longitudinal measurement snapshots to show that the third-party a&t domains
usage was already increasing significantly between 2005 and 2008. Their work
provides a comprehensive overview of how some organizations (e.g., Google)
have increased their tracking coverage both by increased usage of some of the
third-party domains that they own, and through active acquisition of new
domains (e.g., DoubleClick) that provide ads services, analytics, and
tracking. Metwalley et al., [39] did a passive study on the interaction of
users with a&t domains and ad-blockers. They showed that only one in five
internet users have installed ad-blocker plugins on their web browsers, and
nearly four out of five internet users interact with a&t domains as soon as
they begin to browse. Pujol et al. [40] extended their work using traces by
measuring traffic in ISP. They showed that the ad component is significant in
web traffic, and many users install ad-blockers to block advertisements rather
than protect their privacy.

Our work focuses on the evolution of a&t domains and ad-blocking lists. We
perform a longitudinal study to show the prevalence of a&t domains in the
popular web of different countries. We show that websites embed new a&t
domains over time to evade ad-blocking lists and provide targeted
advertisements to the users. The ad-blocking lists continue to detect new a&t
domains each year. We also provide the update speed of ad-blocking lists to
illustrate if they have been reactive or proactive in handling a&t domains.

Lerner et al. found an increase in the prevalence and complexity of tracking,
as well as an increase in the inter-connectedness of the ecosystem by
analyzing Internet Archive data from 1996 to 2016 [36]. Fruchter et al.
studied geographic variations in tracking [31]. More recently, Libert studied
third-party HTTP requests on the top 1 million sites [37], providing a view of
tracking across the web. In this study, Libert showed that Google can track
users across nearly 80% of sites through its various third-party domains.
Englehardt and Narayanan detected third-party fingerprinting services on Alexa
top 1M websites [29].

Characterization of online ads and tracking: Roesner et al., [41] categorize
mainstreaming tracking and analytics services and studied their prevalence in
the wild. Ikram et al., [32] evaluate the usability of five different tracking
prevention plugins and reveal that contemporary blacklists based tracking
prevention tools are inefficient and have high false positives. Ikram et al.
[33] classify the tools for browser plugins, VPN or proxy-based tools, DNS-
based filtering, and stand-alone browsers. Their work also delineate tracking
prevention techniques for mobile as well as desktop environments. Blacklists
are discussed and machine-learning based techniques are evaluated for tracking
prevention [32, 28].

##  VI Limitations and Concluding Remarks

Despite our best efforts to collect the dataset from the Wayback Machine to
perform our study, there are still some limitations. First, the webpage
requested to Wayback Machine sometimes encounters a 301 redirect, meaning that
the page has permanently moved to a new location, or a 302 redirect meaning
that the move is only temporary. Second, websites are archived with varied
frequency on Wayback Machine. Hence, it is possible to miss an a&t domain
embedded in a website for the duration the website was not archived by Wayback
Machine. When we encounter a web page redirect on a given snapshot, we discard
that snapshot and pick another snapshot from the same three month interval.
Similarly, if we are unable to obtain any snapshot in that three month
interval after four tries, we skip that interval altogether and move to the
next three month interval. For a website, if we are unable to capture at least
a single snapshot for any year between 2009 and 2017, we discard that website
from our analysis.

We have presented an approach to measure the longitudinal performance of ad-
blocking blacklists. We have proposed two metrics, stability and diversity,
that measure the rate of change of ad-blocking blacklists, and introduced two
metrics, reactive and proactive, that measure the update speed of ad-blocking
blacklists. Our analysis revealed that since 2011, more than 50% of the
crawled websites embedded a different number of a&t domains compared to the
preceding year. Ad-blocking blacklists also continue to detect new a&t domains
each year (at most 15% in 2017). Among the ad-blocking blacklists, Mahakala
and hpHosts blocked the most number of a&t domains on the web, but these two
ad-blocking blacklists also contain the highest number of black-listed
domains. On resource-constrained devices, leaner ad-blocking blacklists like
EasyList, EasyPrivacy, and EasyList_China, that provide slightly less coverage
are more suitable. Among the countries we evaluated, the US, Canada, and the
UK have the highest number of distinct a&t domains per website in the popular
web and also the highest proactive scores. This suggests that the ad-blocking
lists update by prioritizing the a&t domains reported in the popular websites
from these countries.

## References

  * [1] “Adblock: content filtering and ad blocking browser extension,” https://getadblock.com.
  * [2] “Adblock Plus: open-source browser extension for content-filtering and ad blocking,” https://adblockplus.org.
  * [3] “CyberCrime Statistics,” https://cybercrime-tracker.net/stats.php.
  * [4] “Ghostery: open-source browser extension and mobile browser application to control JavaScript tags,” https://www.ghostery.com.
  * [5] “Internet Archive: Digital Library of Free & Borrowable Books, Movies, Music & Wayback Machine,” http://archive.org.
  * [6] “Privacy Badger: open-source browser extension for blocking advertisements and tracking cookies,” https://www.eff.org/fr/node/99095.
  * [7] “Canonical repository for the Disconnect services file,” https://disconnect.me/trackerprotection, 2018.
  * [8] “Dasient Smart Web Security Q3 2010 Malware Update,” https://web-beta.archive.org/web/20111017194759/http://www.dasient.com/documents/Dasient_3Q_2010_Qtrly_Malware_F.pdf, 2018\.
  * [9] “AdAway Hosts,” https://adaway.org/hosts.txt, 2019.
  * [10] “AdGuard,” https://adguard.com, 2019.
  * [11] “AdZHosts,” https://adzhosts.eu, 2019.
  * [12] “Cameleon,” http://sysctl.org/cameleon, 2019.
  * [13] “Cybercrime-tracker,” https://cybercrime-tracker.net, 2019.
  * [14] “EasyList,” https://easylist.to/easylist/easylist.txt, 2019.
  * [15] “EasyList China+EasyList,” https://easylist-downloads.adblockplus.org/easylistchina+easylist.txt, 2019\.
  * [16] “EasyPrivacy,” https://easylist.to/easylist/easyprivacy.txt, 2019.
  * [17] “Fanboy’s Social Blocking List,” https://easylist-downloads.adblockplus.org/fanboy-social.txt, 2019.
  * [18] “hpHosts Online – Simple, Searchable & FREE!” http://www.hosts-file.net, 2019.
  * [19] “Mahakala Adblocking Hosts,” https://adblock.mahakala.is/, 2019.
  * [20] “Momento: Time Travel,” http://timetravel.mementoweb.org, 2019.
  * [21] “MVPs Hosts Lists,” http://winhelp2002.mvps.org/hosts.txt, 2019.
  * [22] “Pi-hole: A black hole for Internet advertisements,” https://pi-hole.net, 2019.
  * [23] “Privoxy - Home Page,” https://www.privoxy.org, 2019.
  * [24] “Proxomitron.Info …the Webhiker’s Guide to Proxomitron.” https://proxomitron.info, 2019.
  * [25] “QuidsUp,” https://quidsup.net, 2019.
  * [26] “Sa-blacklist,” http://www.stearns.org/sa-blacklist, 2019.
  * [27] A. Barth, C. Jackson, and J. C. Mitchell, “Securing frame communication in browsers,” _Communications of the ACM_ , 2009.
  * [28] S. Bhagavatula, C. Dunn, C. Kanich, M. Gupta, and B. Ziebart, “Leveraging machine learning to improve unwanted resource filtering,” in _Proceedings of the 2014 Workshop on Artificial Intelligent and Security Workshop_. ACM, 2014, pp. 95–102.
  * [29] S. Englehardt and A. Narayanan, “Online tracking: A 1-million-site measurement and analysis,” in _Proceedings of ACM CCS 2016_ , 2016.
  * [30] M. Falahrastegar, H. Haddadi, S. Uhlig, and R. Mortier, “The rise of panopticons: Examining region-specific third-party web tracking.” in _TMA_ , 2014, pp. 104–114.
  * [31] N. Fruchter, H. Miao, S. Stevenson, and R. Balebako, “Variations in tracking in relation to geographic location,” _CoRR_ , vol. abs/1506.04103, 2015. [Online]. Available: http://arxiv.org/abs/1506.04103
  * [32] M. Ikram, H. J. Asghar, M. A. Kaafar, A. Mahanti, and B. Krishnamurthy, “Towards seamless tracking-free web: Improved detection of trackers via one-class learning,” _Proceedings on Privacy Enhancing Technologies_ , vol. 2017, no. 1, pp. 79–99, 2017.
  * [33] M. Ikram and M. A. Kaafar, “A first look at mobile ad-blocking apps,” in _Network Computing and Applications (NCA), 2017 IEEE 16th International Symposium on_. IEEE, 2017, pp. 1–8.
  * [34] M. Ikram, R. Masood, G. Tyson, M. A. Kaafar, N. Loizon, and R. Ensafi, “The chain of implicit trust: An analysis of the web third-party resources loading,” in _WWW_ , 2019.
  * [35] B. Krishnamurthy and C. Wills, “Privacy diffusion on the web: a longitudinal perspective,” in _Proceedings of the 18th international conference on World wide web_. ACM, 2009, pp. 541–550.
  * [36] A. Lerner, A. K. Simpson, T. Kohno, and F. Roesner, “Internet jones and the raiders of the lost trackers: An archaeological study of web tracking from 1996 to 2016,” in _USENIX Security_ , 2016.
  * [37] T. Libert, “Exposing the hidden web: An analysis of third-party HTTP requests on 1 million websites,” _CoRR_ , vol. abs/1511.00619, 2015. [Online]. Available: http://arxiv.org/abs/1511.00619
  * [38] J. R. Mayer and J. C. Mitchell, “Third-party web tracking: Policy and technology,” in _Security and Privacy (SP), 2012 IEEE Symposium on_. IEEE, 2012, pp. 413–427.
  * [39] H. Metwalley, S. Traverso, M. Mellia, S. Miskovic, and M. Baldi, “The online tracking horde: a view from passive measurements,” in _International Workshop on Traffic Monitoring and Analysis_. Springer, 2015, pp. 111–125.
  * [40] E. Pujol, O. Hohlfeld, and A. Feldmann, “Annoyed users: Ads and ad-block usage in the wild,” in _Proceedings of the 2015 ACM Conference on Internet Measurement Conference_. ACM, 2015, pp. 93–106.
  * [41] F. Roesner, T. Kohno, and D. Wetherall, “Detecting and defending against third-party tracking on the web,” in _Proceedings of the 9th USENIX conference on Networked Systems Design and Implementation_. USENIX Association, 2012, pp. 12–12.
  * [42] F. Roesner, C. Rovillos, T. Kohno, and D. Wetherall, “Sharemenot: Balancing privacy and functionality of third-party social widgets,” _Usenix; login_ , 2012.
  * [43] A. Vastel, P. Snyder, and B. Livshits, “Who filters the filters: Understanding the growth, usefulness and efficiency of crowdsourced ad blocking,” _arXiv preprint arXiv:1810.09160_ , 2018.
  * [44] C. E. Wills and D. C. Uzunoglu, “What ad blockers are (and are not) doing,” in _HotWeb_ , 2016.
  * [45] S. Zhu, X. Hu, Z. Qian, Z. Shafiq, and H. Yin, “Measuring and disrupting anti-adblockers using differential execution analysis,” in _NDSS_ , 2018\.