\append{Gradient of ISWET by the adjoint-state method}
The gradient of the objective function can be computed with the Frech$\acute {\rm e}$t derivatives. However, even for 2D applications of ISWET this computation can be very expensive. An efficient way to compute the gradient without using Frech$\acute {\rm e}$t derivatives is the adjoint-state method \cite[]{chavent:476,plessix2006}. \cite{plessix2006} describes two methodologies for computing the gradient of the objective function using the adjoint-state method. One methodology uses the perturbation theory, which states that, at first order, a perturbation of the model parameters causes a perturbation of the objective function. The other uses the augmented Lagrangian. The augmented Lagrangian is formed by the objective function and the scalar product of the adjoint-state variables with general solutions of the forward modeling equations. The adjoint-state variables are, in turn,  solutions of the adjoint-state equations. The adjoint-state equations are defined by equating to zero the derivatives of the augmented Lagrangian with respect to the state variables. For the linear case, the adjoint of the modeling operator applied to the adjoint-state variables gives the gradient of the objective function.
\par 
Although previous studies have computed the gradient of the ISWET objective function by the adjoint-state method (for example, \cite{shen:2132}), they have not provided a detailed derivation. Here, I show a detailed derivation the gradient of the ISWET objective function using the augmented Lagrangian methodology. The derivation is valid whether ISWET uses areal-shot migration or shot-profile migration. A complete description of the forward and adjoint of ISWET operators is given in \cite{yax2008}. A numerical example of slowness optimization on the Marmousi model illustrates the use of the gradient by a quasi-Newton algorithm.
\par



Image-space wave-equation tomography aims to iteratively solve for the slowness model, $s=s({\bf x})$, that minimizes the nonlinear objective function:
\begin{eqnarray}
J(s)=\frac{1}{2} \left | \left | \Delta r(s) \right | \right | ^2 = \frac{1}{2} \left | \left | r(s) - {\bf M} r(s) \right | \right | ^2,
\label{eq:objfun}
\end{eqnarray}
where $\Delta r = \Delta r({\bf x},{\bf h})$ is the image perturbation that measures the accuracy of the slowness model. To compute $\Delta r$, a differential residual-focusing operator $\bf M$ is applied to the image $r = r({\bf x},{\bf h})$ obtained with the current slowness \cite[]{biondi2008}, using either differential residual prestack migration \cite[]{sava2004a,sava2004b} or differential-semblance optimization (DSO) operators \cite[]{shen:VE49}. In this paper, operators are represented by bold capital letters.
\par
If the differential residual-focusing operator $\bf M$ is independent of the slowness, the gradient of this objective function evaluated at the current slowness ${\widehat s}={\widehat s}(\bf x)$ is
\begin{eqnarray}
\nabla J(s)= \left. \left ( \frac{\partial r}{\partial s} \right )'\right |_{s=\widehat s} \left ({\bf I - M}' \right )  {\Delta \widehat r},
\label{eq:gradobjfun}
\end{eqnarray}
where the prime denotes the adjoint, $\bf I$ is the identity operator, and $\Delta \widehat r=\Delta \widehat r({\bf x},{\bf h})$ is the perturbed image obtained with the current slowness model.  The linear operator $\frac{\partial r}{\partial s}$ defines the mapping $\frac{\partial r}{\partial s} \Delta s = \Delta r$ between the slowness perturbation $\Delta s$ and the image perturbation $\Delta r$; it is called the image-space wave-equation tomographic operator.
\par
Because the image-space wave-equation tomographic operator is composed of different operators, it is difficult to envision from equation \ref{eq:gradobjfun} which operations are performed to compute the gradient. Therefore, for a clear explanation of the operators involved, I use the adjoint-state method to derive the gradient of the objective function (equation \ref{eq:objfun}).
\par
In migration with generalized sources or shot-profile migration, the source and receiver wavefields are propagated independently, and the image $r_z = r_z({\bf x},{\bf h})$ at a depth level $z$, is computed by the crosscorrelation
\begin{eqnarray}
r_z({\bf x},{\bf h}) = \sum_{\omega} d_z^*({\bf x}-{\bf h},\omega) u_z({\bf x}+{\bf h},\omega), 
\label{eq:imgcond1}
\end{eqnarray}
where $d_z({\bf x},\omega)$ is the source wavefield for a single frequency $\omega$ at horizontal coordinates ${\bf x}=(x,y)$ ; $u_z({\bf x},\omega)$ is the receiver wavefield, ${\bf h}=(h_x,h_y)$ is the subsurface half-offset, and ``$^*$'' stands for the complex-conjugate. An additional summation over shots is required when migrating more than one shot. Hereafter, letters $d$ and $u$ stand for source and receiver wavefields, respectively, irrespective of the migration scheme.
\par
In a more compact notation, not explicitly writing the dependencies on $\bf x$ and $\bf h$, equation \ref{eq:imgcond1} can be re-written as follows:
\begin{eqnarray}
r_z = {\bf{S}}{\bf{D}}_z'(\omega) u_z(\omega) = {\bf{S}}{\bf{U}}_z(\omega) d_z^*(\omega), 
\label{eq:imgconcompact}
\end{eqnarray}
where ${\bf{D}}$ and ${\bf{U}}$ are convolutional operators composed of $(h_x,h_y)$-shifted versions of $d_z({\bf x},\omega)$ and $u_z({\bf x},\omega)$, respectively. Operator ${\bf{S}}$ corresponds to the summation over frequency.
\par
For subsequent depth levels, $d({\bf x},\omega)$ is computed by means of the recursive downward propagation 
\begin{eqnarray}
\left\{ \begin{array}{l}
d_{z+1}(\omega)  = {\bf T}_z^{\downarrow}(\omega,s)d_z(\omega) \\ 
d_1(\omega) = q(\omega),
\label{eq:downmod}
\end{array}
\right.
\end{eqnarray}
where ${\bf T}_z^{\downarrow}$ is the downward continuation operator, which is a function of the slowness $s$, and $q(\omega)$ is the source wavefield used as a boundary condition. In the case of conventional shot-profile migration, $q(\omega) = f_s(\omega)\delta({\bf x - x_s})$ is the source signature located at ${\bf x}_s=(x_s,y_s,0)$. If using the generalized sources in the image space, $q(\omega)$ represents  the image-space phase-encoded source wavefield of equation~\ref{eq_perm_bD2}.
\par
The downward continuation of the receiver wavefield is performed by
\begin{eqnarray}
\left\{ \begin{array}{l}
u_{z+1}(\omega)  = {\bf T}_z^{\downarrow}(\omega,s)u_z(\omega) \\ 
u_1(\omega) = w(\omega),
\label{eq:upmod01}
\end{array}
\right.
\end{eqnarray}
where $w(\omega)$ is the recorded data at the surface for shot-profile migration. If using generalized sources in the image space, $w(\omega)$ is the phase-encoded areal receiver wavefield of equation~\ref{eq_perm_bU2}. In equations \ref{eq:downmod} and \ref{eq:upmod01}, I omitted the dependencies of the wavefield with respect to ${\bf x}$. The subscript ${\it 1}$ in equations \ref{eq:downmod} and \ref{eq:upmod01} represents the surface for the shot-profile migration and the ``collection'' depth level, $z_c$, for the image-space phase-encoded wavefields.
\par
In the image-space wave-equation tomography problem, the perturbed source and receiver wavefields, and the image perturbation are used to compute the slowness perturbation that updates the current slowness model. From the perturbation theory, we have $d = \widehat d + \Delta d$, $u = \widehat u + \Delta u$, and, consequently, $r = \widehat r + \Delta r$ are physical realizations with $s = \widehat s + \Delta s$, where the $hat$ refers to fields obtained with the background slowness. To the first order (Born approximation), these perturbed fields are given by 
\begin{eqnarray}
\Delta d_{z+1}(\omega) = {\bf T}_z^{\downarrow}(\omega,\widehat s) \Delta d_z(\omega) + \widetilde {\widetilde{\bf D}}_z(\omega) \Delta s_z %\Delta {\bf T}_z^{\downarrow}(\omega,\widehat  s) {\widehat d}_z(\omega),
\label{eq:pertsou}
\end{eqnarray}
and
\begin{eqnarray}
\Delta u_{z+1}(\omega) = {\bf T}_z^{\downarrow}(\omega,\widehat s) \Delta u_z(\omega) + \widetilde {\widetilde{\bf U}}_z(\omega) \Delta s_z. %\Delta {\bf T}_z^{\downarrow}(\omega,\widehat s) {\widehat u}_z(\omega)
\label{eq:pertrec}
\end{eqnarray}
The diagonal operators $\widetilde {\widetilde{\bf D}}_z$ and $\widetilde {\widetilde{\bf U}}_z$ have in the diagonal entries the scattered source and receiver wavefields, respectively. These wavefields are given by the action of the scattering operator $\Delta {\bf T}_z^{\downarrow}$ on the background wavefields:
\begin{eqnarray}
\widetilde {\widetilde {\bf D}}_z(\omega) = \Delta {\bf T}_z^{\downarrow}(\omega,\widehat s) {\widehat d}_z(\omega) = i\frac{\omega^2 \widehat s}{\sqrt{\omega^2 \widehat s^2 - \left | {\bf k} \right | ^2}} dz~ \widehat d_z(\omega)
\label{eq:pscatt}
\end{eqnarray}
and
\begin{eqnarray}
\widetilde {\widetilde {\bf U}}_z(\omega) = \Delta {\bf T}_z^{\downarrow}(\omega,\widehat s) {\widehat u}_z(\omega) = -i\frac{\omega^2 \widehat s}{\sqrt{\omega^2 \widehat s^2 - \left | {\bf k} \right | ^2}} dz~ \widehat u_z(\omega).
\label{eq:uscatt}
\end{eqnarray}
\par
The perturbed image is given by
\begin{eqnarray}
\Delta r_{z} = {\bf S} \left ({\widehat {\bf U}}_z(\omega) \Delta d_z^*(\omega) + {\widehat {\bf D}}_z'(\omega) \Delta u_z(\omega)  \right ).
\label{eq:pertimg}
\end{eqnarray}
The matrix representations of equations \ref{eq:pertsou}, \ref{eq:pertrec}, \ref{eq:pertimg} are
\begin{eqnarray}
\Delta {\underline {\bf d}} = {\bf T}^{\downarrow} \Delta {\underline {\bf d}} + {\widetilde {\widetilde {\bf P}}} {\bf S}' \Delta {\bf s},
\label{eq:pertsou2}
\end{eqnarray}
\begin{eqnarray}
\Delta {\underline {\bf u}} = {\bf T}^{\downarrow} \Delta {\underline {\bf u}} + {\widetilde {\widetilde{\bf U}}} {\bf S}' \Delta {\bf s},
\label{eq:pertrec2}
\end{eqnarray}
and
\begin{eqnarray}
\Delta {\bf {\underline r}} = {\bf S} \left ({\widehat {\bf U}} \Delta {\underline {\bf d}}^* + {\widehat {\bf D}}' \Delta {\underline {\bf u}}  \right ),
\label{eq:pertimg2}
\end{eqnarray}
where ${\bf S}'$ is a spreading operator that replicates the slowness perturbation for every frequency.
\par
Equations \ref{eq:pertsou2}, \ref{eq:pertrec2}, and \ref{eq:pertimg2} are the forward modeling equations of the image-space wave-equation tomography problem using the generalized sources or shot-profile  schemes. They depend on the state variables $\Delta {\underline {\bf d}}$, $\Delta {\underline {\bf u}}$, and $\Delta {\underline {\bf r}}$. \cite{plessix2006} describes how to compute the adjoint states using the augmented functional methodology. By introducing the adjoint-state variables ${\underline \lambda}_d$, ${\underline \lambda}_u$, and ${\underline \lambda}_r$, the augmented Lagrangian reads
\begin{eqnarray}
\begin{array}{l}
\mathcal{L}(\Delta {\underline {\bf d}}, \Delta {\underline {\bf u}}, \Delta {\underline {\bf r}}, {\underline \lambda}_d, {\underline \lambda}_u, {\underline \lambda}_r; \Delta {\bf s}) =  \mathcal{R} \left [ \frac{1}{2} \left | \left | \Delta {\underline {\bf r}} \right | \right | ^2 - \right. \nonumber \\
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ \left \langle {\underline \lambda}_d,\left ( {\bf I} - {\bf T}^\downarrow \right )    \Delta {\underline {\bf d}} - {\widetilde {\widetilde{\bf D}}} {\bf S}' \Delta {\bf s} \right \rangle - \nonumber \\
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ \left \langle {\underline \lambda}_u,\left ( {\bf I} - {\bf T}^{\downarrow} \right )    \Delta {\underline {\bf u}} - {\widetilde {\widetilde{\bf U}}} {\bf S}' \Delta {\bf s} \right \rangle - \nonumber \\
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ \left \langle {\underline \lambda}_r,\Delta {\underline {\bf r}} - {\bf S} \left ({\widehat {\bf U}} \Delta {\underline {\bf d}}^* + {\widehat {\bf D}}' \Delta {\underline {\bf u}} \right )  \right  \rangle. \\
\end{array}
\label{eq:augfct01}
\end{eqnarray}
The adjoint-state variables are computed by taking the derivative of $\mathcal{L}$ with respect to the state variables and equating to zero, which gives
\begin{subequations} 
\label{subeqn:list} 
\begin{align} 
\left ( {\bf I} - {\bf T}^\downarrow \right )' {\underline \lambda}_d = {\widehat {\bf U}} {\underline \lambda}_r, \label{subeqn:adjstvar1} \\
\left ( {\bf I} - {\bf T}^\downarrow \right )' {\underline \lambda}_u = {\widehat {\bf D}} {\underline \lambda}_r, \label{subeqn:adjstvar2} \\
{\underline \lambda}_r = \Delta {\bf {\underline r}}.\label{subeqn:adjstvar3} 
\end{align} 
\end{subequations} 
Notice that 
\begin{eqnarray}
\left ( {\bf I} - {\bf T}^\downarrow \right )' = \left ( {\bf I} - {\bf T}^{\downarrow '} \right ) = \left ( {\bf I} - {\bf T}^{\uparrow} \right )
\end{eqnarray}
corresponds to the recursive upward propagation operator. Therefore, equations \ref{subeqn:adjstvar1} and \ref{subeqn:adjstvar2} can be written as
\begin{subequations} 
\label{subeqn:list} 
\begin{align} 
{\underline \lambda}_p = {\bf T}^\uparrow {\underline \lambda}_d + {\widehat {\bf U}} {\underline \lambda}_r, \label{subeqn:adjstvar10} \\
{\underline \lambda}_u = {\bf T}^\uparrow {\underline \lambda}_u + {\widehat {\bf D}} {\underline \lambda}_r, \label{subeqn:adjstvar20} 
\end{align} 
\end{subequations} 
which correspond to the recursive upward propagation of the perturbed wavefields resulting from the convolution of the wavefields computed with the current slowness and the perturbed image.
\par
Finally, the gradient of $J$ is 
\begin{eqnarray}
\nabla_s J({\bf s}) = {\bf S} \left ({\widetilde {\widetilde{\bf D}}}'{\underline \lambda}_d  +  {\widetilde {\widetilde{\bf U}}}'{\underline \lambda}_u \right ).
\label{eq:adjstgrad}
\end{eqnarray}
To compute the gradient, the adjoint-state wavefields, ${\underline \lambda}_d$ and ${\underline \lambda}_u$, are upward propagated and cross-correlated in time with the scattered wavefields.

