\chapter{Optimal design in $CO_{2}$ injection monitoring}
\label{Ch6}

As we have mentioned in the introduction, in order to reduce the Greenhouse effect, companies and organizations have
piped the liquefied $CO_{2}$ into underground storage, which brings up the need for the research into monitoring techniques
to verify that the $CO_{2}$ remains effectively trapped underground.

It has been shown that time-lapse seismic imaging is an effective technique for underground $CO_{2}$ monitoring 
\cite{Arts2002, Gritto2004}. Borehole seismic techniques can image the change in seismic velocity caused 
by the movement of $CO_{2}$ \cite{Daley2007}. 

However, due to the huge expense in deployment of the monitoring sensors, it is imperative for us to find a cost-effective
arrangement for placing the sensors while the quality of the images is still guaranteed. Therefore, in this chapter, we will 
present an approach to find the optimal design for $CO_{2}$ injection monitoring.

\section{Crosswell array constraints}
\label{SecConst}

We begin this chapter with a brief description of the basic structure of the experimental setting used in $CO_{2}$ injection 
monitoring. Figure \ref{pic:Well} provides a schematic of the crosswell tomography experiment.

\begin{figure}[ht!]
 	\centering
 	\includegraphics[scale=0.5]{Well.jpg}
 	\caption{\label{pic:Well}A schematic of the crosswell tomography experiment.}
\end{figure}

This schematic has brought us various crosswell array constraints that have to be taken into account when designing the 
$CO_{2}$ injection monitoring. We mention some major constraints in the following:

\begin{itemize}
\item{\bf Fixed packers:}
The seismic sensors are deployed on production tubing, both above and below a packer in the observation well.
The locations of packers in both observation wells are fixed.
\item{\bf Jewelry on tubing:}
The pressure compensator is at the bottom of each receiver array and the pre-amplifier is at the top. They both
have their individual sizes.
\item{\bf Well boundary:}
The possible positions of the sensors are bounded from both above and below.
\item{\bf Offset between the wells:}
The distance between the two wells is also fixed.
\item{\bf Equally spaced receivers:}
Although sources can be placed almost anywhere, the two groups of receivers must have equal distances between 
consecutive pairs (within the groups). 
\item{\bf Length of sensors:}
Each sensor has its individual length which cannot be ignored.
\item{\bf Maximum ray length:}
The maximum length of any possible ray path is not allowed to exceed certain limit.
\end{itemize}

All these constraints bring significant difficulties in the optimal design.

\section{Mathematical framework}

In this section, we discuss the mathematical framework of the $CO_{2}$ injection monitoring problem. Before the
formulation is proposed, we will describe a bit about the generation of the forward modeling operator.

\subsection{The forward modeling operator}

If the underground $CO_{2}$ saturation and/or plume thickness increase along a given raypath, the traveltime of a ray would 
decrease. From the time-distance data, a tomographic reconstruction of the wave speed can be carried out by
using sensitivity kernels obtained in the ray approximation.

When we discretize the region of interest, the integrated sensitivity kernel in a given grid cell becomes the
length of the ray path in that cell. Therefore, the quality of the imaging reconstruction depends on the accuracy 
of the sensitivity kernels employed.

We apply the Fresnel-Zone-Based Kernel \cite{Jensen2000} hence, the sensitivity function is formulated as
$$
S(x) = K A(s, x)A(x, r)\cos{(\frac{2\pi \Delta t(x)}{T})}\exp{(-(\frac{\alpha \Delta t(x)}{T/4})^{2})},
$$
where $S(x)$ is the sensitivity at the position vector $x$, $A$ represents the amplitude of the wave field
propagation and $K$ is a normalization constant. The parameter $\alpha$ controls the degree of cancellation
in Fresnel zones beyond the first, $\Delta t(x)$ stands for the delay time between the first arrival and 
the arrival of waves that have been scattered at $x$ and $T$ is the dominant period of the wave.

\subsection{Formulation of the inversion}

The $CO_{2}$ injection monitoring problem can be expressed by this inverse problem, 
$$
Gm + \epsilon = d,
$$
where $G$, is a matrix in which each row represents a tomographic ray path that connects a source to a receiver.
We calculate $G$ in a $2.5$ dimensional geometry designed to accommodate heterogeneous 2D structure and well 
trajectories with out-of-plane deviations. The solution to the eikonal equation is calculated using a finite-difference 
eikonal solver \cite{HansenMosegaard, Jensen99sensitivitykernels}. Our numerical experiments are based on the FAST package 
developed by Zelt \cite{Zelt98}. The vector $m$ is the geophysical model we would like to reconstruct, $\epsilon$ is the 
noise that is assumed to be Gaussian iid with 0 mean and standard deviation $\sigma$ and finally $d$ is the data
we have observed. 

Due to the existence of the complicated spatial constraints that we have raised in the previous section, we cannot hope to
control the quality of the reconstructed images by simply abandoning certain raypaths or adjusting the frequency with which
certain raypaths is chosen. Therefore, the optimal design that uses a weight matrix $W$ as was described in chapter \ref{Ch2}
will not be helpful. Therefore, we intend to explore a method to gradually modify the kernel matrix $G$ in each iteration
as we update the placement of sensors, thereby converging to an optimal kernel $G$ with a corresponding optimal 
placement of sensors.

Since geophysical inverse problems are usually ill-posed, regularization in general is needed. In this problem,
we apply the Tikhonov regularization framework as was described in chapter \ref{Ch4}. Hence, an estimated model is given by
$$
\hat{m} = {\rm arg} {\rm \min} \hspace{2mm} (d - Gm)^{\top}(d - Gm) + \alpha \parallel Lm \parallel^{2}.
$$ 
In this case, we set $L$ to be a discrete derivative operator to introduce smoothness to the recovered model.
The solution of the above problem is given by
$$
\hat{m} = (G^{\top}G + \alpha L^{\top}L)^{-1}G^{\top}d.
$$

In order to determine the quality of the estimates, we calculate the $\ell^{2}$ norm difference between the true model and the 
corresponding tomographic reconstructed one for a given geometry. In many cases, the desired information is usually only within
a particular zone of the whole model, therefore the quality measurement is only considered within the target zone:
\begin{equation}
\label{D}
D = \parallel Q\odot(\hat{m} - m_{true}) \parallel^{2}_{2},
\end{equation}
where $Q$ is a window function that is used to specify the zone that we are interested in.

The optimality criterion is defined by the average MSE that is given in the following
\begin{equation}
\label{MSECO2}
\begin{array}{lll}
{\rm MSE} & = & \alpha^{2} \parallel Q \odot (G^{\top}G + \alpha L^{\top}L)^{-1}L^{\top}Lm_{true} \parallel^{2}_{2} \\
		  & + & \sigma^{2} {\rm trace}[G(G^{\top}G + \alpha L^{\top}L)^{-1} \odot Q^{\top}Q \odot (G^{\top}G + \alpha L^{\top}L)^{-1}G^{\top}]. 
\end{array}
\end{equation}
In most real applications, the true model in the AMSE formulation above is not available.
In its place, we use a suite of reference models that are related to the model to be reconstructed and solve the problem
via a training approach.
\begin{equation}
\label{AMSECO2}
\begin{array}{lll}
{\rm AMSE} & = & \frac{1}{s} \displaystyle \sum_{j=1}^{s} \alpha^{2} \parallel Q \odot (G^{\top}G + \alpha L^{\top}L)^{-1}L^{\top}Lm_{true, j} \parallel^{2}_{2} \\
		  & + & \sigma^{2} {\rm trace}[G(G^{\top}G + \alpha L^{\top}L)^{-1} \odot Q^{\top}Q \odot (G^{\top}G + \alpha L^{\top}L)^{-1}G^{\top}], 
\end{array}
\end{equation}
where $s$ is the number of reference models.
 
Therefore, the main purpose of this work is to determine an optimal placement of sources and receivers in 
order to provide maximum information about the target geophysical structure at minimum financial cost under 
the spatial constraints. In particular, the best possible estimated model is obtained by minimizing the AMSE through 
optimization approaches and the optimal reconstructed model is obtained based on the resulting optimal design thereafter.

\section{Numerical optimization through DIRECT algorithm}

Obviously standard gradient-based optimization methods are not very effective for solving our problem due to the existence 
of those constraints on the solution domain. In stead, Ajo-Franklin \cite{AjoFranklin2008} has chosen two non gradient-based
optimization methods: the Nelder-Mead downhill simplex method \cite{Nelder} and the multilevel coordinate search algorithm
\cite{Huyer}. However, the work did not provide a way to find the optimal locations of sensors. Also the discussion of
location constraints was not presented.

\subsection{The DIRECT algorithm}

In order to accommodate those location constraints for placing sensors, we apply the DIRECT algorithm, which was proposed 
by Daniel E. Finkel \cite{Finkel}. The DIRECT optimization algorithm was first introduced in \cite{Jones} based on the 
Lipschitzian Optimization for solving difficult global optimization problems with bound constraints and a real-valued 
objective function. 

The algorithm converges to the global minimum of the objective functional without requiring its gradient. 
It samples points in rectangular domains by taking the midpoints of the searching spaces. This overcomes the 
shortcoming of the Lipschitzian Optimization when solving high dimensional problems. The DIRECT algorithm uses all
values in the current rectangle to determine if a region of the domain should be broken into sub-rectangles. 
Therefore, it does not require the estimate of the Lipschitz constant or the continuity of the objective functional \cite{Jones}.

In a very brief description, DIRECT algorithm starts with a unit hyper-cube and divides it into smaller hyper-rectangles with
the restriction that the division is only done along the longest dimension(s) of the hyper-rectangles. This ensures that 
the rectangles will be divided along all directions and no dimension will be ignored. 

For our application of the DIRECT algorithm, we only need to provide the evaluation of the objective functional and the 
upper and lower bounds of the variables we intend to estimate. Thus, the crosswell ray constraints can be easily satisfied 
if we carefully choose the variable bounds.

\subsection{Discussion of the constraints}

We now discuss details about finding the appropriate bounds for locations of sources and receivers. First, we need to 
determine the size of the kernel $G$ in the inverse problem, which represents all possible locations of sources and receivers.  

Based on the knowledge of the maximum ray length (we denote as $l_{\rm MaxRay}$) and the offset ($d_{\rm off}$) between 
the two wells, we obtain the maximum distance that a raypath can reach along the well direction 
$d_{\rm max} = \sqrt{(l_{\rm MaxRay})^{2} - (d_{\rm off})^{2}}$. 

In the source well, we assume that the lowest source above the packer can be placed right on top of the packer and the highest
source below the packer can be placed right under it. Similar assumptions work for the receiver well as well.
Also we denote $l_{\rm source}$ to be the length of the sources and $l_{\rm receiver}$ to be the length of the receivers. 

Thus, in order to guarantee that all rays from the sources are able to reach all receivers, we must not place any sensors in
areas that a ray cannot reach. Therefore, the upper bound on the lowest receiver below the packer (we denote as $R_{\rm bottom}$)
is at: top of the packer in the source well - $l_{\rm source}/2$ + $d_{\rm max}$. Similarly, the lower bound on the highest 
receiver above the packer ($R_{\rm top}$) is at: bottom of the packer in the source well + $l_{\rm source}/2$ - $d_{\rm max}$. 
Following the same approach, the upper bound on the lowest source below the packer ($S_{\rm bottom}$) is determined to be at:
top of the packer in the receiver well - $l_{\rm attach}$ - $l_{\rm receiver}/2$ + $d_{\rm max}$ and the lower bound on 
the highest source above the packer ($S_{\rm top}$) is at: bottom of the packer in the receiver well + $l_{\rm attach}$ + 
$l_{\rm receiver}/2$ - $d_{\rm max}$, where $l_{\rm attach}$ is the length of the attached materials on receivers such as
the pressure compensator and the pre-amplifier.

The size of the kernel $G$, or in other words, the size of the space for all possible locations of sources and receivers
is decided by the above four bounds. 

\subsubsection{Boundaries for sensor location constraints}

Since the sources can be placed at any location within the space provided by the kernel $G$, we must determine individual 
bounds for each of them. We assume that the number of sources above the packer is $n_{\rm Stop}$ and the number below the
packer is $n_{\rm Sbottom}$. Similarly, we assign numbers to $n_{\rm Rtop}$ and $n_{\rm Rbottom}$ as well.

We start with the possible locations for the highest source. Its lower bound is simply $S_{\rm top}$ and 
its upper bound can be calculated as: top of the packer in the source well - $l_{\rm source}/2$ - ($n_{\rm Stop} - 1$)
$l_{\rm source}$. With respect to the second highest source, the lower bound is $S_{\rm top}$ + $l_{\rm source}$ because the space to hold the
highest source must be kept. Similarly, its upper bound is $l_{\rm source}$ more than the upper bound
of the highest source. Following the same routine, it is straightforward to obtain bounds for all sources.

For the receivers, since each group is equally spaced, we only need the four bounds on the top and bottom receiver of the 
group above the packer and below it. The lower bound for the top receiver in the top group is simply $R_{\rm top}$ and its 
upper bound is calculated as: top of the packer in the receiver well - $l_{\rm attach}$ - $l_{\rm receiver}/2$ - 
($n_{\rm Rtop} - 1$)$l_{\rm receiver}$. The other bounds can be determined in the same fashion. Finally, after translating 
the absolute bounds into the coordinate system of the kernel $G$, we put the bounds for all sources and receivers in a 
$n_{\rm Stop} + n_{\rm Sbottom} + n_{\rm Rtop} + n_{\rm Rbottom}$ by $2$ matrix.

\subsubsection{Special treatment for receivers}

Above, we have decided the general boundaries for placing sources and receivers. However, due to the complicated constraints
we have listed in section \ref{SecConst}, there are certain places inside the solution interval between upper and lower
bounds where we cannot place sensors. Hence, special treatments need to be considered when finding the optimal placement.
Since sources and receivers have different spatial constraints, different treatments need to be carried out in the two cases.

Let us start with the easy case: the receivers. We use the top group of receivers as an example. There are two things we need 
to check every time we update to a new arrangement of receivers. First, if the distance between the first and last receivers 
in this group is less than $l_{\rm receiver}(n_{\rm Rtop}-1)$, we stretch the whole group until this constraint is satisfied. 
Here, distance between two sensors means the distance between the centers of the two. Since, in our problem, we stretch the 
receivers downwards, we need to check the position of the last receiver. If it is below the top of the packer, we will move 
all receivers up until this conflict is resolved. Of course, it is possible to stretch the receivers upwards as well. However,
since the target zone is usually in the middle area between the two wells, it is natural for us to always try to gather
sensors close to the target zone. Therefore, for the top group of receivers, we stretch downwards.

A symmetric process is performed on the bottom group of receivers, for which we stretch the conflicting receivers upwards.
In this case, we need to check the position of the first receiver in this group. If it is above the bottom of the packer,
we will move all receivers down until this conflict is resolved.

\subsubsection{Special treatment for sources}
 
Things are more complicated for sources. If the DIRECT algorithm generates the distance between two sources to be 
less than $l_{\rm source}$, the resulting system is not feasible because these two sources will overlap each other.
In order to fix this problem, we suggest a perturbation process. Again, we use the top set of sources as an example. 
A symmetric process can be applied to the group below the packer.

If two (or more) sources are conflicting, we arbitrarily choose one and move it to the closest position, within the space 
between the packer and the ceiling of the working space, which does not conflict with another source. This process is repeated 
until no conflicts remain among all sources. 

Note that no source will be moved more than the distance of $(n_{\rm Stop}-1)(2l_{\rm source}-1)$. The proof of this claim is
given as well.
\begin{proposition}
\label{eachsourcemoving}
Each source moves at most $(n_{\rm Stop} - 1)(2l_{\rm source} - 1)$.
\end{proposition}
\proof
First note that the DIRECT algorithm will not place any sources in conflict with the packer. Thus, if $n_{\rm Stop} = 2$ 
and the two sources are in conflict, we may move one source to the opposite side of the other, thereby moving the source by at 
most $2l_{\rm source} - 1$.  This would be the case if source $1$ is located right next to the packer while source $2$ is 
at a distance of $l_{\rm source} - 1$ away.

If $n_{\rm Stop} > 2$, the largest possible distance that a source could move is clearly $(n_{\rm Stop} - 1)(2l_{\rm source} - 1)$,
which is the situation when the first source to be considered is still located right next to the packer, the second is 
$l_{\rm source} - 1$ away, the third is $2(l_{\rm source} - 1) + l_{\rm source}$ away and so on. In this situation, 
the $i^{th}$ source would be located $(i - 1)(l_{\rm source} - 1) + (i - 2)l_{\rm source}$ from the packer for all 
$i \geq 2$. \qed

Unfortunately, this process does not always yield a solution that is as close as possible to the output of the DIRECT algorithm, 
but we can easily bound the difference. We claim that the total movement of the sources is less than $n_{\rm Stop}^{2}l_{\rm source}$. The
proof is given in the following.
\begin{proposition}
The sum of the movements of all the sources is at most \\
$\frac{n_{\rm Stop}(n_{\rm Stop} - 1)(2l_{\rm source} - 1)}{2}$.
\end{proposition}
\proof
We prove this result by induction on $n_{\rm Stop}$.  First suppose $n_{\rm Stop} = 2$. Then, by 
Proposition~\ref{eachsourcemoving}, one source will move at most $2l_{\rm source} - 1$.  If one source moves, the other will clearly 
remain unmoved so the total movement is at most $2l_{\rm source} - 1 = \frac{n_{\rm Stop}(n_{\rm Stop} - 1)(2\l_{\rm source} - 1)}{2}$.

Now assume $n_{\rm Stop} > 2$ and suppose the result holds for fewer than $n_{\rm Stop}$ sources.  Suppose we have $n_{\rm Stop} - 1$ 
sources placed and we would like to add a new source.  By Proposition~\ref{eachsourcemoving}, this source could potentially 
have to move at most $(n_{\rm Stop} - 1)(2l_{\rm source} - 1)$.  The total of the other sources movements, by induction, is at 
most $\frac{(n_{\rm Stop} - 1)(n_{\rm Stop} - 2)(2l_{\rm source} - 1)}{2}$ for a total of

\begin{eqnarray*}
& ~ & \frac{(n_{\rm Stop} - 1)(n_{\rm Stop} - 2)(2l_{\rm source} - 1)}{2} + (n_{\rm Stop} - 1)(2l_{\rm source} - 1) \\
& = & \frac{n_{\rm Stop}(n_{\rm Stop} - 1)(2l_{\rm source} - 1)}{2},
\end{eqnarray*}
as desired. \qed

Of course, there are cases where other algorithms will produce results that are closer to the output than our method. However, 
the computational complexity of these algorithms is orders of magnitude larger than the complexity of our algorithm.

At each step in our correction process, we move a source one distance unit in both directions and then check its position 
relative to all the other sources along with the packer and the outer boundary of the well (either the surface of the ground 
if these sources are above the packer or the bottom of the well if they are below). Thus, for each unit moved, we check a 
total of $2(n_{\rm Stop} + 1)$ pairs for conflicts. This makes the computational complexity of our algorithm exactly
\begin{eqnarray*}
& ~ & 2(n_{\rm Stop} + 1)\frac{n_{\rm Stop}(n_{\rm Stop} - 1)(2l_{\rm source} - 1)}{2} \\
& = & n_{\rm Stop}(n_{\rm Stop} + 1)(n_{\rm Stop} - 1)(2l_{\rm source} - 1) \\
& = & O(n_{\rm Stop}^{3}l_{\rm source}).
\end{eqnarray*} 

There are two obvious limitations in the method that we have chosen. One drawback of our algorithm is that, in some cases, 
the order in which the sources are considered causes the movement to increase. Hence, one may consider all possible 
orderings of the sources and choose the ordering which minimizes the movement. Unfortunately, the number of such orderings 
is $n_{\rm Stop}!$ so this would greatly increase the computational complexity to $O(n_{\rm Stop}!n_{\rm Stop}^{3}l_{\rm source})$.
Also, consider the following extreme example. Suppose all sources start in the same position, right next to the packer.  
Call this example $A$. In this case, our correction process finds the best possible configuration regardless of the order 
in which the sources are considered. This example shows that considering all orderings of the sources would likely be a 
waste of effort.

Another drawback of our approach is that there are cases where moving one source at a time cannot produce the best possible 
result. For example, suppose $l_{\rm source} = 3$, source $1$ starts at point $0$, source $2$ starts at point $4$, source 
$3$ starts at point $5$ and source $4$ starts at point $9$ (see Figure~\ref{badconf}). Clearly sources $2$ and $3$ are in 
conflict but moving only one source at a time, as in our algorithm, cannot place these two in their optimal positions 
(clearly $3$ and $6$ respectively). Thus, a slightly more sophisticated algorithm would allow pairs of sources to move 
at the same time, but this would have a computational complexity of at least $O(n_{\rm Stop}^{4}l_{\rm source})$ since 
we would have to consider all possible pairs of sources which introduces an $n_{\rm Stop}^{2}$ term as the number of pairs 
where we considered only the $n_{\rm Stop}$ sources before. Furthermore, as seen in example $A$, this more complicated 
algorithm could not always provide better results than our approach. 

\begin{figure}[ht!]
 	\centering
 	\includegraphics[scale=0.9]{badsources.jpg}
 	\caption{\label{badconf}An example of source placement.}
\end{figure}

The process we propose is very fast even when the number of sources gets large.
Therefore, we feel that our repair process provides a sufficient balance between computational cost and the result accuracy.

\section{Numerical experiments}

For our experiment we have generated four time-lapse images of a $CO_{2}$ flood progressing through a permeable layer as 
is shown here. We use the top three images as our reference models and we want to recover the one on the bottom row.

\begin{figure}[ht!]
 	\centering
 	\includegraphics[scale=0.23]{b1.jpg}
 	\includegraphics[scale=0.23]{f1.jpg} 
 	\includegraphics[scale=0.23]{f3.jpg} \\
 	\includegraphics[scale=0.35]{f2mtrue.jpg}
 	\caption{\label{CO2true}Four time-lapse images of a $CO_{2}$ flood progressing through a permeable layer.
The top three images are our reference models and the goal is to recover the one on the bottom row.}
\end{figure}

We designed a crosswell continuous active-source seismic monitoring (CASSM) experiment in which data would be acquired 
continuously during injection along a set of fixed raypaths.

We give a list of parameters we have used in our experiments. We assume that the length of the attached materials on the 
receiver groups is $2$ feet. The maximum length of raypaths is $150$ feet and the offset between the two wells is $80$ feet.
Hence, we obtain that the maximum distance that a raypath can reach along the well direction is $d_{\rm max} = 126$ feet. 
The length of both sources and receivers is $2$ feet. The model images have the size of $81$ feet by $83$ feet. 
The length of packers is $2$ feet and they are located at pixel $38$ to $40$. We discretize the well space into grids of 
size $1$ foot by $1$ foot.

After the generation of the bounds for all sources and receivers, we apply the DIRECT algorithm and terminate the iteration when
the difference between the two consecutive objective functional values is less than a certain tolerance.
We set the maximum number of iterations to be $100$, the maximum number of functional evaluations to be $2000$  
and the maximum number of rectangle divisions to be $1000$. The Jones parameter is $10^{-5}$.

The homogeneous background velocity model is provided by the Cranfield CASSM experiment. It is based on the extrapolation 
of the CFU 31F-1 sonic extrapolated along a $-1.9$ degree dipping to the $F_{2}$ and $F_{3}$ well section. All experiments 
are done using the infinite frequency kernel.

It is evident that the kernel $G$ is different in each iteration of the DIRECT algorithm due to the change in the locations 
of sources and receivers. Thus, the goal of the optimal experimental design in this work is to generate the optimal kernel 
$G$ that gives us the best reconstructed image with reasonable financial cost.

\subsection{Numerical Experiment 1: $S/R = 0.5$}

We first show the case in which the ratio of numbers of sources and receivers is fixed. In our case, we set the ratio to be 
$0.5$. 

\begin{figure}[ht!]
 	\centering
 	\includegraphics[scale=0.35]{Ltemp.jpg}
 	\caption{\label{LCO2}L-curve for $CO_{2}$ injection monitoring.}
\end{figure}

We have tested our algorithm on different numbers of sources, with the number ranging from $2$ to $32$. The corner of the 
L-curve happens when we use eight sources and sixteen receivers. It is obvious that this placement has cut out half of the 
financial cost while the result quality is almost as good.

\begin{figure}[ht!]
 	\centering
 	\includegraphics[scale=0.35]{trueWin.jpg}
 	\includegraphics[scale=0.35]{mrWin.jpg}
 	\caption{\label{LCO2}The true image (left) and the reconstructed one within the interest zone based on the optimal 
number of source and receivers (right).}
\end{figure}

We display the true image and the reconstructed one within the interest zone based on the optimal number of source and
receivers. The reconstruction quality is quite satisfying. The relative error is only about $1.5\%$. 

\begin{figure}[ht!]
 	\centering
 	\includegraphics[scale=0.35]{Resultray.jpg}
 	\caption{\label{Resultray}The resulting raypaths based on the optimal sources and receivers.}
\end{figure}

The bounds for placing the four sources above the packer are $[1, 31]$, $[3, 33]$, $[5, 35]$, $[7, 37]$ and $[41, 72]$, 
$[43, 74]$, $[45, 76]$, $[47, 78]$ for the four sources below the packer. For receivers, the bounds are $[1, 21]$, $[15, 35]$ 
for the top group and $[43, 64]$, $[57, 78]$ for the bottom group. In Figure \ref{Resultray} we show the resulting raypaths 
based on the optimal placement of sources and receivers. The position for the four top sources are at $8$, $16$, $23.3$, $35.3$ feet
from the top of the image and $42.7$, $48.2$, $70.8$, $52.2$ feet away for the bottom four sources. For receivers, the highest
one in the top group is located at $4.3$ feet away from the top of the image and the lowest receiver is located at $33.9$ feet
away. For the bottom group, the highest one is at $44.2$ feet away and the lowest one is at $60.5$ feet away. 

\subsection{Numerical Experiment 2: $S/R$ unfixed}

In general, there is no fixed ratio between numbers of sources and receivers; this means that the relative error is a 
function of two independent variables: the number of sources and the number of receivers. In fact, each of these variables 
can be split into top and bottom groups. This will give us a 4D hyper-curve, which is very difficult to determine where
the corner point of the curve is.

In order to still obtain an L-like curve so that the optimal solution is obvious, we consider plotting the relative error
as a function of the financial cost of the sources and receivers. Since sources and receivers have different costs and
these costs also differ based on the placement relative to the packer, we assign values and simply total the cost.
For each total cost, we choose the smallest relative error among the trials we ran.

In the absence of real cost data, we have assigned a cost of 10 monetary units for a source above the packer and a cost of
20 units for a source below the packer. Also we assign costs of 5 and 10 units for receivers above and below the packer,
respectively. The resulting L-curve is shown in Figure \ref{Costcurve}.

\begin{figure}[ht!]
 	\centering
 	\includegraphics[scale=0.35]{costcurve.jpg}
 	\caption{\label{Costcurve}The L-curve shows the relative error as a function of the financial cost of the sources 
and receivers.}
\end{figure}

According to the L-curve, we choose the point corresponding a monetary value of 110 units, which yields a relative error
of $2.21\%$. In our experiment, this corresponds to placing two sources above the packer and three below, also two recivers
both above and below the packer. 

This curve is very useful in practice because the amount of funding available is frequently bounded and thus, one can 
easily find a point on this curve which corresponds to the best reconstruction quality possible within the budget.

Figure \ref{mr2322} displays the true model together with the reconstructed one based on the optimal placement we have chosen
within the interested zone. From the result, we observe that the optimal placement has cut out more than $2/3$ of the financial 
cost while still maintaining a quite good reconstruction quality. 

\begin{figure}[ht!]
 	\centering
 	\includegraphics[scale=0.35]{trueWin.jpg}
 	\includegraphics[scale=0.35]{mr2322.jpg}
 	\caption{\label{mr2322}The true image (left) and the reconstructed one within the interest zone based on the optimal 
number of source and receivers (right).}
\end{figure}