## 2 Materials and Methods
### Reflectivity Spectrum Extrapolation Method
Any missing frequency outside the original bandwidth can be computed using Equation (7) below [39], assuming that \(K\) is the number of impulse pairs in the analysis window (including zero), while \(n\) varies from zero to the number \(N\) to accommodate all bed thicknesses and \(dt\) is the sample rate. The broadband spectrum which can range from  
\[R_{c}=\frac{A_{i2}-A_{i1}}{A_{i2}+A_{i1}} \tag{4}\]1 Any arbitrary pair of reflection coefficients \(r_{1}\) and \(r_{2}\) (Figure 7) can be represented as the sum of even and odd components. The even pair has equal magnitude and sign, and the odd pair has equal magnitude and opposite sign. The even and odd pair can be further decomposed into unit even and odd pairs respectively. zero to the Nyquist frequency (i.e., the maximum frequency that can be observed in the data) of any reflectivity series within the analysis window can be clearly expressed as:  
\[R_{c}=\frac{A_{i2}-A_{i1}}{A_{i2}+A_{i1}} \tag{4}\]2  
In this study, the coefficients (\(r_{e}\)) and (\(r_{o}\)) components were computed using the basis pursuit inversion algorithm proposed by [17,21]. Many methods for spectral decomposition have been previously described and analyzed by [28,46,47,48]. Spectral decomposition is non-unique, and each method has advantages and disadvantages. To use the spectral analysis for layer thickness estimation, precision is needed in both the frequency domain (in order to resolve the thickness) and in the time domain (in order to predict its location in time). Thus, the inversion-based basis pursuit decomposition (BPD) was used for this analysis. This technique proposes that since the sample rate is \(dt\), each trace of even reflectivity consists of a pair of equal impulse functions (spikes) with interval \(n.dt\)[17]. The reflector kernel matrix for the reflectivity pair is constructed by shifting the reflectivity pair along the time axis with \(m.dt\), where \(m\) ranges from 1 to the maximum number of samples in the seismic trace [17]. So, each even (\(r_{e}\)) and odd (\(r_{o}\)) reflectivity are given in time domain as Equations (11) and (12) respectively. \[r_{e}\left( \begin{smallmatrix}t&m,\ n,\ \Delta t\end{smallmatrix} \right)= \delta\left(t-m.dt\right)+\delta(t-m.dt+n.dt) \tag{11}\] \[r_{o}\left( \begin{smallmatrix}t&m,\ n,\ \Delta t\end{smallmatrix} \right)= \delta\left(t-m.dt\right)-\delta(t-m.dt+n.dt) \tag{12}\]  
For our purpose, this same operation was applied in the frequency domain to compute the even and odd impulse pairs. A sum of frequency sinusoids multiplied by the wavelet spectrum is the narrowband data spectrum, which is associated with a sum of reflector pairs. The convolution of reflectivity with seismic wavelets in the time domain is equivalent to multiplication with wavelet spectrum in the frequency domain [49]. Therefore, to obtain reflectivity in the frequency domain, the seismic data spectrum (in the frequency domain) is divided by a known or pre-estimated wavelet spectrum. Since the wavelet is known, the reflectivity spectrum of the bandlimited data is obtained by dividing it by the wavelet spectrum [50]. This produces a flattened spectrum, which corresponds to the wavelet spectral limit section of the broadband reflectivity spectrum obtained using Equation (10), and a best-fit sinusoid is readily fitted to these amplitude points. It is then easy to extrapolate the layer reflectivity spectrum outside the original seismic band by extending the spectra sinusoid (Figure 9). Figure 9: Reflectivity amplitude spectrum for a simple layer of 50 ms thickness, with equal and opposite reflection coefficients of unit amplitude at the top and base. Modified after [19].