TITLE: Generic Question Regarding modular arithmetic
QUESTION [2 upvotes]: I was going over a paper regarding linear algebra and its relation to chaos.  One particular section focuses on Arnold's cat map. It briefly mention modular arithmetic.
In this case, it dealt with the generic form of $x \hspace{1mm} mod \hspace{1mm} 1$.
More specifically, $x \hspace{1mm} mod \hspace{1mm} 1$ denotes the unique number in the interval $[0,1)$ that differes from $x$ by an integer.
That said, the paper gave a few examples:
$2.3 \hspace{1mm} mod \hspace{1mm} 1 = 0.3$
$0.9\hspace{1mm}mod \hspace{1mm}1 = 0.9$
$-3.7\hspace{1mm}mod\hspace{1mm}1 = 0.3$
$2.0\hspace{1mm}mod\hspace{1mm}1 = 0$
$(2.3, -7.9)\hspace{1mm}mod\hspace{1mm}1 = (0.3, 0.1)$
Based on those examples, what is the fundamental reason behind $0.9\hspace{1mm}mod \hspace{1mm}1 \hspace{1mm} equaling \hspace{1mm} 0.9$ and $0.3\hspace{1mm}mod \hspace{1mm}1 \hspace{1mm} equaling \hspace{1mm} 0.3$?  In more detail, why couldn't the latter statement be equal to $0.7$?
Obviously, that isn't the answer, but is it a concept similar to why we round 7.5 up to 8 (i.e. the 0.5)?

REPLY [1 votes]: $$x \equiv a \pmod{1}$$
is sometimes used a short hand notation to denote the fractional part of $x$, where $a$ denotes the fractional part of $x$ i.e. $$x = \lfloor x \rfloor + \{x\}$$ where $\lfloor x \rfloor$ is the greatest integer less than or equal to $x$ and $\{x\} \in [0,1)$ is the fractional part of $x$.
Another interpretation of  $x \pmod{1}$ is the distance of $x$ from the largest integer no larger than $x$ i.e. the largest integer that falls to the left of $x$ or $x$ itself on the real number line. I have tried to illustrate this with couple of diagrams below.
The first picture indicates $x\pmod{1}$ when $x > 0$.

The second picture indicates $x\pmod{1}$ when $x <0$.


TITLE: Relationship between different types of correlation coefficients
QUESTION [0 upvotes]: Let,
$r_{1(2.34...p)}$ = Correlation between $x_1$ and $x_{2.34...p}$. The latter being the residuals after regressing $x_2$ on $x_3 , x_4 ....x_p$.
$r_{1.234..p}$ = Multiple correlation coefficient of regressing $x_1$ on $x_2 , x_3, x_4....x_p$
Prove that -
${r_{1.23...p}}^2 = {r_{1p}}^2 + {r_{1(p-1.p)}}^2 + ...... + {r_{1(2.34...p)}}^2$
I tried writing the $correlation^2$ coefficients first in terms of $covariance^2$ by variance*variance. Variance of $x_1$ will cancel out from both the sides. Then I tried substituting all the residuals/fitted values in the covariances with linear combinations of ${x_i}'s$, but to no avail. How to prove this equality?

REPLY [0 votes]: ${r_{1.23...p}}^2 = {r_{1p}}^2 + {r_{1(p-1.p)}}^2 + ...... + {r_{1(2.34...p)}}^2$
$r_{12.34...p}$ = Partial Correlation between 1 and 2 removing the effects of 3,4,...p.
$x_{1.34...p}$ = Residuals of 1 after regressing on 3,4,...p
$s_{11.34....p}$ = Variance of residuals of 2 after regressing on 3,4...p
$r_{12.34...p}^2$ = $\left({Cov(x_{1.34...p},x_{2.34...p})}\over {\sqrt(s_{11.34....p}) \sqrt(s_{22.34...p})}\right)^2$
$r_{12.34...p}^2$ = $\left({Cov(x_{1},x_{2.34...p})}\over {\sqrt(s_{11.34....p}) \sqrt(s_{22.34...p})}\right)^2$
Because, normal equations, residuals of $x_1$ with $x_3,x_4...x_p$ will give zero when multiplied with $x_3,x_4....x_p$
Multiplying dividing with $(\sqrt(s_{11}))^2$
$r_{12.34...p}^2$ = $(r_{1(2.34...p)})^2 \times s_{11} \over s_{11.34....p}  $

Now,
$x_{1.23....p}$ = Values of $x_1$ regressed on $x_2, x_3....x_p$
We look at
$\sum_{i} ((x_{1.23....p})_i)^2 = \sum_{i} ((x_1)_i)\times((x_{1.23....p})_i)$
= $\sum_{i} ((x_{1.34...p})_i)\times ((x_{1.23....p})_i)$
Writing $((x_{1.23....p})_i) = ((x_1)_i) - \sum_{j=2}^{p}b_j \times ((x_j)_i)$
= $\sum_{i} ((x_{1.34...p})_i)\times((x_1)_i) - \sum_{j=2}^{p}((x_{1.34...p})_i)\times b_j \times ((x_j)_i)$
= $\sum_{i} ((x_{1.34...p})_i)\times((x_1)_i) - ((x_{1.34...p})_i)\times b_2 \times ((x_2)_i)$
We know, $b_2 = b_{12.34...p}$ ( Coefficient of $x_2$ when $x_1$ is regressed on $x_2,x_3...x_p$ is same as partial relation coefficient between residuals of $x_1$ and $x_2$ after removing the effects of $x_3,x_4...x_p$
= $\sum_{i} ((x_{1.34...p})_i)\times((x_1)_i) - ((x_{1.34...p})_i)\times b_{12.34...p} \times ((x_2)_i)$
= $\sum_{i} ((x_{1.34...p})_i)\times((x_1)_i) - ((x_{1.34...p})_i)\times b_{12.34...p} \times ((x_{2.34...p})_i)$
= $\sum_{i} (x_{1.34...p})_i) \times (((x_1)_i) - b_{12.34...p} \times (((x_{2.34...p})_i)$
= $\sum_{i} ((x_{1.34...p})_i))^2 - ((x_{1})_i)\times b_{12.34...p} \times ((x_{2.34...p})_i)$
So,
$s_{11.23...p} = s_{11.34...p} - b_{12,34,,,p} \times \sum_{i} ((x_1)_i) \times ((x_{2.34...p})_i)$

Using
$b_{12.34...p}$ = $r_{12.34...p} \sqrt{s_{11.34...p}} \over \sqrt{s_{22.34...p}}$
1 - ${r_{12.34...p}}^2$ = $s_{11.23...p} \over s_{11.34..p}$
and
1 - ${r_{1.23...p}}^2$ = $s_{11.23...p} \over s_{11}$
in the two equations derived above, cancelling and manipulating, we will get the desired result.


TITLE: Finding distribution of $\frac{1}{\sqrt X}$ when $X$ is uniform on $[1,9]$
QUESTION [2 upvotes]: I haven't been able to figure this question out and have been googling for about a good hour so thought I should post this here.

Let $X$ be a $\mathrm{uniform}(1,9)$ random variable, i.e. it has the probability density function
  $$
  f_X(x) =
\begin{cases}
1/8,  & \text{if 1 $\le$  $x$ $\le$ 9,} \\
0, & \text{otherwise.}
\end{cases}
$$
  Find the Cumulative Distribution function and Probability density function of
  $$Y = \frac {1}{\sqrt{X}}.$$

My wrong attempt went as follows:
First I determined the range of $ Y $ to be $y \in [\frac13,1] $ by filling in the boundrairies 1 and 9 into $ Y $
Calculate the CDF by taking the close integral over the interval $ [1,x]$
$$ CDF: \int_1^x \frac18 du =  \frac18(x-1)$$
Then I find CDF of $ Y $ $$  F_y = P(Y \le y )  =  P( \frac{1}{\sqrt{X}} \le y ) = P(X \le \frac{1}{y^2} ) = \frac18(\frac{1}{y^2} -1) $$
But this solutions seems to be wrong since it doesn't add up to 1 when integrated over its range $ [\frac13,1]$
If someone could give me an explanaiting where I made an error that be great!

REPLY [2 votes]: Mistake is in the line $\mathbb P(\frac{1}{\sqrt{X}} \le y) = \mathbb P(X \le \frac{1}{y^2}),$ because due to $\sqrt{X},y > 0$, we have:
$ \frac{1}{\sqrt{X}} \le y  \iff 1 \le \sqrt{X}y \iff \frac{1}{y} \le \sqrt{X} \iff \frac{1}{y^2} \le X,$ so it should be $P(\frac{1}{\sqrt{X}} \le y) = \mathbb P(X \ge \frac{1}{y^2})$
As you calculated $F_X(t) = \frac{x-1}{8}\chi_{[1,9]}(t) $ and correcting your error, we get $($for $y \in [\frac{1}{3},1])$:
$F_Y(y) = \mathbb P(\frac{1}{\sqrt{X}} \le y) = \mathbb P(X \ge \frac{1}{y^2}) = 1 - \mathbb P(X < \frac{1}{y^2}) = 1 - \mathbb P(X \le \frac{1}{y^2}) = 1 - F_X(\frac{1}{y^2}) = 1 - \frac{1}{8}(\frac{1}{y^2}-1) = \frac{1}{8}(9-\frac{1}{y^2}) $
So we get the CDF of $Y$: $F_Y(t) = \frac{1}{8}(9-\frac{1}{t^2})\chi_{[\frac{1}{3},1]}(t)$
Now PDF, as we know since CDF is continuosly differentiable (piecewise at least) it holds, that $f_Y(y) = F_Y'(y)$ for every point in which $F_Y'$ exists (we can put whatever we want in the rest, since $PDF$ is up to the set of measure zero).
Differentiating gives us PDF : $f_Y(y) = \frac{1}{4y^3} \chi_{(\frac{1}{3},1)}(y)$ (Note that you could have included points $\frac{1}{3}$ and/or $1$ in characteristic function.


TITLE: Asymptotic expansion of integral (Bessel function really)
QUESTION [2 upvotes]: The integral
$$I = \int_{-\infty}^\infty \frac{e^{-\varepsilon x^2}} { \sqrt{1+x^2} } dx$$
is convergent for $\varepsilon > 0$ and can even be given in terms of the Bessel function $K_0$. As $\varepsilon \to 0$ it is divergent and $I \sim -\log \varepsilon$. What would be the simplest way to derive the above leading term in an asymptotic $\varepsilon \to 0$ expansion directly in terms of the above integral?
Clearly, if one uses the exact result in terms of $K_0$ and then for instance uses the second order differential equation satisfied by $K_0$ it is quite simple to derive the $\log\varepsilon$ form of the divergence. But I'm looking for a way to derive this directly from the above integral. The reason is that I have a much more complicated integral to analyze which can not be given in a closed form but the above simple integral captures its difficulty so would like to understand this one first.
Another related question: the integral
$$J = \int_0^{2\pi} e^{-\sin(x)^2/\varepsilon^2}dx $$
can also be evaluated exactly in terms of the Bessel function $I_0$ which result will imply $J \sim \varepsilon$ as $\epsilon\to 0$. But again, directly from the integral what is the simplest way to see this leading term in the $\varepsilon\to 0$ expansion?

REPLY [2 votes]: Differentiate $I(\varepsilon):=2\int_0^\infty e^{-\varepsilon x^2}(1+x^2)^{-1/2}dx$    with respect to $\varepsilon$ under the sign of integral; change variable putting $u:=\varepsilon x^2$. We get
$$I'(\varepsilon) = -\frac{1}{\varepsilon} \int_0^\infty e^{-u}\sqrt{\frac{u}{u+\varepsilon}}du= -\frac{1}{\varepsilon}\big(1+o(1)\big)\, ,$$
by the dominated convergence theorem, and integrating
$$I(\varepsilon)=-\log(\varepsilon)\big(1+o(1)\big),\,\,  \mathrm{as }\,  \varepsilon\to0\, .$$
For the latter integral, change variable with $x:=\varepsilon u $, so
$$J(\varepsilon):=2\int_{-\frac{\pi}{2}}^{+\frac{\pi}{2}}e^{-{\sin(x)^2}/{\varepsilon^2}}dx=2\varepsilon\int_\mathbb{R} e^{-\sin(\varepsilon u)^2/\varepsilon^2}\chi_{ [-\frac{\pi}{2\varepsilon},\, +\frac{\pi}{2\varepsilon} ] }(u) du\, .$$
The integrand converges pointwise to $e^{-u^2}$ as $\varepsilon\to0$, and it is dominated by $e^{-u^2/4}$ for any $u\in\mathbb{R}$ (just because  $\sin(x)\ge x/2\ge 0$ for any $0\le x \le \pi/2$ ). By the dominated convergence theorem,
$$J(\varepsilon)=2\varepsilon\int_\mathbb{R} e^{-u^2}du\big(1+o(1)\big)=2\sqrt{\pi}\varepsilon\big(1+o(1)\big),\, \,  \mathrm{as }\,  \varepsilon\to0\, .$$


TITLE: Superscript Function Notation(variable as base) $x^{\phi ^{-1}}$
QUESTION [0 upvotes]: I apologize if this is a silly question but I cannot find the necessary information elsewhere on the internet: If $\phi : u \mapsto ux$ is a function from $U$ to $U$ , does the notation $x^{\phi ^{-1}}$ mean the inverse function of $\phi $ where the domain variable is $x$??
For full context of the notation Ill put a link to Kurzweil-Theory book on groups where it is on page 5 (page 18 of pdf):
https://www.google.com/url?sa=t&source=web&rct=j&url=http://www.math.ku.dk/~olsson/manus/GruFus/Kurzweil-Stellmacher_Theory%2520of%2520finite%2520groups.pdf&ved=0ahUKEwiYmsTnrODPAhVB5SYKHcXHAMcQFggbMAA&usg=AFQjCNH-EJPLVWDku5ZBY2zwQqlEGCsz3g&sig2=MuavjBUWoHaq38F-F17FtA

REPLY [1 votes]: It is not uncommon in group theory to write $x^g$ for $g^{-1}xg$: i.e., what you get if you apply conjugation by $g$ to $x$, where $x$ and $g$ are members of some group $G$ (this can save quite a lot of ink). The conjugation function $x \mapsto x^g$ is called an inner automorphism on $G$. By extension of this notational idea, it isn't uncommon to write $x^{\alpha}$ for the result $\alpha(x)$ of applying an arbitrary automorphism $\alpha$ on $G$ to $x$.
The authors of this book have clearly taken this further and write $x^{f}$ for $f(x)$ for functions $f$ in general, as in in your quotation (where $f$ is $\phi^{-1}$). The notation $\phi : x \mapsto ux$ is a fairly standard way of saying "the function $\phi$ defined by $\phi(x) = ux$".


text:
\begin{document}

\title{ Adjoint for advection schemes on the sphere in ICON model}

\author{Ramaz Botchorishvili}
\address{Faculty of Exact and Natural Sciences, Ivane Javakhishvili Tbilisi State University, Tbilisi, Georgia}
\email[Corresponding author]{ramaz.botchorishvili@tsu.ge (Corresponding author)}

\author{Hendrik Elbern}
\address{Rhenish Institute for Environmental Research, University of Cologne, Cologne, Germany;
         Institute of Energy and Climate research 8: Troposphere, Research Center J\"ulich, J\"ulich, Germany }
\email{h.elbern@fz-juelich.de}

\author{Tamari Janelidze}
\address{Faculty of Exact and Natural Sciences, Ivane  Javakhishvili Tbilisi State University, Tbilisi, Georgia;
         Institute of Energy and Climate research 8: Troposphere, Research Center J\"ulich, J\"ulich, Germany }
\email{t.janelidze@fz-juelich.de}

\subjclass[2010]{65M32, 65M08, 65Z99}
\keywords{adjoint solver, linear advection equation, ICON model, variational data assimilation,  inverse modeling}

\date{\today}

\begin{abstract}
		Among the most advanced and sophisticated methods for state analysis of an atmospheric
		system is the four dimensional variational data assimilation. The numerically challenging
		task of this approach is the development and application of the adjoint model components.
		For tracer transport in fluid dynamics accuracy of numerical advection schemes is vital.
		It is even more important for applications in space-time variational data assimilation with
		adjoint model version. We propose novel straightforward and efficient approach - artificial
		source term method - for adjoint advection solver development. It has several benefits
		compared to traditional adjoint model building technique. One of the attractive features of
		the new approach is that it reuses existing advection solver code, thus resulting into
		significant reduction of time needed for adjoint solver development. The stability, accuracy
		and convergence of the adjoint schemes are investigated. The method is implemented and evaluated
		for the linear advection equation on the sphere in the Icosahedral Nonhydrostatic Model (ICON).
		Adjoint solvers developed by the conventional and the methods are compared against each other on a
		collection of standard advection test cases and variational data assimilation test cases developed
		here. The advantages of the artificial source term method is especially obvious in case of monotonic
		advection equation solvers, granting, e.g. absence of oscillations and nonphysical negative
		concentrations.
\end{abstract}

\maketitle

\section{Introduction}
	Chemistry transport models (CTM) are widely used for air pollution modelling.  These models
	are complex and they are usually solved numerically using operator splitting methods when
	equations describing chemical reactions, horizontal advection and vertical advection  are treated separately
	at different splitting steps, see e.g. ~\cite{ElbernSchwingerBotchorishvili:Chemical}.
	Incorporating observations from both in situ and remote sensing devices in CTMs, e.g. satellites etc., is
	essential for the quality of air pollution forecast. The method used for this purpose is called data
	assimilation that first appeared in meteorology ~\cite{LeDimetTalagrand1986}. Four dimensional variational
	data assimilation integrates into the model observational data distributed in space and finite time interval.
	For achieving this goal a constrained minimization problem is formulated where governing equations serve
	as constraints and a quadratic cost function measures misfits between observations and model forecast.
	Solving this minimization problem results in optimal initial data that are consistent with observations and
	with model dynamics. These data are used for improved forecasting on longer time interval,
	see  ~\cite{TalagrandCourtier1987},~\cite{ElbernSchmidtEbel1997},~\cite{Wang2001},~\cite{SanduChai2011}
	for detailed exposition on the subject, see section ~\ref{sec:variational} for the formulation of variational
	problem for the linear advection equation.

	A central module in CTMs as part of a 4D-var system is  the advcetion algorithm and its ajoint.  \cite{vukicevic2001properties}
	studied the influence of linear and non--linear numerical advection algorithm properties on variational data assimilation
	results in a 2D idealized scalar advection framework. Results suggested that  exact the same scalar advection algorithm in
	forward and adjoint computations obtains, at lower cost, an optimal solution accuracy that is consistent with the forward
	model accuracy.

	\cite{thuburn2001adjoints} investigated the impact of switches in non-oscillatory advection schemes for their adjoints.
	They  showed that there is no possibility of smoothing the switches in nonoscillatory advection schemes to remove
        the discontinuities while retaining an obvious and desirable scaling property. On the basis of established
        equivalence between Eulerian backtracking or retro-transport and adjoint transport with respect to an
        air-mass-weighted scalar product, \cite{hourdin2006eulerian} studied the question which arises as to whether it is
        preferable to use the exact numerical adjoint, or the retro-transport model for a model that is not time-symmetric.
        They concluded that the presence of slope limiters in the Van Leer advection scheme can produce in sonic circumstances
        unrealistic, even negative adjoint sensitivities. The retro-transport equation, on the other hand, generally produces
        robust and realistic results.

        In the contex of the GEOS-Chem \cite{henze2007development} tested the accuracy of the adjoint model by comparing adjoint to
        finite difference sensitivities, which are shown to agree within acceptable tolerances. They explore the robustness
        of these results, noting how discontinuities in the advection routine hinder, but do not entirely preclude, the use of such
        comparisons for validation of the adjoint model.

        A  comprehensive study on the consistency of the discrete adjoints of upwind numerical schemes was provided by
        \cite{liu2008properties}. Both linear  and nonlinear  discretizations of the one-dimensional advection equation are considered,
        representing finite differences or finite volumes and slope or flux-limited techniques, respectively. \cite{gou2011continuous}
        studied the effect of using discrete and continuous adjoints of the advection equation in chemical transport modeling numerically.
        They concluded that discrete advection adjoints are more accurate in point-to-point comparisons against finite differences,
        whereas the continuous adjoints of advection perform better as gradients for optimization in 4D-Var data assimilation.

        Considering  a variational implementation of the RETRO-TOM model, \cite{haines2014adjoint} demonstrated a time symmetric forward
        advection scheme with second-order moments to be efficiently exploited in the backward adjoint calculations, at least for
        problems in which flux limiters in the advection scheme are not required. For this case the authors found the flexibility and
        stability of a 'finite difference of adjoint' formulation with the accuracy of an 'adjoint of finite difference' formulation.

        \cite{holdaway2015assessing} examined the tangent linear and adjoint versions of NASA's Goddard Earth Observing System version 5
        (GEOS-5). Tests exhibited that piecewise parabolic methods with flux limiters development of unrealistically large perturbations
        within the tangent linear and adjoint models, and that using a linear third-order scheme for the linearised model produces better
        behaviour.

	Iterative methods are used for minimizing the quadratic cost function. Quasi-Newton limited memory methods are
	among most popular minimization algorithms  ~\cite{LiuNocedal1989}. A numerically challenging task in variational
	data assimilation is the development of the adjoint model for the gradient computation that is used by quasi
	Newton algorithm.  For the linear advection equation, as demonstrated in ~\cite{thuburn2001nonoscillatory},
	the results depend on the implementation of adjoint advection schemes especially in case of nonlinear schemes.
	Several different approaches can be used for adjoint solver development, see e.g. ~\cite{LeDimetTalagrand1986},
	\cite{nodet2016variational} for  derivation of adjoint equations in continuous setting for which numerical solver
	has to be developed, see e.g. ~\cite{thuburn2001adjoints}(სანახავია) for the development of adjoint in discrete
	setting - i.e. adjoint of a given numerical advection scheme,  and see ~\cite{hascoet2013tapenade} for the
	development of adjoint code based on powerful automated differentiation tool Tapenade.

	Here we propose simple and powerful approach for the adjoint development for linear advection equation in
	conservative form. Our approach is based on introducing artificial source term that ensures consistency of
	the adjoint numerical scheme with the adjoint of the linear advection equation under consideration.
	Implementation of our approach is easy and time efficient: it essentially reuses given linear advection solvers
	with minor modifications, maintains accuracy and stability of the given parent solver and it uses the same
	parallelization routines. The proposed method is general, though here we consider advection schemes on triangular
	mesh of the sphere in the ICON model. See subsection \ref{sec:ICONmesh} for short description
	and ~\cite{wan2013icon} for details, as part of broader effort on development of a next generation 4D-Var data
	assimilation system with CTM.

	The paper is organized as follows: in section \ref{sec:variational} the linear advection equation on the sphere
	is given, the variational problem on data assimilation is formulated, and the adjoint and cost functions are given;
	in section \ref{sec:NumSchemes} a new method for building adjoint scheme is proposed that is applied for building
	the adjoint solver using linear advection schemes in the ICON model, the properties  of the developed adjoint schemes
	are investigated theoretically; in section \ref{sec:NumTests} the new method is studied numerically on test problems
	in the context of data assimilation for tracer transport models.

\section{Variational problem formulation for linear advection}
\label{sec:variational}
\subsection{Linear advection equation on the sphere}
The linear advection equation in conservative form writes:
\begin{equation}\label{eq:LinearAdvection}
\frac{\partial(\rho q)}{\partial t}+\nabla\cdot(\rho q\vec{v})=0,
\end{equation}
where  $\rho$ and $q$ are fluid density and mixing ratio, respectively. We are interested in their advection on the surface of the sphere $\Omega$.
Therefore $\vec{v}=(v_\lambda,v_\theta)^T$ is the 2D horizontal wind vector and $\nabla\cdot$ is a spherical horizontal divergence operator given by
\begin{equation}\label{eq:NablaHorizontal}
	\nabla\cdot v=\frac{1}{\sin\theta}\big[\frac{\partial v_\lambda}{\partial \lambda}+\frac{\partial (v_\theta\sin\theta)}{\partial\theta}\big],
\end{equation}
with  $0\leq\lambda\leq 2\pi$ longitude and $0\leq\theta\leq \pi$ latitude. For simplicity, we take the radius $R$ of the sphere to be one.

Throughout this paper in (\ref{eq:LinearAdvection})  $\rho=\rho(t,\lambda, \theta)>0$ and $v=v(t,\lambda, \theta)$ are sufficiently smooth given functions,
$q=q(t,\lambda, \theta)$ is the unknown function that should be determined from the equation (\ref{eq:LinearAdvection}) subject to the initial
condition that writes:
\begin{equation}\label{eq:InitCondition}
q(0,\lambda,\theta)=q_0(\lambda,\theta),
\end{equation}
where $q_0(\lambda,\theta)$ is a sufficiently smooth given function.

Notice that initial condition (\ref{eq:InitCondition}) is given on the surface of the sphere. Because of this initial condition is sufficient and no
boundary conditions are needed, the problem (\ref{eq:LinearAdvection})-(\ref{eq:InitCondition}) is well posed for any finite time interval $[0,T]$.

\subsection{Variational problem, adjoint equation, gradient of cost function}

Consider the following cost function:
\begin{equation}\label{eq:CostFunc}
J(q_0)=f^b(q_0,q^b)+\int_{0}^{T}\int_{\Omega}f^o(q,q^o,t,\lambda,\theta)d\lambda d\theta dt,
\end{equation}
where   $q_0$ is initial value function from the initial condition (\ref{eq:InitCondition}); $q^b=q^b(\lambda, \theta)$ is given, it is the so called
background function that is also used as first guess for iteration method minimizing $J(q_0)$ with respect to $q_0$; in $q^o$ the superscript "o" stands
for observations and is a given function; the functions $f^b$ and $f^o$ will be detailed in next subsection. Function $f^b$ measures the misfit between initial and
background functions, while $f^o$ measures misfit between mixing ratio and observations. Here we assume $f^b$ and $f^o$ are smooth enough, we also assume that
gradients of $f^b$ and $f^o$ with respect to $q_0$, denoted here via $\nabla_{q_0} f^b$ and $\nabla_{q_0} f^o$ respectively, can be computed analytically.

Now we can state the variational data assimilation problem in the following way:

\begin{problem}[Variational data assimilation]
Find a initial function $q_0$ minimizing cost function (\ref{eq:CostFunc}) subject to constraints (\ref{eq:LinearAdvection}),(\ref{eq:InitCondition}).
\end{problem}

The gradient of cost function $J(q_0)$ writes:
\begin{equation}\label{eq:GradCostFunc}
\nabla_{q_0} J(q_0)=\nabla_{q_0}f^b(q_0,q^b) - \rho(0,\lambda, \theta)q^*(0,\lambda, \theta),
\end{equation}
where $q^*(0,\lambda, \theta)$ is solution at $t=0$ of the following adjoint equation
\begin{equation}\label{eq:AdjointLinearAdvection}
\rho (\frac{\partial q^*}{\partial t}+v\cdot \nabla q^*)=\nabla_{q}f^o(q,q^o,t,\lambda,\theta),
\end{equation}
which is subject to the initial condition at $t=T$, in particular,
\begin{equation}\label{eq:AdjointInitCondition}
q^*(T,\lambda,\theta)=0.
\end{equation}

Notice that in the adjoint equation (\ref{eq:AdjointLinearAdvection}) on the right hand side the function $q$ is the solution of the conservative linear advection
equation (\ref{eq:LinearAdvection}) with initial condition (\ref{eq:InitCondition}). Notice also that the problem  (\ref{eq:LinearAdvection}),(\ref{eq:InitCondition})
is integrated forward in time and the problem (\ref{eq:AdjointLinearAdvection}),(\ref{eq:AdjointInitCondition}) is integrated backward in time.
Because of this reason the problem  (\ref{eq:LinearAdvection}),(\ref{eq:InitCondition}) is often referred as forward problem and the problem (\ref{eq:AdjointLinearAdvection}),
(\ref{eq:AdjointInitCondition}) is referred as adjoint problem. The same applies to the numerical solvers, i.e. forward solver and adjoint solver. Using these terms we
can easily formulate procedure for the gradient computation. In particular, we have:

\begin{procedure}[Computing the gradient of the cost function] \hspace{2cm}
\begin{enumerate}
\item Solve the forward problem (\ref{eq:LinearAdvection}),(\ref{eq:InitCondition}) and store $q(t,\lambda,\theta)$.  \\
\item Put $q(t,\lambda,\theta)$ on the right hand side of the equation (\ref{eq:AdjointLinearAdvection}) and solve the adjoint problem (\ref{eq:AdjointLinearAdvection}),
      (\ref{eq:AdjointInitCondition}), store  $q^*(0,\lambda,\theta)$.\\
\item Insert $q^*(0,\lambda,\theta)$ in (\ref{eq:GradCostFunc}) and compute the gradient of the cost function $J(q_0)$ defined by (\ref{eq:CostFunc}). \\
\end{enumerate}
\end{procedure}

\subsection{Cost function}

Cost functions of variational data assimilation problem are mainly given in discrete or semi-discrete form in scientific literature,
see e.g.~\cite{ElbernSchwingerBotchorishvili:Chemical, Wang2001, SanduChai2011}  for fully discrete cost functions and
see  ~\cite{ElbernSchmidtEbel1997, elbern1999four, Wang2001} for semi-discrete cost functions.  Under semi-discrete version
we mean cost function where continuous time integration is used. Here we establish a relationship between discrete,
semi-discrete and continuous cost functions that is needed later for studying convergence of developed numerical schemes.
For this purpose we have to make a discretization in space and time of dependent and independent variables.
Since finite volume schemes are of interest of the present paper, we denote by $\Omega_j$ finite volume cell $j$,
$ \Omega_j \subset \Omega$, $\cup_{j=1}^{N_c} \Omega_j= \Omega$, where $|\Omega_j|$ is area of the cell $\Omega_j$,
and $N_c$ is number of cells. The average value of function q on $\Omega_j$ at time $t$ is denoted as $q_j(t)$.
The same value at time $t=t_n$ is denoted as $q_j^n$, $t_n=n\Delta t, n=0,1,..,N_T$, where $\Delta t= T/N_T$ is
time discretization step, and $N_T$ is number of nodal points in time. In the data assimilation context all discretized function
values at time $t$ are interpreted as vectors, e.g. $\vec{q}(t)=(q_1(t),q_2(t),..,q_N(t))^T$, $\vec{q}^n=(q_1^n,q_2^n,..,q_{N_c}^n)^T$,
since such notations are convenient for formulating cost function. Adopting our notations the semidiscrete cost function
from ~\cite{elbern1999four} writes:
\begin{equation}\label{eq:CostFuncSemidiscrete}
J(\vec{q}_0))=\frac{1}{2}(\vec{q~}^b-\vec{q}_0)^TB^{-1}(\vec{q~}^b-\vec{q}_0)+ \frac{1}{2}\int_{t_0}^{t_{N_T}}(\vec{q~}^o(t)-\textit{H}[\vec{q}(t)])^TR^{-1}(\vec{q~}^o(t)-\textit{H}[\vec{q}(t)])dt,
\end{equation}
where $B$ and $R$ matrices are the background and observation error covariance matrices, respectively,
$\vec{q~}^o(t)=(q^o_{i_1}(t),q^o_{i_2}(t),..,q^o_{i_{N_o}}(t))^T$, $1 \leq i_k \leq N_c, k=1,2,..,N_o$,  $N_o$
is number of observations, $N_o < N_c$. Operator $\textit{H}$ projects the model state to observation space that practically
means the following: $\textit{H}[\vec{q}(t)]=(q_{i_1}(t),q_{i_2}(t),..,q_{i_{N_o}}(t))^T$.
We denote by $\mathcal{O}$ the set of indices corresponding to observations, $\mathcal{O}=\{ i_1,i_2,..,i_{N_o} \}$.

The fully discrete version of the cost function from ~\cite{ElbernSchwingerBotchorishvili:Chemical} writes:
\begin{equation}\label{eq:CostFuncDiscrete}
J(\vec{q}_0))=\frac{1}{2}[\vec{q~}^b-\vec{q}_0]^TB^{-1}[\vec{q~}^b-\vec{q}_0]+
\frac{1}{2} \sum_{n=0}^{N_T}[\vec{q~}^{on}-\textit{H}(\vec{q~}^n)]^TR^{-1}[\vec{q~}^{on}-\textit{H}(\vec{q~}^n)].
\end{equation}

If we assume that covariance matrices $B$ and $R$ are the same in (\ref{eq:CostFuncSemidiscrete}) and (\ref{eq:CostFuncDiscrete}) then
for the compatibility of semi-discrete and discrete cost functions the following modification is necessary:
\begin{equation}\label{eq:CostFuncDiscreteM}
	J(\vec{q}_0))=\frac{1}{2}[\vec{q~}^b-\vec{q}_0]^TB^{-1}[\vec{q~}^b-\vec{q}_0]+
	\frac{T}{2N_T} \sum_{n=0}^{N_T}[\vec{q~}^{on}-\textit{H}(\vec{q~}^n)]^TR^{-1}[\vec{q~}^{on}-\textit{H}(\vec{q~}^n)],
\end{equation}

Notice that (\ref{eq:CostFuncSemidiscrete}) and (\ref{eq:CostFuncDiscreteM}) are compatible in the sense that
when $N_T \longrightarrow \infty$ limit of discrete cost function (\ref{eq:CostFuncDiscreteM}) coincides with
semidiscrete cost function (\ref{eq:CostFuncSemidiscrete}).

Studying covariance matrices is out of scope of this paper. Therefore, without loss of generality and for the convenience
of further exposition, we replace inverse covariance matrices by some positive definite kernels
denoted by $K_b(\xi,\xi')$ and $K_o(\xi,\xi')$, $\xi, \xi' \in \Omega$, $\xi=(\lambda, \theta)$, $\xi'=(\lambda', \theta')$, and we
define the function $f^b$ from (\ref{eq:CostFunc}) in the following way:
\begin{equation}\label{eq:CostFunc:f^b}
	f^b(q_0,q^b)=\frac{1}{2|\Omega|^2} \int_{\Omega} \int_{\Omega} K_b(\xi,\xi')[q_0(\xi)-q^b(\xi)] [q_0(\xi')-q^b(\xi')] d\xi d\xi'.
\end{equation}

After standard discretization of (\ref{eq:CostFunc:f^b}) we have:
\begin{equation}\label{eq:CostFunc:J^b}
J^b =  \sum_{i=1}^{N_c} \sum_{j=1}^{N_c} \frac{|\Omega_i||\Omega_j|}{2|\Omega|^2}   K_b(\xi_i,\xi'_j)[q_{0j}-q^b_i] [q_{0j}-q^b_j]= \frac{1}{2}[\vec{q~}^b-\vec{q}_0]^TK^b[\vec{q~}^b-\vec{q}_0],
\end{equation}
where $K^b$ is a symmetric positive definite matrix with elements
\begin{equation*}
K^b_{ij}=\frac{|\Omega_i||\Omega_j|}{2|\Omega|^2}   K_b(\xi_i,\xi'_j), ~i,j=1,2,..,N_c.
\end{equation*}

Thus consistency is ensured between discrete and continuous versions of background terms (\ref{eq:CostFunc:f^b}) and
(\ref{eq:CostFunc:J^b}), i.e. $\lim_{N_c \rightarrow \infty} J^b = f^b$.
The formulation will allow us to study the impact of mesh refinement on data assimilation in section ~\ref{sec:NumTests}.

By analogy with (\ref{eq:CostFunc:J^b}) we define $J^o$ and its continuous version:
\begin{multline}\label{eq:CostFunc:J^o}
J^o =  \frac{T}{2N_T} \sum_{n=0}^{N_T} \sum_{i \in \mathcal{O}} \sum_{j \in \mathcal{O}} \frac{|\Omega_i||\Omega_j|}{2|\Omega^o|^2}   K_o(\xi_i,\xi'_j)[q_i^{on}-q^n_i] [q_j^{on}-q^n_j]=
\frac{T}{2N_T} \sum_{n=0}^{N_T}[\vec{q~}^{on}-\textit{H}(\vec{q~}^n)]^TK^o[\vec{q~}^{on}-\textit{H}(\vec{q~}^n)],
\end{multline}

\begin{equation}\label{eq:f^o}
f^o(q,q^o,t, \xi)=
\begin{cases}
\frac{[q(t,\xi)-q^o(t,\xi)]}{2|\Omega^o|^2}  \int_{\Omega^o} K_o(\xi,\xi') [q(t,\xi')-q^o(t,\xi')] d\xi', & ~\xi \in \Omega^o,\\
0,  & \text{otherwise},
\end{cases}
\end{equation}
where
\begin{equation*}
K^o_{ij} = \frac{|\Omega_i||\Omega_j|}{2|\Omega^o|^2}   K_o(\xi_i,\xi'_j), ~{i,j \in \mathcal{O}},~\Omega^o=\cup_{i \in \mathcal{O}} \Omega_i.
\end{equation*}
Consequently we have
\begin{equation}\label{eq:Nabla_qf^o}
\nabla_q f^o(q,q^o,t, \xi)=
\begin{cases}
\frac{1}{|\Omega^o|^2}  \int_{\Omega^o} K_o(\xi,\xi') [q(t,\xi')-q^o(t,\xi')] d\xi', & ~\xi \in \Omega^o,\\
0,  & \text{otherwise}.
\end{cases}
\end{equation}

Notice that in (\ref{eq:f^o}) we assume that observations are given on a subdomain $\Omega^o$, $|\Omega^o| >0.$

\section{Numerical schemes}
\label{sec:NumSchemes}
\subsection{ICON grids}
\label{sec:ICONmesh}
The name "ICON"  stands for the joint project of the Max Plank Institute for Meteorology (MPI-M) and the German Weather Service (DWD) on development
of \textbf{ICO}sahedral \textbf{N}onhydrostatic models ~\cite{wan2013icon}. In the earlier version of ICON model triangular and hexagonal discretizations
were tested and later only the triangular version of the icosahedral mesh is used~\cite{wan2013icon}. In this paper we use the triangular version of icosahedral
grids from "Published list of DWD (EDZW) ICON grids"  by MPI-M and DWD given at http://icon-downloads.zmaw.de/. Meshes are referred as "$Rn_rBn_b$ grids" and
from this name one can recover the algorithm how the grid was constructed. In partciular, starting point in any "$Rn_rBn_b$ grid" is an icosahedron with vertices
on the sphere that is entirely projected onto the sphere, then edges along a great arc are devided in $n_r$ equal parts and then each edge of the obtained triangles
are recursively devided $n_b$ times. See ~\cite{Sadourny1968integration, wan2013icon} for detailed description of the algorithm, for different optimization
approaches and for grid characteristics. Some important characteristics of these meshes are given in the  Table 1 in ~\cite{wan2013icon} that shows non uniformity
of $Rn_rBn_b$ grids in terms of triangle area ratios that can reach values as high as $1.53$ for $R2B7$ grid.  The same table also shows that the ratio is increasing
together with mesh refinement. Here we also give some additional characteristics of $Rn_rBn_b$ grids in the  ~\cref{tab:R2BnGrids}.
When grids are refined  the edge ratio is also  increasing and its maximum can reach $1.34$ for $R2B7$ grid. Though inside each triangle the edge ratio is
almost constant and it increases very slowly, for example in particular, from $1.1734$ for $R2B2$ grid to $1.1761$ for $R2B7$ grid. Close to 1 edge ratio also
means that triangular mesh changes smoothly on the surface of the sphere. $R2B7$ is the finest grid we consider in this paper since.
Notice that our cost function accounts for the variability of grid.

\begin{table}[h]
	\centering
	\caption{ICON grids, $R2Bn_b$ }
	\begin{adjustbox}{width=1\textwidth}
		\begin{tabular}{r r r r r r r r}
			\hline
			Grid & Number of  & Number of & Min triangle             & Max:min       & Max:min        & Min edge \\
			     & triangular & triangle  & cell area, $\text{km}^2$ & edge length   & edge length    & length, km \\
			     & cells      & edges     &                          & ratio, global & ratio, triangle & \\
			\hline

				R2B0  & 80      & 120     & 6010381.55 & 1.14 & 1.1350 & 3526.95 \\
				R2B1  & 320     & 480     & 1440873.32 & 1.17 & 1.1660 & 1737.06 \\
				R2B2  & 1280    & 1920    & 333434.84  & 1.21 & 1.1734 & 836.96  \\
				R2B3  & 5120    & 7680    & 78835.01   & 1.24 & 1.1750 & 407.12  \\
				R2B4  & 20480   & 30720   & 18777.28   & 1.27 & 1.1754 & 198.71  \\
				R2B5  & 81920   & 122880  & 4507.50    & 1.30 & 1.1755 & 97.36   \\
				R2B6  & 327680  & 491520  & 1089.56    & 1.32 & 1.1756 & 47.87   \\
				R2B7  & 1310720 & 1966080 & 265.08     & 1.34 & 1.1761 & 23.61   \\
				\hline

		\end{tabular}
	\end{adjustbox}
        \label{tab:R2BnGrids}
\end{table}

\subsection{Linear advection schemes in ICON model}
\label{sec:ICONschemes}
Tracer transport schemes in the ICON triangular version are based on ~\cite{miura2007upwind} for finite volume discretization
and on ~\cite{ollivier2002high} for high order reconstruction, both methods adapted to triangular grids on a sphere. The scheme is
referred as ICON-FFSL (ICON-Flux Form Semi-Lagrangean) in ~\cite{lauritzen2014geoscientific} for its second order version.
In this paper we will also use the same abbreviation for ICON tracer transport schemes. These schemes can be used with or without limiters.
Here we build adjoint schemes for both cases. Adjoint without limiter is given in the section \ref{sec:Adjoint:IntegrByParts} and the adjoint
with limiter is introduced in the section \ref{sec:Adjoint:ArtSource}.  Limiters in ICON-FFSL
schemes are based on  Zalesak's \textbf{F}lux \textbf{C}orrected \textbf{T}ransport (FCT) by ~\cite{zalesak1979fully} the limiter and its
positive definite modifications by  ~\cite{schar1996synchronous} and ~\cite{harris2011flux}.  Details including numerical results on
ICON tracer transport schemes  are given in  http://www.cgd.ucar.edu/cms/pel/transport-workshop/2011/16-Reinert.pdf,
~\cite{wan2013icon} and ~\cite{lauritzen2014geoscientific}. For the convenience of further exposition on constructing adjoint schemes the above
finite volume advection schemes ICON-FFSL can be written in the following flux form:
\begin{equation} \label{fv:FFSL4LA}
	\frac{\rho_j^{n+1} q_j^{n+1}-\rho_j^n q_j^n}{\Delta t}+\frac{1}{| \Omega_j |} \sum_{i\in I_j} F_{ji}(\{\bar{\rho}_k^n, \bar{v}_k^n,q_k^n\}_{k\in K_{ji}},\vec{n}_{ji}, l_{ji})=0,
\end{equation}
where $I_j$ is a set of reference numbers to cell interfaces surrounding cell $\Omega_j$, i.e. $I_j$ consists of three elements in
case of triangular cell.  $F_{ji}$ is a numerical flux function on the cell interface between cells $\Omega_j$ and $\Omega_i$,
set $K_{ji}$ contains reference numbers to cells that are used for computing the numerical flux function $F_{ji}$, $\vec{n}_{ji}$ is a unit
outward normal of the cell interface, $l_{ji}$ is length of the edge shared by triangles $j$ and $i$, $\bar{\rho}_k^n$ and $\bar{v}_k^n$
stands for time average of $\rho_k$ and $\vec{v}_k$ on $[]t_n,t_{n+1}]$. Notice that $K_{ji}=K_{ij}$, $l_{ji}=l_{ij}$, $\vec{n}_{ji}=-\vec{n}_{ij}$, $F_{ji}=-F_{ij}$.
The latter ensures mass conservation of the scheme (\ref{fv:FFSL4LA}).

\subsection{Integration by parts method for adjoint schemes}
\label{sec:Adjoint:IntegrByParts}
Here we consider adjoint scheme for ICON-FFSL with numerical flux function without limiters. In this case ICON-FFSL scheme (\ref{fv:FFSL4LA}) can
be equivalently written in the following form:

\begin{equation} \label{fv:FFSL:NonFluxForm}
\frac{\rho_j^{n+1} q_j^{n+1}-\rho_j^n q_j^n}{\Delta t} +
\frac{1}{| \Omega_j |} \sum_{i\in S_j} \alpha_{ji}(\{\bar{\rho}_k^n, \bar{v}_k^n\}_{k \in K_{ji}})q_{i}^n
=0, ~~S_j=\cup_{i\in I_j} K_{ji},
\end{equation}
where $S_j$ contains a reference number to cells on the stencil of the scheme, see ~\cref{fig:stencil}. The scheme (\ref{fv:FFSL:NonFluxForm})
is linear with respect to $\vec{q~}^n$ and it is the starting point for building the adjoint, i.e. for deriving the scheme that is consistent
with equation (\ref{eq:AdjointLinearAdvection}). In the latter equation the approximation of the time derivative and the right hand side
is straightforward. The spatial derivatives can be approximated, e.g. by means of rewriting second term in (\ref{fv:FFSL:NonFluxForm}) in matrix
form and then using transpose matrix. Equivalent but easy to use approach is using discrete analogue of integration
by parts formula, i.e. discrete analogue of the approach which was used for the derivation of (\ref{eq:AdjointLinearAdvection}).
For the convenience of further exposition we set: $\alpha_{ji}=0, ~\text{if}~ i \notin S_j$. Thus we have

\begin{multline} \label{fv:SumByParts:NonFluxForm}
\sum_{j} | \Omega_j | q^{*,n}_j \frac{1}{| \Omega_j |} \sum_{i\in S_j} \alpha_{ji}(\{\bar{\rho}_k^n, \bar{v}_k^n\}_{k \in K_{ji}})q_{i}^n
=
\sum_{j}  \sum_{i} \alpha_{ji}(\{\bar{\rho}_k^n, \bar{v}_k^n\}_{k \in K_{ji}})q_{i}^nq^{*,n}_j
= \\
 \sum_{i}  \sum_{j} \alpha_{ji}(\{\bar{\rho}_k^n, \bar{v}_k^n\}_{k \in K_{ji}})q_{i}^nq^{*,n}_j =
\sum_{i} | \Omega_i | q^{n}_i \frac{1}{| \Omega_i |} \sum_{j\in S^*_i} \alpha_{ji}(\{\bar{\rho}_k^n, \bar{v}_k^n\}_{k \in K_{ji}})q^{*,n}_{i},
\end{multline}
where $ S^*_i$ is set of indexes for which $\alpha_{ji} \neq 0$. On account of (\ref{fv:SumByParts:NonFluxForm}), (\ref{fv:FFSL:NonFluxForm})
and (\ref{eq:AdjointLinearAdvection}) adjoint numerical scheme writes:
\begin{equation} \label{fv:Adjoint:NonFluxForm}
\bar{\rho}_j^{n}\frac{ q^{*,n+1}_j- q^{*,n}_j}{\Delta t} +
\frac{1}{| \Omega_j |} \sum_{i\in S^*_j} \alpha_{ij}(\{\bar{\rho}_k^n, \bar{v}_k^n\}_{k \in K_{ij}})q^{*,n}_i
=\nabla_{q}f^o(q^{*,n}_j,q^{o,n}_j,t_n,\xi_j).
\end{equation}
Coefficients $\alpha_{ji}$ are given in the appendix \ref{A:Coef4Adjoint}.

\subsection{Artificial source term method for adjoint schemes}
\label{sec:Adjoint:ArtSource}
In this section we introduce a simple and efficient new method for developing the adjoint scheme and we apply it for building the
adjoint based on ICON-FFSL. The method can be used when the numerical flux function is  limited with or without flux limiters.
Our starting point is ICON-FFSL in flux form (\ref{fv:FFSL4LA}). Assuming the scheme is consistent with equation (\ref{eq:LinearAdvection})
the goal is constructing a numerical scheme that is consistent with the equation (\ref{eq:AdjointLinearAdvection}). The latter can be equivalently written
\begin{equation}\label{eq:LA+ArtSource}
\rho \frac{\partial q^*}{\partial t}+
\nabla\cdot(\rho q^*\vec{v})=
q^*\nabla\cdot(\rho \vec{v}) +
\nabla_{q}f^o(q,q^o,t,\xi),
\end{equation}
where $\nabla_{q}f^o(q,q^o,t,\xi)$ is defined by (\ref{eq:Nabla_qf^o}).
We will use this equivalent formulation of the adjoint equation as starting point of our method. The equation in the
form  (\ref{eq:LA+ArtSource}) is a good choice because we can effortlessly reuse discretization schemes developed
for the linear advection equation in conservative form (\ref{eq:LinearAdvection}). In particular, the spatial discretization
of ICON-FFSL (\ref{fv:FFSL4LA}) can be reused for the second and first term respectively, both in the left hand side and in the
right hand side of the equation (\ref{eq:LA+ArtSource}). As a result we have the following numerical scheme for the gradient computation
\begin{multline} \label{fv:adjointFFSL4LA}
\bar{\rho}_j^n \frac{q^{*,n+1}_j- q^{*,n}_j}{\Delta t}+\frac{1}{| \Omega_j |} \sum_{i\in I_j} F_{ji}(\{\bar{\rho}_j^n, \bar{v}_k^n,q^{*,n}_k\}_{k\in K_{ji}},\vec{n}_{ji}, l_{ji})= \\
\frac{1}{| \Omega_j |} \sum_{i\in I_j} F_{ji}(\{\bar{\rho}_j^n, \bar{v}_k^n,q^{*,n}_j\}_{k\in K_{ji}},\vec{n}_{ji}, l_{ji}) + \nabla_{q}f^o(q^{n}_j,q^{o,n}_j,t_n,\xi_j).
\end{multline}
Notice that the first sum in the numerical scheme (\ref{fv:adjointFFSL4LA}) is exactly the same as the sum in (\ref{fv:FFSL4LA}). The
second sum in (\ref{fv:adjointFFSL4LA}) is in principle the same with $q^{*,n}_j$ instead of $q^{*,n}_k$ that makes the computation of
the sum even easier. Therefore the particularly interesting feature of the numerical scheme (\ref{fv:adjointFFSL4LA}) is the
following: it makes easy reusing of the source code of the parent scheme, in our case ICON-FFSL scheme, and that dramatically reduces
time needed for adjoint code development, e.g. it took just a couple of days in our case.

\subsection{Properties of the adjoint scheme with artificial source term}

It is expected that the adjoint scheme constructed by artificial source term method will have good properties in the case if the parent
scheme has the same good properties as well. According to Lax's equivalence theorem the minimal set of properties of the schemes ensuring
convergence are consistency and stability. Therefore we will assume that parent scheme has these properties and then we prove that its
descendant adjoint scheme also enjoys the same properties.

\begin{theorem}[Consistency theorem]\label{Th:Consistency}
	Suppose $q,\rho, \vec{v}, q^*$ are sufficiently smooth and numerical scheme (\ref{fv:FFSL4LA}) is consistent with the linear
	advection equation (\ref{eq:LinearAdvection}) in the sense of local truncation error. Then adjoint scheme
	with artificial source term (\ref{fv:adjointFFSL4LA}) is also consistent with adjoint equation (\ref{eq:LA+ArtSource}).
\end{theorem}
\begin{proof}
	From the consistency requirement of numerical scheme (\ref{fv:FFSL4LA}) with equation  (\ref{eq:LinearAdvection}) we easily
	obtain consistency between corresponding terms in the scheme and equation. In particular it is clear that the expression
\begin{equation}\label{consistency1}
\frac{\rho_j^{n+1} q_j^{n+1}-\rho_j^n q_j^n}{\Delta t}
														~\text{is consistent with}~
\frac{\partial(\rho q)}{\partial t}
\end{equation}
at $(t_n, \xi_j)$ and therefore the expression
\begin{equation}\label{consistency2}
\frac{1}{| \Omega_j |} \sum_{i\in I_j} F_{ji}(\{\bar{\rho}_k^n, \bar{v}_k^n,q_k^n\}_{k\in K_{ji}},\vec{n}_{ji}, l_{ji})
~\text{is consistent with}~
\nabla\cdot(\rho q\vec{v}).
\end{equation}
at $(t_n, \xi_j)$. Assuming in the above expression $q_k=q_j$ means that $q$ is constant on the stencil and therefore
\begin{equation}\label{consistency3}
\frac{1}{| \Omega_j |} \sum_{i\in I_j} F_{ji}(\{\bar{\rho}_k^n, \bar{v}_k^n,q_j^n\}_{k\in K_{ji}},\vec{n}_{ji}, l_{ji})
~\text{is consistent with}~
q_j \nabla\cdot(\rho \vec{v}).
\end{equation}
Similarly with (\ref{consistency1}) it is evident that the first and fourth term in the numerical scheme with artificial source
term (\ref{fv:adjointFFSL4LA}) are consistent at $(t_n, \xi_j)$ with first and fourth term in the adjoint equation  (\ref{eq:LA+ArtSource}).
Substituting $q_j^n$ with $q_j^{*n}$ in (\ref{consistency2}) and (\ref{consistency3}) yields consistency at $(t_n, \xi_j)$ of the
second and third terms in the numerical scheme (\ref{fv:adjointFFSL4LA}) with  the second and third term in the equation  (\ref{eq:LA+ArtSource})
that concludes the proof.
\end{proof}
\begin{theorem}[Stability theorem]\label{Th:Stability}
Suppose
\begin{enumerate}
	\item $\rho$ is time independent.
	\item $\rho \geq \rho_{min} >0$.
	\item $\rho, \vec{v},F, q, \tilde{q}$ are smooth enough.
	\item $\vec{q^n}, \vec{\tilde{q}}^{n}$ are numerical solutions constructed by numerical scheme (\ref{fv:FFSL4LA}) with initial values $\vec{q_0}, \vec{\tilde{q}}_{0}$ respectively.
	\item Under Courant-Friedrichs-Levy (CFL) condition with CFL  number $CFL_0$  the numerical scheme (\ref{fv:FFSL4LA}) is stable in some norm:
	\begin{equation}\label{stab:FFSL4LA}
\|  \vec{q}^{~n+1} - \vec{\tilde{q}}^{~n+1} \| \leq (1+ C_0 \Delta t) \|  \vec{q}^{~n} - \vec{\tilde{q}}^{~n} \|, ~~0\leq n < N_t,
	\end{equation}
	constant $C_0$ is independent of $n$.
\end{enumerate}
If $1-5$ are valid the numerical scheme  with artificial source term (\ref{fv:adjointFFSL4LA}) is also stable under CFL condition with CFL
number $CFL_*=2CFL_0$ in the same norm,
 	\begin{equation}\label{stab:adjointFFSL4LA}
 	\|  \vec{q}^{~*,n} - \vec{\tilde{q}}^{~*,n} \| \leq C_1 \|  \vec{q}^{~*}_0 - \vec{\tilde{q}}^{~*}_{0} \|, ~~0<n\leq N_t,
 	\end{equation}
where 	$C_1$ is some constant independent of $n$.
\end{theorem}
\begin{proof}
	Numerical scheme (\ref{fv:adjointFFSL4LA}) equivalently writes:
\begin{multline} \label{fv:adjointFFSL4LA:split1}
\frac{1}{2}\bar{\rho}_j^n \frac{q^{*,1,n+1}_j- q^{*,n}_j}{\Delta t}+\frac{1}{| \Omega_j |} \sum_{i\in I_j} F_{ji}(\{\bar{\rho}_j^n, \bar{v}_k^n,q^{*,n}_k\}_{k\in K_{ji}},\vec{n}_{ji}, l_{ji})= 0,\\
\frac{1}{2}\bar{\rho}_j^n \frac{q^{*,2,n+1}_j- q^{*,n}_j}{\Delta t}=
\frac{1}{| \Omega_j |} \sum_{i\in I_j} F_{ji}(\{\bar{\rho}_j^n, \bar{v}_k^n,q^{*,n}_j\}_{k\in K_{ji}},\vec{n}_{ji}, l_{ji}) + \nabla_{q}f^o(q^{n}_j,q^{o,n}_j,t_n,\xi_j),\\
q^{*,n+1}_j=\frac{1}{2}(q^{*,1,n+1}_j+q^{*,2,n+1}_j).
\end{multline}
The first numerical scheme in (\ref{fv:adjointFFSL4LA:split1}) can be obtained from (\ref{fv:adjointFFSL4LA}) by replacing $\bar{\rho}_j^{n+1}$
with $\bar{\rho}_j^n$, see requirement $1$ of the theorem, and by replacing $\Delta t$ with $2\Delta t$.
Therefore, similar to (\ref{stab:FFSL4LA}), the estimate is valid for $q^{*,1,n+1}$ under the CFL condition with CFL number $CFL_*$ and we have:
\begin{equation}\label{stab:FFSL4LA:2}
\|  \vec{q}^{~*,1,n+1} - \vec{\tilde{q}}^{~*,1,n+1} \| \leq (1+ 2C_0 \Delta t) \|  \vec{q}^{~*,n} - \vec{\tilde{q}}^{~*,n} \|, ~~0\leq n < N_t.
\end{equation}
The first term in the right hand side of the second equation in (\ref{fv:adjointFFSL4LA:split1}) is consistent with $q_j \nabla\cdot(\rho \vec{v})$
according to theorem ~\ref{Th:Consistency}. Therefore, on account of the requirements of theorem ~\ref{Th:Stability},
we can assume that there exists a constant $C_2$ independent of $n$ such that the following inequality holds true:
\begin{multline}\label{fv:adjointFFSL4LA:split2}
 \frac{2}{\bar{\rho}_j^n|\Omega_j |} | \sum_{i\in I_j}[
F_{ji}(\{\bar{\rho}_j^n, \bar{v}_k^n,q^{*,n}_j\}_{k\in K_{ji}},\vec{n}_{ji}, l_{ji})-
F_{ji}(\{\bar{\rho}_j^n, \bar{v}_k^n,\tilde{q}^{*,n}_j\}_{k\in K_{ji}},\vec{n}_{ji}, l_{ji})] | \leq
2C_2 |  q_j^{~*,n} -  \tilde{q}_j ^{~*,n} |.
\end{multline}
The second term on the right hand side of the second equation in (\ref{fv:adjointFFSL4LA:split1}) does not depend on $q^*$. Therefore,
by analogy with explicit Euler time integration scheme on account of (\ref{fv:adjointFFSL4LA:split2}) we obtain the following stability inequality:
\begin{equation}\label{fv:adjointFFSL4LA:split3}
\|  \vec{q}^{~*,2,n+1} - \vec{\tilde{q}}^{~*,2,n+1} \| \leq (1+ 2C_2 \Delta t) \|  \vec{q}^{~*,n} - \vec{\tilde{q}}^{~*,n} \|, ~~0\leq n < N_t.
\end{equation}
Putting together the third equation in (\ref{fv:adjointFFSL4LA:split1}), (\ref{stab:adjointFFSL4LA}) and  (\ref{fv:adjointFFSL4LA:split3})
we obtain the following  stability estimate for $q*$ for one time step:
\begin{equation}\label{fv:adjointFFSL4LA:split4}
\|  \vec{q}^{~*,n+1} - \vec{\tilde{q}}^{~*,n+1} \| \leq (1+ [C_0+C_2] \Delta t) \|  \vec{q}^{~*,n} - \vec{\tilde{q}}^{~*,n} \|, ~~0\leq n < N_t.
\end{equation}
From (\ref{fv:adjointFFSL4LA:split4}) we arrive to (\ref{stab:adjointFFSL4LA}) with $C_1=\exp((C_0+C_2)T)$ which concludes the proof.
\end{proof}
\begin{remark}
	$CFL_*=2CFL_0$ is needed for the theoretical estimates only. In practice when calculating numerical tests the adjoint numerical
	scheme with artificial source term is stable under the same CFL condition as a parent scheme.
\end{remark}
\begin{remark}
	Theorem \ref{Th:Stability} deals with time independent $\rho$. The case with time dependent $\rho$  can be also be easily treated
	by a different equivalent formulation of the adjoint scheme with artificial source term. In particular a numerical scheme of the form
\begin{equation*}
\bar{\rho}_j^n \frac{q^{*,n+1}_j- q^{*,n}_j}{\Delta t}=RHS^n
\end{equation*}
can be equivalently written as the following two step scheme:
\begin{equation*}
 \frac{\bar{\rho}_j^{n+1}q^{*,1,n+1}_j- \bar{\rho}_j^nq^{*,n}_j}{\Delta t}=RHS^n, ~~
 \frac{q^{*,2,n+1}_j- q^{*,1,n+1}_j}{\Delta t}= \frac{1}{\bar{\rho}_j^n}\frac{\bar{\rho}_j^{n+1}-\bar{\rho}_j^n}{\Delta t} q^{*,2,n+1}_j.
\end{equation*}
The second scheme is an implicit Euler scheme and obtaining stability estimate is trivial. For obtaining stability estimate for the first scheme
theorem ~\ref{Th:Stability} can be invoked by similar arguments.
\end{remark}
\begin{remark}
	Theorems \ref{Th:Consistency} and \ref{Th:Stability} are sufficient for concluding convergence of artificial source term scheme since
	consistency and stability ensure convergence according to Lax equivalence theorem.
\end{remark}

\section{Numerical tests}
\label{sec:NumTests}
\subsection{Input for tests}

Here we consider two different types of tests: linear advection test cases and data assimilation test cases. Linear advection test
cases are used for comparing adjoint schemes with their parent scheme. Therefore numerical calculations are done with three schemes:
ICON-FFSL, adjoint scheme obtained by integration by parts formula,  and an adjoint scheme by introducing artificial source term.
Two different flux limiters are used together with ICON-FFSL and with the adjoint with artificial source term. In particular flux
limiters from ~\cite{zalesak1979fully}, \cite{schar1996synchronous} and from ~\cite{zalesak1979fully}, \cite{harris2011flux} are used.
Data assimilation test cases are used for comparing two adjoint schemes presented in this paper. In both test cases the true solution
is compared with numerical solutions computed with the above mentioned numerical schemes. Errors in numerical solutions are measured
in three different absolute and relative  norms. In  particular the following norms are used:
\begin{multline*}
l_{1,rel}=\frac{\sum_{i=1}^{N_c}|\Omega_i||q_i-q_i^{true}|}{\sum_{i=1}^{N_c}|\Omega_i||q_i^{true}|}, ~~~~
l_{1,abs}=\sum_{i=1}^{N_c}|q_i-q_i^{true}|, \\
l_{2,rel}=\frac{\sqrt{\sum_{i=1}^{N_c}|\Omega_i|{(q_i-q_i^{true})}^2}}{\sqrt{\sum_{i=1}^{N_c}|\Omega_i|(q_i^{true})^2}}, ~~~~l_{2,abs}=\sqrt{\sum_{i=1}^{N_c}{(q_i-q_i^{true})}^2},\\
l_{\infty,rel}=\frac{\text{max}_{i=\overline{1,N_c}}|q_i-q_i^{true}|}{\text{max}_{i=\overline{1,N_c}}|q_i^{true}|}, ~~~~
l_{\infty,abs}=\text{max}_{i=\overline{1,N_c}}|q_i-q_i^{true}|.
\end{multline*}

Notice that $l_1$ and $l_2$ measure integral characteristics and $l_{\infty}$ measures the maximum local deviation from the true
solution, but none of them measures either oscillations in the numerical solution directly or shape preservation. Therefore
where appropriate these two criteria will be also used.  We also consider magnitude of cost function and  norms of its gradient as most
important criteria for evaluating performance of numerical schemes in case of data assimilation tests cases that are given in the
subsection ~\ref{sec:assimilation}.

Standard test cases are available in scientific literature, such as ~\cite{nair2008moving}, \cite{nair2010class}, \cite{nair2002mass},
\cite{williamson1992standard}, \cite{zalesak1979fully}. We collected them in the appendix \ref{A:Tests:Wind,Shape},
in  ~\cref{tab:InitCond}, and in ~\cref{tab:Velocity} as collection of initial conditions and velocity fields.
In particular we define five different initial scalar fields with references in  ~\cref{tab:InitCond} and four different velocity
vector fields as given in ~\cref{tab:Velocity}. Combination of in initial scalar fields and velocity vector fields are used for
defining different test cases below.

\subsection{Linear advection test cases}
\subsubsection{ICON-FFSL, standard adjoint and adjoint with artificial source term }
First we compare the parent scheme ICON-FFSL with  two derived adjoint schemes presented in this paper. Notice that ICON-FFSL
is consistent with equation (\ref{eq:LinearAdvection}) and adjoint schemes are consistent with equation (\ref{eq:AdjointLinearAdvection}).
These two equations coincide with each other if $\rho=const$, $\nabla_q f^o=0$ and  $div~\vec{v}=0$. The latter condition is
satisfied by wind fields 1,2 and 4 in the ~\cref{tab:Velocity}. We also set $\rho=1$. With this selection of parameters we
can compare all three numerical schemes against each other on solutions of the problem (\ref{eq:LinearAdvection}),(\ref{eq:InitCondition}).
In particular we consider the following test problems
\begin{itemize}
	\item Solid body rotation with cosine bell.
	\item Solid body rotation with slotted cylinder.
	\item Deformational flow with two cosine bells.
	\item Deformational flow with two slotted cylinders.
\end{itemize}

Numerical results for all these schemes are summarized in  ~\cref{tab:SolidBodyRotationcos,tab:SolidBodyRotationcyl,tab:DeformationalFlow1cos,tab:DeformationalFlow1cyl}
and  on figures ~\cref{fig:CosineBellContours,fig:CosineBellCurve,fig:SlottedCylinderContours,fig:SlottedCylinderCurve}.
From the numerical results we conclude the following:
\begin{enumerate}
	\item All three schemes produce almost similar numerical results if flux limiters are not used in ICON-FFSL and in the adjoint
	      scheme with artificial source term.
	\item If flux limiters are used then ICON-FFSL scheme and  adjoint scheme with artificial source term produce practically similar
	      numerical results. Using different flux limiters has almost  no affect on accuracy of computation.
	\item The standard adjoint scheme produces larger errors compared to the ICON-FFSL scheme and the adjoint scheme with artificial
	      source term, if flux limiters are used, see ~\cref{fig:CosineBellCurve} and ~\cref{fig:SlottedCylinderCurve}.
	\item Tests are calculated on R2B4 grid with 20480 nodal points. If flux limiters are not applied then the maximum principle is
	      violated in almost in a half of the nodal points, overshoots and undershoots result in nonphysical negative numerical values.
	      ICON-FFSL and the adjoint with artificial source are free of these drawbacks when flux limiters are applied.
\end{enumerate}

\subsubsection{Standard adjoint and adjoint with artificial source term }
In this subsection the adjoint scheme with artificial source term with flux limiters is compared with the standard adjoint scheme.
Test problems featuring  deformational flow and moving vortices are considered. Numerical results are given
in ~\cref{tab:DeformationalFlow2cos,tab:DeformationalFlow2cyl,tab:MovingVortices}, and
in ~\cref{fig:DefFl1:2CosBellContour,fig:DefFl1:2CosBellCurve,fig:DefFl1:2SlottedCylinderContour,fig:DefFl1:2SlottedCylinderCurve,fig:DefFl2:2CosBellContour,fig:DefFl2:2CosBellCurve,fig:DefFl2:2SlottedCylinderContour,fig:DefFl2:2SlottedCylinderCurve,fig:MovingVortexContour,fig:MovingVortexAtT,fig:MovingVortexAtT/2}.
The analysis of the numerical results suggests that the adjoint scheme with artificial source term is much more accurate then
the standard adjoint scheme. In particular we observe the following properties:
\begin{enumerate}
	\item Standard adjoint scheme does not maintain the shape of the contours while the adjoint with artificial source is almost
	      indistinguishable from the exact solution, see e.g. ~\cref{fig:MovingVortexContour}.
	\item The standard adjoint scheme produces larger errors compared to adjoint scheme with artificial source term.
	\item Maximum principle is violated by the standard adjoint scheme resulting in overshoots and undershoots and negative
	      concentrations in almost half of nodal points of the ICON grid.
	\item For some test problems the adjoint scheme with artificial source term with flux limiter ~ \cite{zalesak1979fully},
	      \cite{schar1996synchronous} also gives negative concentrations in almost a third of nodal points of
	      the ICON grid, though negative values are of $10^{-8}$ magnitude at most. When flux limiter ~\cite{zalesak1979fully},
	      \cite{harris2011flux} is used then negative concentrations appear just in 14 nodal points and it's  magnitude
	      is very small - $10^{-15}$.
	\end{enumerate}

\subsection{Data assimilation test cases}
\label{sec:assimilation}
\subsubsection{Setup for tests}
\label{sec:setup}
In this subsection we study established test methods in the framework of passive tracer data assimilation.
We compare against each other standard adjoint, artificial source term adjoint with limiter from ~\cite{zalesak1979fully},
\cite{schar1996synchronous} and without limiter. Numerical results of advection tests given in previous subsections have
shown that standard adjoint and artificial source term adjoint without limiter produce almost similar results. Therefore
for some tests one case of them only will be considered; artificial source term method with limiter is considered in all
test cases.  The following three combinations of initial scalar fields and velocity vectors are selected:
\begin{enumerate}
\item Moving vortices that correspond to initial scalar field 3 from ~\cref{tab:InitCond} and to the vector field
4 from ~\cref{tab:Velocity}.
\item Deformational flow with cosine bells   that correspond to initial scalar field 4 from ~\cref{tab:InitCond} and
      to the vector field 3 from  ~\cref{tab:Velocity}.
\item Deformational flow with slotted cylinder that correspond to initial scalar field 5 from ~\cref{tab:InitCond} and
      to the vector field 3 from ~\cref{tab:Velocity}.
\end{enumerate}

Notice that for selected test problems with moving vortices $ div(\vec{v})=0$ and for the deformational
flow $ div(\vec{v}) \neq 0$.

Cost function ~(\ref{eq:CostFunc}) with discretization of the background term according to ~(\ref{eq:CostFunc:J^b}) and with
discretization of the observation term corresponding to ~(\ref{eq:CostFunc:J^o}) is used in numerical tests.
In this paper we do not study covariance models and the goal is development of adjoint solvers.
Therefore operators $K_b$ and $K_o$ are set to identity operators in numerical tests.  The cost function also contains the
functions $q_b$ and $q_o$ as input. Therefore for finalizing data assimilation related test problems we define observations
and background initial condition. For different test problems this is done differently, in particular, observations $q_o$
are defined as follows:
\begin{itemize}
	\item For numerical tests with moving vortices exact solution is known for an arbitrary time moment $t$ and
	therefore the exact solution is used for defining observations in selected points in space and in time,
	i.e.  $q_o(t_k, {\xi}_i)=u_{exact}(t_k,{\xi}_i), ~k=1,2,..N_T,~i\in \mathcal{O}$ .
	\item For numerical tests with deformational the flow exact solution is not known for an arbitrary time
	moment $t$ and therefore ICON-FFSL will be used for computing a reference solution which is used for
	defining observations in selected points in space
	and in time, i.e.  $q_o(t_k, x_i)=u_{ICON-FFSL}(t_k,x_i), ~k=1,2,..N_T,~i\in \mathcal{O}$.
\end{itemize}

The number of observation points and their location is different for different numerical tests and therefore they are presented
together with concrete numerical tests in the next subsection.

The background initial condition is constructed as an error in the true initial condition. This is done in the following way:
\begin{itemize}
	\item For numerical tests with moving vortices the background initial condition is defined as the
	true initial condition plus 	10 $\%$ of error in each nodal point.
	\item For numerical tests with deformational flow in the true initial condition the error is introduced
	in one half of computational domain only. For the second half of computational domain background initial
	condition and true initial condition coincide with each other. These domains are selected such that one
	cosine bell or slotted cylinder remains unchanged in the background initial condition. In those nodal
	points where the error is introduced and where the true initial condition is not zero background initial condition
	is defined as true initial condition plus 10 $\%$ of error. For the rest of nodal points of the half of
	computational domain 1$\%$ of maximum of the true initial condition is added.
\end{itemize}

The above procedure results in differences between true and background initial condition, details of which are
given in the ~\cref{tab:Datass_norms_weights_0_1}.

For minimizing the cost function the quasi-Newton method LBFGS ~\cite{LiuNocedal1989} is used here. The developed solvers
are applied for supplying gradient of the cost function in LBFGS. ICON-FFSL is used for forward
runs in order to supply needed input data to adjoint solvers. Magnitude of the cost function is used as most
suitable measure for comparing different adjoint solvers against each other. The difference between the true initial
condition and the one found by minimization procedure can also be used where appropriate.
Details on results of numerical tests are given in the next subsections.

\subsubsection{Convergence of iteration process}
Here we study behavior of developed adjoint solvers by means of performing 300 LBFGS iterations. We suppose 300 iterations
will be enough for demonstrating convergence of iteration process. Calculations are done on ICON R2B4 grid with
5120 observation points.
Observation points are distributed on the grid evenly. In case of 5120 observation points, every fourth point is observed.
For the test problem with moving vortices the initial cost function is around $2.5\cdot10^6$ and after
final iterations the cost is reduced to $131.9$ for standard adjoint, to $30.7$ for artificial source term
adjoint with limiter and to $19.2$ for artificial source term adjoint without limiter. For all three methods the progress of minimisation
is illustrated on ~\cref{fig:MV5120_cost_it300_1-10,fig:MV5120_cost_it300_11-20,fig:MV5120_cost_it300_295-300}.
For different number of iterations situation is different: for 10 iterations best reduction of the cost function is given
by standard adjoint, for 100 iterations best result is given by artificial source term adjoint and for 300 iterations best
result is given by artificial source term adjoint without limiter.
Oscillations visible in ~\cref{fig:MV5120_cost_it300_11-20} are due to the restart in LBFGS method.
Restart in minimisation is done when LBFGS algorithm reaches 5 attempts to find $\alpha$ step satisfying Wolfe
conditions ~\cite{LiuNocedal1989}. After the restart background condition is updated by the last
optimized initial condition and the minimisation process starts again.
The figure shows that after restart the cost the function increases first and then it decreases again. Similar behavior is observed in case
of all adjoint solvers considered.

Background and observation costs are given in
~\cref{fig:MV5120_obscost_it300_295-300} and ~\cref{fig:MV5120_bcost_it300_295-300}
respectively for the iterations 295-300.  We see that the standard adjoint solver is better at minimizing
the observation term, while the artificial source term adjoint solvers are much better in minimizing background
term of the cost function. This example shows the importance of the right adjoint solvers on the
minimization process: gradient defines descent direction and it is calculated using adjoint solver; therefore
different adjoint solvers can lead to different priorities, e.g. which term to be reduced in the cost function.

Standard adjoint and artificial source term adjoint  are numerically investigated on convergence  on test problem
of deformational flow with two cosine bells on ICON R2B4 grid with 5120 observation points.  Initial cost function is
around $3.5\cdot10^4$ and after 300 iterations the cost is reduced to $2.2$ for standard adjoint solver
and to $4\cdot10^{-2}$ for artificial source term adjoint solver. Both methods ensure convergence, the process is
illustrated on  ~\cref{fig:DF3cos5120_cost_it300_1-10,fig:DF3cos5120_cost_it300_11-20,fig:DF3cos5120_cost_it300_295-300}.
Before 30 iterations standard adjoint gives better results and after 30 iterations cost function corresponding to
artificial source term adjoint is smaller. In this test problem we do not observe a similar effect as in previous test
problem behavior in the quest for prioritisation of the observation or background term. Both methods minimize
background and observation terms of the cost function in the similar way.

These numerical results suggest that for the selected test problems all considered adjoint solvers ensure convergence of
iteration process, initial cost function is reduced approximately $10^5$ times after 300 iterations and
those adjoint solvers that give best reduction of the cost function after 300 iterations are different from those
which give better results for smaller number, e.g. 10 iterations.

\subsubsection{Impact of number of observations}
We study numerically three test problems given in the subsection \ref{sec:setup}.
Calculations are done on ICON R2B4 grid with 2560, 5120, 10240 and 20480 observation points for 50 LBFGS iterations.
Notice that R2B4 grid contains 20480 nodal points, i.e. we also consider the idealized case when observations are given in
all nodal points of the grid. Numerical results are given in ~\cref{fig:MV_R2B4_cost_it50} for tests with moving vortexes,
in the ~\cref{fig:DF3cos_R2B4_cost_it50} for tests with deformational flow and cosine bells,
in the  ~\cref{fig:DF3cyl_R2B4_cost_it50} for tests with deformational flow and slotted cylinders. These  results suggest
the following:
\begin{itemize}
	\item  In case of using standard adjoint solver the cost function monotonically increases together with the number of
	observation points for tests with deformational flow. For tests with moving vortices cost function first decreases
	for 5120 observation points and then it increases monotonically together with number of observation points.
	\item In case of using artificial source term adjoint solver with limiter the cost function decreases together with
	the number of observation points for tests with deformational flow and cosine bells. For a deformational flow
	with slotted cylinders the cost function increases for 20480 observation points though the growth is not as dramatic
	as in case of using standard adjoint solver. For tests with moving vortices cost function also decreases though we
	observe oscillations and peaks for 10240 observation points.
\end{itemize}

Numerical results suggest that the efficiency of iterations decreases when increasing number of observation points, i.e.
for the same number of iterations we get a larger cost function in case of more observation points. Using different adjoint
solvers influences minimization process differently and the adjoint solver with artificial source term offers more reliable
behavior when increasing number of observation points.

\subsubsection{Effect of mesh refinement}
Mesh refinement is a standard procedure for studying numerically convergence of numerical methods. As a result of mesh
refinement accuracy increases for convergent methods. Namely ICON FFSL and its descendant adjoint solvers should produce
more accurate solutions on refined meshes. Here we study the effect of mesh refinement on selected test problems in the
context of variational data assimilation. In numerical tests given in previous subsections the $R2B4$ grid was used.
Here we consider sequence of grids $R2B4,R2B5,R2B6,R2B7$ with 20480 observation points.
Corresponding cost functions are given in \cref{tab:Datass_mesh_ref}. We observe that the final cost
function increases together with the mesh refinement since we keep same number of observation points for all meshes.
Notice that for each test problem the starting value of cost function is almost the same for all meshes.
Though after 50 iterations we observe large difference in cost function evolution as impact of using different adjoint solvers.
Namely for the same number of iterations the artificial source term adjoint solver with limiter yields cost function values which are
from 8 to 660  times smaller than cost functions obtained with the adjoint solver without limiter.

\subsubsection{Manipulating wieghts of background and observation terms}
Weights can be used for prioritizing the background or observation term in the cost function. In numerical tests
considered in previous subsections  the weight was set to $0.5$ for both terms. Here in numerical tests  we gradually increase
the relative
weight of observation term up to $1$ and at the same time we decrease weight of background term down to $0$. Corresponding
cost functions on $R2B4$ grid for $50$ LBFGS iterations and $5120$ observation points are
given on ~\cref{fig:MV_weights,fig:DF3cos_weights,fig:DF3cyl_weights}. These figures suggest the following:
\begin{itemize}
\item For tests with moving vortices cost functions decrease when weight of observation term increases.
Cost function values calculated by artificial source term solver without limiter is
smaller then cost function with the same solver using a limiter, namely when background term disappears, i.e. weight of
observation term is $1$, we have  $J_{o,Art.S.NoLimiter}/J_{o,Art.S.WithLimiter}= 0.248777957$.
\item  For tests with deformational flow and cosine bells cost function valuess decrease when the weight of observation
term increases.  Cost function values calculated by artificial source term solver with limiter is smaller then cost
function produced by standard adjoint solver, we have  $J_{o,Std.Adjoint}/J_{o,Art.S.WithLimiter}= 1.878180425$.
\item For tests with deformational flow and slotted cylinders cost functions increase together with weight of observation
term. If the standard adjoint solver is used the growth of the cost function
is faster and  its value is also bigger compared to the case when the artificial source term adjoint with limiter
is used, $J_{o,Std.Adjoint}/J_{o,Art.S.WithLimiter}= 58.210774364$.
\end{itemize}
In case when the background term disappears in the cost function the true initial scalar field can be compared with those
found by data assimilation. The results are given in the
 ~\cref{tab:Datass_norms_weights_0_1}. They suggest the following
\begin{itemize}
	\item For tests with moving vortices some error norms are approximately two times smaller in case of using
	artificial source term adjoint without limiter compared to the same solver with limiter.
	\item For tests with deformational flow and cosine bells the standard adjoint solver gives better result compared
	to artificial source term adjoint solver with limiter,
	namely  in $l_1,l_{1,rel},l_2,l_{2,rel}$ norms error is slightly smaller and for $l_{\infty},l_{\infty,rel}$ norms
	the error is approximately 3 times smaller.
	\item For tests with deformational flow and slotted cylinders using the adjoint solver with artificial source term
	with limiter gives approximately 7 times smaller error norms compared
	to the case when standard adjoint solver is used.
\end{itemize}
Putting together ~\cref{fig:MV_weights,fig:DF3cos_weights,fig:DF3cyl_weights} and \cref{tab:Datass_norms_weights_0_1} we
conclude that value of cost function is a good
indicator for comparing adjoint solvers in the sense that if the cost function  is smaller, then the solution of the variational
data assimilation better approximates the true initial condition. We observe this behavior in case of
numerical tests considered with moving vortices and slotted cylinders. For the test problem with cosine bells this conclusion is
not valid for 5120 observations and it is true for 20480 observations.

\section{Conclusions}
We considered two approaches for developing adjoints of ICON-FFSL schemes for the linear advection equation in the
ICON model. The first approach is standard and it is used for building the adjoint of ICON-FFSL scheme without flux
limiters. Another approach is new and it is used for constructing the adjoint scheme of ICON-FFSL with or without
flux limiters. The new approach is based on rewriting the adjoint equation in flux form using an artificial source term and
then using the same discretization method as its parent ICON-FFSL scheme.  Because of this feature the development
of adjoint solver is easy and very fast if,as it is typically the case, the parent forward solver is available.
Another advantage is that flux limiters can be applied in our artificial source term method in a straightforward way
thus ensuring
stability regardless of smoothness of the solution. Stability and consistency of the adjoint scheme is proved
under assumption that the parent scheme has these properties. ICON-FFSL and its decsendant adjoint solvers are compared
against each other on a collection of standard advection test cases and variational data assimilation test cases
developed here. For some test cases when the solution is smooth, the adjoint solvers without limiter can provide competitive
and even smaller errors. For other and more realistic of test cases the artificial source term solver with limiter is the
clear winner ensuring reduction of the cost function even in cases where adjoint solvers without limiter fail to
perform minimization.

\section{Appendicies}
   \input{apendadjoint4icon_advection_test}
   \input{advection_numerical}
   \input{assimilation_numerical}
   \input{apendadjoint4icon_adjoint}

\bibliography{adjoint4icon_arXiv}

\end{document}

config: arxiv


TITLE: Product Sigma Algebra of two sets
QUESTION [1 upvotes]: Say I have two sigma algebras: $ \sigma(C) = \{E,\{\},\{1,2\},\{3,4,5\},\{1,2,3,4,5\},\{6\},\{3,4,5,6\},\{1,2,6\}\} $
where $E= \{1,2,3,4,5,6\}$
and $\sigma(Y) = \{E,\{\},\{1,3,5\},\{2,4,6\}\}$
How would one go about finding $\sigma(C) ⊗\sigma(Y) ?$
I know $\sigma(C) ⊗\sigma(Y) = \sigma(\{A\times B :  A\in \sigma(C) , B\in \sigma(Y)   \})$ is intuitively the smallest sigma algebra generated by the cartesian product.
Would this product algebra just give a set containing every tuple of sets from both sigma algebras eg: $(\{1,2\}\times\{1,3,5\})$?
Thanks.

REPLY [1 votes]: Answer to your second question:
These sets will belong to the product $\sigma$-algebra, but not every set in it will be such a set.
Realize that the complement of e.g. $\{1,2\}\times\{1,3,5\}$ is not a set that can be written as $A\times B$.
As you know $\sigma$-algebras are closed under complementation.

Answer to your first question:
In this special case there are partitions $\mathcal P=\{\{1,2\},\{3,4,5\},\{6\}\}$ and $\mathcal Q:=\{\{1,3,5\},\{2,4,6\}\}$ on set $E$. A set is element of the $\sigma(C)$  if and only if it can be written as a finite union of elements of $\mathcal P$. This union is also allowed to be empty. Likewise a set is element of the $\sigma(Y)$  if and only if it can be written as a finite union of elements of $\mathcal Q$.
This partition induces a partition $\mathcal R:=\mathcal P\times\mathcal Q=\{A\times B\mid A\in\mathcal P,B\in\mathcal Q\}$ on set $E\times E$, and a set is element of $\sigma(C)\otimes\sigma(Y)$ if it can be written as a finite union of elements of $\mathcal R$.

Note: because set $E$ is finite the class of algebras on it coincides with the class of $\sigma$-algebras on it.


text:
\begin{document}
\title{The Mean Interference-to-Signal Ratio and its \\Key Role in Cellular and Amorphous Networks}
\author{Martin Haenggi \\Department of Electrical Engineering\\{University of Notre Dame}
\thanks{Manuscript date \today. The support of the NSF (grants CNS 1016742 and CCF 1216407) is gratefully acknowledged.}
}

\maketitle
\begin{abstract}
We introduce a simple yet powerful and versatile analytical framework to approximate
the SIR distribution in the downlink of cellular systems. It is based on the
{\em mean interference-to-signal ratio} and yields the horizontal gap (SIR gain) between
the SIR distribution in question and a reference SIR distribution.
As applications, we determine the SIR gain for base station silencing, cooperation,
and lattice deployment over a baseline architecture that is based on a Poisson deployment
of base stations and strongest-base station association. The applications demonstrate that
the proposed approach unifies several recent results and provides a convenient framework
for the analysis and comparison of future network architectures and transmission schemes,
including {\em amorphous networks} where a user is served by multiple base stations and,
consequently, (hard) cell association becomes obsolete.
\end{abstract}
\section{Introduction}
\subsection{Motivation and contribution}
The SIR distribution is a key metric in interference-limited wireless systems.
Due to high capacity demands and limited spectrum, current- and next-generation
cellular systems adopt aggressive frequency reuse schemes, which makes
interference the main performance-limiting factor.
To overcome coverage and capacity problems due to interference, many
sophisticated transmission schemes, including base station cooperation and
silencing, successive interference cancellation, multi-user MIMO, and multi-tier
architectures have recently been proposed. However, a simple evaluation and
comparison of their effect on the SIR distribution has been elusive.

In this paper, we propose a novel technique that provides tight approximations
of the SIR gain of advanced downlink architectures and cooperation
schemes over a baseline scheme. It is based on the mean
 interference-to-signal ratio (MISR), which is used to quantify
the {\em horizontal gap} between two SIR distributions.
To account for the spatial irregularity of current and future cellular system,
we use point process models  for the positions of the base stations (BSs)
\cite{net:Dhillon12jsac,net:ElSawy13tut}.

\subsection{The horizontal gap in the SIR distribution}
We focus on the complementary cumulative distribution (ccdf)
$\bar F_{\sir}(\theta)\triangleq\P(\sir>\theta)$ of the SIR\footnote{The ccdf is often referred
to as the transmission success probability, while its complement, the cdf, is the
outage probability.}. There are two ways to compare
SIR distributions, vertically or horizontally, see \figref{fig:generic_sir} for an
illustration. Using the vertical gap,
\ie, the gain in the coverage probability, has several disadvantages:
(1) it depends strongly on the value of $\theta$ where it is evaluated; (2) it is often unclear whether the gain is measured in absolute or relative terms (for example, at -10 dB, the gap is 0.058, or 6.4\%;
at 0 dB, the gap is 0.22, or 39\%, and
at 20 dB, the gap is 0.05, or 78\%---and sometimes even the absolute gain is expressed in percentages);
(3) the gain also depends heavily on the path loss law and fading models.

In contrast, the horizontal gap (SIR gain)
is often quite insensitive to the probability where it is evaluated and the path loss models.
Formally, the gap between the distributions of $\sir_1$ and $\sir_2$ is defined as
\begin{equation}
G(p)\triangleq\frac{\bar F_{\sir_2}^{-1}(p)}{\bar F_{\sir_1}^{-1}(p)},\quad p\in(0,1),
\label{gp}
\end{equation}
where $\bar F_{\sir}^{-1}$ is the inverse of the ccdf of the $\sir$, and $p$ is the target
success probability.
In \figref{fig:generic_sir}, for example, $G(p)=5 $dB, irrespective of $p$.
In the following, we will illustrate that this behavior is commonly observed.
We also define the asymptotic gain (whenever the limit exists) as
\begin{equation}
G\triangleq G(1)=\lim_{p\to 1} G(p).
\label{gp2}
\end{equation}

\begin{figure}
\centerline{\epsfig{file=generic_sir,width=\figwidth}}
\caption{Example SIR ccdfs for a baseline and an improved scheme.
The vertical gap between the distributions depends
strongly on the value of $\theta$ where it is evaluated, while the horizontal gap is almost
constant.}
\label{fig:generic_sir}
\end{figure}

\section{The Mean Interference-to-Signal Ratio and its Significance}
\subsection{Definition}
We first define the interference-to-average-signal ratio $\isr$.
\begin{definition}[$\isr$]
The {\em interference-to-average-signal ratio} $\isr$ is defined as
\[ \isr\triangleq\frac{I}{\E_h(S)} ,\]
where $I$ is the sum power of all interferers and
$\bar S=\E_h(S)$ is the signal power averaged over the fading.
Its mean is denoted by $\misr\triangleq \E(\isr)$.
\end{definition}
The bar over the $S$ in the $\isr$ indicates averaging over the fading.
The $\isr$ is a random variable due to the random positions of the BSs relative to the
typical user.
For the following discussion, we assume
a power path loss law $\ell(r)=r^{-\alpha}$ with a path loss exponent $\alpha$
and (power) fading with unit mean, \ie, for all fading random variables,
$\E(h)=1$.  We also assume that the desired signal comes from a single BS at
distance $R$, while the interferers are located at distances $R_k$
 and their transmit powers (relative to the
one of the serving BS) are $P_k$. In this case, the $\isr$ is given by
\[ \isr=R^\alpha \sum_{k\in\mathcal{I}} h_kP_kR_k^{-\alpha}, \]
where $\mathcal{I}$ is the index set of the interferers and $h_k$ denotes the
channel (power) gain.
The mean follows as
\[  \misr\triangleq \E(\isr)=\sum_{k\in\mathcal{I}} P_k\E \left(\frac {R^\alpha}{R_k^\alpha}\right).\]
So the $\misr$ is a function of the {\em distance ratios} $R_k/R$ between the desired and interfering base stations,
scaled by the relative transmit powers.

\subsection{The asymptotic gap for Rayleigh fading}
The SIR distribution can be expressed using the $\isr$ as
\[ F_{\sir}(\theta)=\P(h\bar S<\theta I)=\P(h < \theta\, \isr). \]
For exponential $h$ and $\theta\to 0$,
\[ \P(h< \theta\, \isr\mid \isr )\sim \theta\, \isr ,\]
thus
\[  \P(h< \theta\,\isr)\sim \theta\, \E(\isr) .\]
So $F_{\sir}(\theta)\sim \theta \,\misr$, and
$\bar F_{\sir}^{-1}(p)\sim (1-p)/\misr$, $p\to 1$.
Consequently, the asymptotic gain between two SIR ccdfs
\eqref{gp2} can be expressed as
\begin{equation}
  G= \frac{\misr_1}{\misr_2},
  \label{asym_gain}
\end{equation}
and if it is finite, we have $\bar F_{\sir_1}(\theta)\sim \bar F_{\sir_2}(\theta/G)$,
$\theta\to 0$.

We will demonstrate in the
next section that this relationship
provides an accurate approximation
for the gain also at non-vanishing values of $\theta$, \ie, that
$\bar F_{\sir}(\theta)\approx \bar F_{\sir}(\theta/G)$ for all
practical values of $\theta$.

Other types of fading will be discussed in Sec.~\ref{sec:naka}.

\subsection{The HIP model and the baseline MISR}
The homogeneous independent Poisson (HIP) model was first introduced as a model
for cellular networks in \cite{net:Dhillon12jsac}.
\begin{definition}[HIP Model]
A {\em homogeneous independent Poisson} (HIP) model with $n$ tiers
consists of $n$ independent Poisson point processes (PPPs) $\Phi_k\subset\R^2$
with intensities $\lambda_k$, $k\in[n]$ and power levels $P_k$.
$\Phi_k$ is the set of locations of the base stations of the $k$-th tier.
\end{definition}
{\em Remarks:}
\begin{itemize}
\item Alternatively, the HIP model can be defined as follows:
Let $\lambda=\sum_{k\in[n]}\lambda_k$.
Starting with a homogeneous Poisson point process (PPP) $\Phi\subset\R^2$ of intensity $\lambda$,
randomly assign each point $x\in\Phi$ to one of the $k$ tiers $\Phi_k$,
where $\P(x\in\Phi_k)=\lambda_k/\lambda$, independently for each $x$.
\item
 Although the HIP model is a model for a heterogeneous network, we call it
{\em homogeneous}, since it is based on homogeneous PPPs, \ie, the base stations form a spatially homogeneous point process.
\item
The HIP model is doubly independent, since it exhibits neither intra-tier nor inter-tier dependence.
This makes it highly tractable but also makes it less accurate in situations where base stations
are deployed in a repulsive fashion (\ie, with a certain minimum distance) \cite{net:Deng14twc} or
where base stations of different tiers are not placed independently \cite{net:Deng14globecom}.
\item Quite remarkably, for the power path loss law with Rayleigh fading and with strongest-BS
association (on average, \ie, not considering small-scale fading),
the SIR distribution for the HIP model does not depend on the number
of tiers $n$, their densities $\lambda_k$, or their power levels $P_k$ \cite{net:Nigam14tcom}.
For $\alpha=4$, the SIR distribution is given by the extremely simple expression
\begin{equation}
 \bar F_{\sir[4]}(\theta)=\frac1{1+\sqrt\theta\arctan\sqrt\theta} .
 \label{arctan}
 \end{equation}
Also, $\E(S)=\infty$ and $\E(\sir)=\infty$ for all values of $\alpha$
due to the proximity of the nearest BS, and  $\E(I)=\infty$ for $\alpha\geq 4$.
\end{itemize}
Due to its tractability, the HIP model with strongest-base station association
is the perfect candidate for a baseline model against which the gains of other schemes
can be measured. Since the SIR distribution does not depend on the density or number of
tiers, we use a single-tier model in the following to calculate the MISR for the HIP model.

Let $R_k$ be the distance from the typical user to the $k$-th nearest BS.
Its distribution is given by \cite{net:Haenggi05tit}. The
distribution of the distance ratio $\nu_k=R_1/R_k$ is \cite[Lemma 3]{net:Zhang14twc}
\[ F_{\nu_k}(x)=1-(1-x^2)^{k-1} ,\quad x\in [0,1],\]
and the
$\alpha$-th moments are
\begin{equation}
   \E(\nu_k^\alpha)=\frac{\Gamma(1+\alpha/2)\Gamma(k)}{\Gamma(k+\alpha/2)} .
   \label{nu_alpha}
 \end{equation}
For equal powers $P_k\equiv 1$, the MISR follows as\footnote{
The parameter $\kappa^{\rm PPP}$ calculated as the limit
$\lim_{\theta\to 0} F_{\sir}(\theta)/\theta$ in \cite[Cor.~1]{net:Guo14tcom} is identical to the MISR for the PPP.}
\begin{equation}
 \E(\isr)=\sum_{k=2}^\infty \frac{\Gamma(1+\alpha/2)\Gamma(k)}{\Gamma(k+\alpha/2)}=\frac{2}{\alpha-2},\quad\alpha>2.
 \label{misr}
\end{equation}
For $\alpha=4$, $\misr=1$, which implies $F_{\sir}(\theta)\sim \theta$,
$\theta\to 0$.

\section{Applications}
\subsection{Base station silencing}
We consider the (single-tier) HIP model and
let $\isr^{(!n)}$ denote the $\isr$ if the $n$ strongest interfering BSs (on average) are silenced,
and all BSs transmit at the same power.

If the nearest interfering BS is silenced, the MISR is obtained by subtracting
$\E(\nu_2^\alpha)$ from \eqref{misr}, which yields
\[ \E(\isr^{(!1)})=\frac{2}{\alpha-2}-\frac{2}{\alpha+2}  =\frac{8}{\alpha^2-4} .\]

For general $n$,
\[ \E(\isr^{(!n)})=\frac{2\Gamma(1+\alpha/2)}{\alpha-2} \,\frac{\Gamma(n+2)}{\Gamma(n+1+\alpha/2)}.\]
The same result has been obtained in \cite[Prop.~1]{net:Zhang14twc} by calculating the limiting outage
probability $F_{\sir}(\theta)$ as $\theta\to 0$.

For $\alpha=4$, $\E(\isr^{(!n)})=\frac{2}{n+2}$, and
the asymptotic gain per \eqref{asym_gain} is simply
\[ G_{\rm silence[4]}=\frac{1}{\E(\isr^{(!n)})}=1+\frac n2 .\]

\figref{fig:bs_silence} shows the SIR distributions for the HIP model without silencing,
for the HIP model with silencing of one BS, and the MISR-based approximation.
The approximation is tight for success probabilities above $3/4$; after that, it is
pessimistic.
\begin{figure}
\centerline{\epsfig{file=bs_silence_a4,width=\figwidth}}
\caption{Gain from silencing one base station in HIP model for $\alpha=4$.}
\label{fig:bs_silence}
\end{figure}

\subsection{Base station cooperation for worst-case users}
We focus on {\em worst-case users} in the single-tier HIP model, which are the ones located at the vertices of the Voronoi
tessellation \cite{net:Jung13cl,net:Nigam14tcom}. These locations are marked by $\times$ in
\figref{fig:voro_worst}, and the SIR ccdf is denoted as $\bar F_{\sir}^\times$ accordingly.
Worst-case users are at a significant disadvantage if they are served by a single BS since
they have two other BSs at the same distance\footnote{In contrast to what is claimed in
\cite{net:Jung13cl}, these worst-case locations do {\em not} form a PPP.}.

With (non-coherent) joint transmission\footnote{The amplitudes of the three signals are adding up, and the combined received signal is still subject to Rayleigh fading.} from the 3 equidistant BSs and $\alpha=4$,
the ccdf follows from \cite[Thm.~2]{net:Nigam14tcom} as
\begin{equation}
   \bar F_{\sir[4]}^{{\times,{\rm coop}}}(\theta)=\bar F_{\sir[4]}^{{2}}(\theta/{3})=\Big(1+\sqrt{\theta/{3}}\arctan(\sqrt{\theta/{3}})\Big)^{-2}.
   \label{coop4}
 \end{equation}
where $\bar F_{\sir[4]}$ is the SIR ccdf for the typical user in the HIP model given in \eqref{arctan}.
The factor of $3$ is due to the gain in signal power, while the
exponent of 2 is due to the larger distance of the nearest BS than in the case of the
typical user.

Without cooperation, $\misr^\times=2+4/(\alpha-2)$, and for
$n\in\{2,3\}$ BSs cooperating, it follows from \cite[Thm.~4]{net:Nigam14tcom} that
\[ \misr^\times_{n-{\rm coop}}=\frac{4+(3-n)(\alpha-2)}{n(\alpha-2)} .\]
So for $n=3$, the gain is
\[ G_{3-{\rm coop}}=\frac{{\rm MISR}^\times}{\rm MISR_{3-coop}^\times}=3+\frac32(\alpha-2) .\]
\figref{fig:worst} shows the
SIR distribution for worst-case users without cooperation, with cooperation from
the 3 nearest BSs, and the MISR-based approximation.

\begin{figure}
\centerline{\epsfig{file=voro_worst,width=.9\figwidth}}
\caption{Illustration of worst-case user locations. Base stations are marked by $\odot$, and the crosses $\times$ are the vertices of the Voronoi tessellation and mark the locations of
the worst-case users. These users have the same distance to the 3 nearest BSs.}
\label{fig:voro_worst}
\end{figure}

\begin{figure}
\centerline{\epsfig{file=worst_case_isr,width=\figwidth}}
\caption{SIR ccdf for worst-case users without cooperation and with 3-BS cooperation (from
\eqref{coop4}) and MISR-based approximation for $\alpha=4$. Here $G=6$ (7.8 dB).}
\label{fig:worst}
\end{figure}

\subsection{Non-Poisson deployment}
An SIR gain can also be obtained by deploying the BSs more regularly (repulsively) than a PPP.
This gain has been termed {\em deployment gain} in
\cite{net:Guo13twc,net:Guo14tcom}. Exact closed-form results for the SIR distribution for
non-Poisson deployments are impossible to derive. However, the MISR-based
approximation, relative to the HIP model, is fairly easy to evaluate and quite accurate.

Simulations show that the MISR of the square
lattice is quite exactly half of that of the PPP, irrespective of the path loss exponent, \ie,
the deployment gain is 3 dB.
\figref{fig:square} shows that the resulting approximation is accurate over a wide
range of $\theta$. As a result,
\[ \bar F_{\sir[4]}^{\rm sq}(\theta) \approx \bar F_{\sir[4]}(\theta/2) ,\]
where $\bar F_{\sir[4]}$ is given in \eqref{arctan}, is a very good analytical approximation
to the SIR distribution in the square lattice. For the triangular lattice (hexagonal cells),
the gain is slightly larger, about 3.4 dB, which is the maximum achievable.

\begin{figure}
\centerline{\epsfig{file=ppp_square_a4,width=\figwidth}}
\caption{SIR ccdf for the square lattice and MISR-based approximation for $\alpha=4$ (Rayleigh fading).}
\label{fig:square}
\end{figure}

\section{General fading and diversity}
\label{sec:naka}
So far we have discussed the case of Rayleigh fading. The MISR framework
easily extends to other types of fading or transmission schemes with diversity
(\eg, coherent BS cooperation, MIMO, retransmission). The diversity under interference (DUI) is defined as \cite{net:Haenggi13twc,net:Zhang14twc}

\[ d\triangleq\lim_{\theta\to 0} \frac{\log F_{\sir}(\theta)}{\log\theta} .\]

When comparing SIR
distributions, the two schemes need to provide the same diversity gain $d$ (\ie, have the same
asymptotic slope in the outage curve when plotted on a log-log scale),
otherwise the asymptotic gain may be undefined or not helpful in approximating one ccdf
with the other.

For example, if the fading distribution satisfies $F_h(x)\sim a x^m$, $x\to 0$, (as, \eg, in Nakagami-$m$ fading)
\[ F_{\sir}(\theta)\sim a \theta^m \E(\isr^m), \]
the diversity order is $m$---if the $m$-th moment of the $\isr$ is finite.\footnote{
For the PPP, it can be shown that all moments of the $\isr$ are finite. Whether it holds
for {\em all} stationary point processes is under investigation.}
The asymptotic gain follows as
\[ G^{(m)}= \left(\frac{\E(\isr_1^m)}{\E(\isr_2^m)}\right)^{1/m}
\approx G^{(1)},\]
where the approximation by $G^{(1)}$ holds since the factor $\E(\isr^m)^{1/m}/\misr$ is
about the same for both schemes and thus cancels approximately. This is illustrated in
\figref{fig:square_naka2} for Nakagami-2 fading and the square lattice. The shift by
$G^{(1)}=3$ dB still yields a very good approximation.
\begin{figure}
\centerline{\epsfig{file=ppp_square_naka2,width=\figwidth}}
\caption{SIR ccdf of square lattice and MISR-based approximation for $\alpha=4$ and Nakagami-2 fading.}
\label{fig:square_naka2}
\end{figure}

\section{Conclusions}
The SIR distributions of two transmission schemes or deployments in cellular networks that provide
the same diversity gain are, asymptotically, horizontally shifted version of each other, and the asymptotic
gap between them is quantified by the ratio of the MISRs of the two schemes. We demonstrated
that this asymptotic gap also provides a good approximation for the gap at finite $\theta$.

If the MISR cannot be calculated analytically, it is relatively easy to determine by simulation since
it only depends on the BS and user locations and the transmit power levels, but not on the fading.
Due to its tractability, the HIP model is the prime candidate as a baseline model against which
other schemes can be evaluated. So even if it is not be accurate in all situations, it serves the
important role of a reference model.

We anticipate that future networks will not be based on a strict cellular architecture
but will become {\em amorphous} due to cooperation between BSs at different levels, relays,
and distributed antenna systems. Since an exact analytical evaluation for the SIR distribution
for these sophisticated and cognitive architectures seems hopeless, we believe that the proposed
MISR framework will  play an important role in the analysis of such emerging  {\em amorphous networks}.
\vspace*{-1mm}

 \bibliographystyle{IEEEtr}