\subsection{Static optimization}
In the second step of the algorithm, realizations of each set are iteratively combined using the gradual deformation (GD) method \citep{Hu_2000} as a optimization technique. The GD is a parametrization method that allow to iteratively perturb a set of realizations until obtention of a new realization that minimizes the mismatch with some observed data. This method was first developed for history matching purposes \citep{Roggero_1998} and involves the linear combination of Gaussian random fields with weights that are adjusted to minimize the data mismatch while preserving the spatial covariance. Many variations of this technique have since been presented \citep{Roggero_1998,Hu_2000,Hu_2002,Hu2001,LeRavalec2000,LeRavalec2012}. The key idea behind the GD is that the sum of two Gaussian random field (Y1 and Y2) is also a Gaussian random field:  
\begin{equation}  
Y(r)=Y_1\ cos(r) + Y_2\ sin(r)  
\end{equation}  
where $r$ are the gradual deformation coefficient. It can be shown that $Y(r)$ is a realization with the same mean and covariance as $Y_1$ and $Y_2$ whatever the value of the deformation parameter $r$. This primary gradual combination scheme involves independent realizations only, meaning that the $Y_1$ and $Y_2$ realizations are unconditional. To combine conditional realization, a variant of the GD named gradual conditioning (GC) method was developed by \citet{Ying_2000}, \citet{Hu_2002}. In this case, the  weights have to satisfy a new constraint. Their sum has to be one while their sum to the square is also one. Within this conditions, the realization built from two conditional realizations is also conditional to the measured data and reproduce the overall covariance. \citet{Hu_2002} proposes the weights parametrization as follows:
\begin{equation}
\begin{cases}
    \alpha_1 = \dfrac{1}{3} + \dfrac{2}{3}\ cos(r)\\
    \alpha_1 = \dfrac{1}{3} + \dfrac{2}{3}\ sin\bigg(-\dfrac{\pi}{6}+r\bigg)\\
    \alpha_1 = \dfrac{1}{3} + \dfrac{2}{3}\ sin\bigg(-\dfrac{\pi}{6}-r\bigg)
    \end{cases}
\end{equation}
with parameter $r \in [-\pi,\pi]$. In comparison with the combination of independent realizations, the above procedure avoids the conditioning step during the optimization process.

When the GCM method is incorporated into the optimization processes, the deformation parameter $r$ becomes the parameter to adjuste to reduce the objective function.  At this stage, for each combination (e.g., for each value of $r$) we obtain a new $V_p$, $V_s$ $\rho$ and $\phi$ field model. Afterwards, the full-wave synthetic seismic response of the earth models $d_{synth}$ is computed using an elastic finite-difference time-domain approach \citep{Bohlen_2002} and compared to the observed data $d_{obs}$. For this study, the objective function is defined as the root mean square error between $d_{obs}$ and $d_{synth}$:
\begin{equation}
J(r) = \sqrt{\frac{1}{N}\sum\limits_{i=1,N}\big(d_{synth}(r)-d_{obs}\big)^2}
\end{equation}
where N is the number of the samples per trace and $i$ the sample index.
The combination for which $r$ minimize the mismatch between synthetic and observed data is retained and combined with two initial SGCS realizations. Ultimately, we obtain $V_p$, $V_s$, $\rho$ and $\phi$ fields that have the best seismic match with the observed data. 
    
  