A Novel Methodology for Hydrocarbon Depth Prediction in Seabed Logging: Gaussian Process-Based Inverse Modeling of Electromagnetic Data

: Seabed logging (SBL) is an application of electromagnetic (EM) waves for detecting potential marine hydrocarbon-saturated reservoirs reliant on a source–receiver system. One of the concerns in modeling and inversion of the EM data is associated with the need for realistic representation of complex geo-electrical models. Concurrently, the corresponding algorithms of forward modeling should be robustly efﬁcient with low computational effort for repeated use of the inversion. This work proposes a new inversion methodology which consists of two frameworks, namely Gaussian process (GP), which allows a greater ﬂexibility in modeling a variety of EM responses, and gradient descent (GD) for ﬁnding the best minimizer (i.e., hydrocarbon depth). Computer simulation technology (CST), which uses ﬁnite element (FE), was exploited to generate prior EM responses for the GP to evaluate EM proﬁles at “untried” depths. Then, GD was used to minimize the mean squared error (MSE) where GP acts as its forward model. Acquiring EM responses using mesh-based algorithms is a time-consuming task. Thus, this work compared the time taken by the CST and GP in evaluating the EM proﬁles. For the accuracy and performance, the GP model was compared with EM responses modeled by the FE, and percentage error between the estimate and “untried” computer input was calculated. The results indicate that GP-based inverse modeling can efﬁciently predict the hydrocarbon depth in the SBL.


Introduction
Seabed logging (SBL) is an application of the marine controlled-source electromagnetic (CSEM) surveying technique for discovering and characterizing high-resistive bodies (e.g., hydrocarbon-saturated reservoirs) remotely in the deep marine environment based on contrasts of electrical resistivity of the subsurface. Its capability of identifying the fluid content inside the potential reservoirs is undoubted. This emerging application uses a powerful electric dipole transmitter to continuously emit low-frequency electromagnetic (EM) waves, normally between 0.1 and 10 Hz [1], in all directions through the seawater column and subsurface beneath the seabed. The SBL exploits the fact that hydrocarbonsaturated reservoirs have higher electrical resistivity (approximately 30-500 Ωm) than its surroundings, such as seawater (approximately 0.5-2.0 Ωm) and sedimentary rocks (approximately 1.0-2.0 Ωm) [2]. The conductivities (i.e., inversely proportional to resistivities) of hydrocarbon-filled reservoirs can reach about 0.01-0.1 S/m or lower in certain depths underneath the seabed. The main component of an SBL is the utilization of a towed EM source known as a horizontal electric dipole (HED) and stationary EM receivers for Numerical modeling is a crucial task in geophysics problems since it provides prior knowledge and information of the potential reservoirs. Many powerful computational modeling techniques have been exploited for processing the EM responses, such as finite element (FE), finite difference (FD), finite volume (FV), integral equations (IE), etc. These numerical techniques are capable of modeling the subsurface realistically in the marine environment by defining grids or meshes. According to [11], the FE technique is highly preferred due to its ability to support unstructured meshes. These meshes are able to provide a quality geological representation by placing small grid elements in the areas where better resolution is highly needed, such as in large conductivity contrast areas [11]. However, as the consequence, the mesh-based algorithms require a very lengthy computational time to solve the linear system; in addition, these forward modeling techniques involve very complicated deterministic mathematical equations. According to [12,13], the most time-consuming task in the algorithms is solving the integral and multiple linear equations, not to mention for the inversion scheme that requires multiple EM forward solutions [13]. Here, the corresponding forward modeling should be robust and capable of operating the EM simulation rapidly, especially for repeated use in the large number of iterations involved in the inversion scheme as well as for multiple source-receiver positions in the SBL application [14].
To address these problems, this work aimed to develop a novel Gaussian processbased inversion methodology to provide meaningful information where the hydrocarbon potentially resides underneath the seabed. For the optimization procedure, this inverse modeling utilized the gradient descent (GD) method to find the best estimate (i.e., hydrocarbon depth) based on the forward Gaussian process (GP) regression model. As mentioned before, GP acts as the alternative forward modeling technique as it is capable of being a surrogate model for the complicated mathematical model that requires high computational time to solve [15]. Researchers in [16] also stated that GP is an elegant and simple approach, where it performs well in modeling nonlinear problems. Mean squared error (MSE) was used as the objective function (i.e., loss function) in the optimization. The main focus of this work was to find the depth of hydrocarbon of observational EM responses, whose hydrocarbon layer was assumed to be present but with an unknown depth. This work is novel in two ways. First, it presents the integration of a stochastic process (the GP regression model acts as the EM forward model) and an indirect search algorithm (the iterative GD method as the acquirer of minimizer) in one inversion methodological flow. Second, applying this proposed inverse modeling for predicting the input parameter of the potential high-resistive bodies (e.g., hydrocarbon depth) is a new procedure or approach in SBL data processing. This attempt could be very beneficial to de-risk hydrocarbon exploration in the offshore environment. Details of the GP and GD approaches are discussed in the next section.

Mathematical Background
Gaussian process (GP) and gradient descent (GD) algorithms are commonly exploited in machine learning (ML) problems. In order to seek the most favorable fairness between the performance of the modeling and the need of low computational time, considering the strengths and integrating these two mathematical approaches could be very advantageous to any applications that struggle with the same concern. Here, Figure 2 shows the common ML process to predict outputs given a set of inputs, its corresponding outputs, and desired/query inputs. It represents the general idea on which this work is focusing. The proposed inversion methodology is believed capable of providing meaningful and useful information to the SBL application with an acceptable accuracy and low computational efforts. In this section, the background and examples from the literature of GP and GD approaches are thoroughly discussed in Sections 2.1 and 2.2, respectively.

Gaussian Process
Mathematically, the Gaussian process (GP) can be defined as an infinite collection of random variables where any of the finite subsets can be treated as Gaussian (i.e., normal) distribution [17,18]. A GP model is fully expressed or specified by a mean function m(x) in a vector and a covariance k(x, x ) in the form of a matrix. GP is a probabilistic and nonparametric method for fitting any functional forms based on observations of the domain. This method differs from most of the other black-box identification approaches. GP does not try to approximate the model by trying to fit the parameters of the chosen basis function, but it searches the relationship among the given knowledge or information. This stochastic approach can provide predictive mean and variance. The mean represents the most likely outputs, whereas the variance plays a role as the uncertainty's quantifier. The magnitude of variance is incorporated into the GP prediction and inference. It does reflect the amount and quality of the data distribution, and it is an essential numerical feature when it comes to differentiating the GP regression from any other computational intelligence methods.
GP regression is a well-known technique and has been adopted in many scientific and engineering fields, such as machine learning, computer experiments, geostatistics, electronics, etc. Its convenient properties in modeling make GP successful in demonstrating outstanding performances in numerous applications. Researchers in [19] employed GP regression for predicting the concentration of ozone in an area of Bourgas, Bulgaria. The prediction was based on the measurements of the concentration of ozone, phenol, sulphur dioxide, and benzene in the air (hourly). Researchers in [20] studied the capability of GP in estimating load-bearing pile capacity. This research work compared the performance of GP prediction with other approaches, and GP worked very well in predicting the pile capacity. In a wireless networking application, researchers in [21] inferred that, besides resulting in low-cost outcomes, GP was also capable of obtaining low complex proximity reports for a received-signal-strength (RSS) threshold. Next, GP was also used by researchers in [22] to explain the uncertainty in estimating the state of health (SOH) of lithium-ion batteries. Additionally, according to researchers in [23], GP can successfully decrease the complexity of the computer vision with a better performance of reconstructions.
In [24], GP was hybridized with particle swarm optimization (PSO) to do advanced prediction of roadway broken road zone based on data collected in China. The study inferred that the hybrid model showed the best performance with higher variance accounted for (VAF), higher coefficient of determination (R 2 ), and lower root mean squared error (RMSE) compared to the other methods such as multiple linear regression (MLR) and artificial neural network (ANN). For geophysics problems, researchers in [25] mentioned that GP was able to yield reliable information of permeability and porosity in well analysis. GP was also used by researchers in [26] to predict the presence of a thin hydrocarbon layer by utilizing its predictive variance for quantifying the uncertainties of EM responses. Furthermore, GP was exploited for seismic signal detection as well. According to [27], the rate of detection when using GP was higher compared to other methods and the rate of false alarm was lower especially for the weak seismic signals. Researchers in [28] proved that GP was successfully capable of being exploited in processing EM responses by calibrating the stochastic and the computer experiment models for magnitude versus offset (MVO) analysis.

Gradient Descent
Gradient descent (GD) is a popular iterative process used to update the model parameters of a function by minimizing the objective function (i.e., loss function) based on a first-derivative (i.e., slope) basis. The derivative of a function is fully employed to determine whether to increase or decrease the model parameters (i.e., direction to move). Moreover, finding the slope can also tell how big of a step should be taken to reach the optimal point. Conceptually, this approach optimizes the function by iteratively moving in the negative direction of the gradient, also known as steepest descent. The smoothness of the moving process to reach the minimum depends on the size of steps, which in the GD method is called the learning rate, i.e., α. A higher learning rate can extend over a larger function area, but there is a risk of overshooting the optimal point. Meanwhile, a smaller learning rate is more precise since the slope is frequently re-calculated. However, if it is not wisely parameterized, the optimization needs more time to converge and obtain the minima. Thus, choosing the suitable step size implies a high efficiency of the optimization algorithm.
This well-known method has been employed and exercised in many scientific applications, mostly in optimization problems such as in machine learning [29]. Researchers in [30] exploited GD to establish an operation model of nodes for the classification of heart conditions based on a deep neural network. Furthermore, researchers in [31] adopted the GD method to produce phase-only computer-generated holograms (CGHs). This study inferred that the GD required low time to optimize the phase-only CGH distribution and provided higher-precision images compared to the Gerchberg-Saxton (GS) algorithm. Next, researchers in [32] used the GD method to optimize the vector centers of the consequent layer functions and receptive field matrices in a neuro-fuzzy model based on the standard criterion of mean square error. Furthermore, GD was also used by [33] for minimizing the equation of error used in the study by updating the W 1 and W 2 matrices in order to solve the predicting student performance (PSP) problem. Researchers in [34] also utilized stochastic GD as one of the machine learning methods to be compared for the inversion of a stochastic skin optical model. The researchers also used multi-layer perceptron (MLP), linear support vector regressor (SVR), convolutional neural network (CNN), and lasso regressor in the comparative study.

Methodology
This work aimed to predict the depth of hydrocarbon in the seabed logging (SBL) application by using a novel methodology, namely Gaussian process (GP)-based inverse modeling. Here, four subsections are discussed: 3.1. Preprocessing computer outputs, 3.2. Evaluating electromagnetic (EM) profile using the GP, 3.3. Predicting depth of hydrocarbon, and 3.4. Estimate and GP model validation.

Preprocessing Computer Outputs
Computer simulation technology (CST) software, which uses finite element (FE), is exploited in this work to model and evaluate prior information of EM responses (e.g., magnitudes of electric (E)-field) in the SBL for the proposed GP-based inversion methodology. CST is part of SIMULIA, a Dassault Systèmes brand, sourced from Darmstadt, Germany. According to [35], CST Studio Suite ® exclusively uses the FE method as the EM simulation solver for low-frequency application. Nine SBL models with the same input properties, except the hydrocarbon depth, were simulated using the CST software. Details of the SBL data acquisition are presented in Appendix A. For demonstration, one SBL model replicated by the CST is depicted in Figure A1. The SBL models, known as target models, were replicated with three background layers such as air, seawater, and sediment. The nine different "tried" depths were between 200 m and 1000 m from the seafloor with an increment of 100 m each. All the computer input properties and parameters used in the SBL modeling are indicated in Table 1. As discussed by researchers in [28], interpretation of the SBL data relies on an assumption that the normalized magnitude versus offset (MVO) between the target responses (i.e., with presence of hydrocarbon) and reference responses (i.e., without presence of hydrocarbon) can pre-identify the presence of high-resistive bodies. The idea is that the magnitudes of target response are divided by the magnitudes of reference response along the offsets (i.e., source-receiver separation distances). Ratio of the normalization should be greater than one, which is associated with the existence of the potential reservoirs. Thus, all the generated CST outputs were preprocessed before being used for the forward GP modeling. By considering the ratio measurements and the symmetrical property of the developed SBL models, we chose nine different datasets with 73 datapoints per set, which are from an approximate distance of 11,840.80 m to an approximate distance of 19,004.97 m, for the data processing presented in the next subsection.

Evaluating Electromagnetic Profile Using the Gaussian Process
This work uses Gaussian process machine learning (GPML) MATLAB codes for the GP regression, previously written by Rasmussen and Nickisch, which can be found in [36]. Note that the MATLAB version used by this work is R2019b. As mentioned before, the Gaussian process (GP) acts as the alternative approach to the numerical forward modeling in the SBL application for evaluating EM profiles at various "untried" hydrocarbon depths. Suppose all nine CST datasets, D 1 , D 2 , . . . , where N = 73, are treated as the training datasets for the GP modeling. It is also possible to simply write the d-dimensional input variable vector and its corresponding output, which is magnitude of E-field vector, as x i ∈ X and y i ∈ Y, respectively. Generally, in this work, X is the N-by-d matrix and Y is the N-by-1 column vector, where x ∈ R d=2 and y ∈ R d=1 . In other words, x = (s, h), where s and h represent the offset (i.e., source-receiver separation distance) and the depth of hydrocarbon, respectively.
In order to evaluate or estimate the EM profiles at various "untried" depth inputs conforming to the testing input, x * = (s * , h * ), the posterior distribution of the testing output y * is defined in Equation (1). This work's interest was to predict the depth parameter of hydrocarbon along the offsets from 11,840.80 m to 19,004.97 m. Thus, here, only the testing input variable of hydrocarbon depth, h * , is varied accordingly, while the s * remains unchanged. Note that the interested hydrocarbon depths in the GP regression are from 200 m to 1000 m with an increment of 20 m each (including the "tried" depths).
where the predictive mean, m * , and the predictive variance, κ * , are expressed in Equations (2) and (3), respectively. and Based on Equations (2) and (3), K(X, X) is the covariance matrix with a dimension of Nby-N amongst the X, K(X, x * ) is the covariance matrix with a dimension of N-by-1 between the X and x * , k(x * , x * ) is the covariance of the x * , I N is an N-dimensional unit matrix, and σ 2 N I N is the noise covariance matrix. All these covariance matrices are composed of squared exponential (SE) terms, where each element of the matrices is computed using Equation (4) as follows: where σ f represents the signal variance and represents the isotropic length scale.These two parameters are known as the hyperparameters of the SE covariance function, θ = σ f , . SE is a very famous covariance function used in GP regression. According to [37], GP with SE is infinitely differentiable. This feature is very important for the derivativebased approach.
Here, it must be noted that all the computations of Equations (1)-(4) require the estimate of θ. These hyperparameters are numerically estimated by minimizing the negative log marginal likelihood (NLML) using a gradient-based optimizer. The NLML equation is defined in Equation (5) as follows: In Equation (5), a logarithmic identifier is used to simplify the integral expression involved in the marginal likelihood method [25]. Once the hyperparameters are successfully estimated, the estimation of the GP regression is then defined using Equations (2) and (3). The m * is a mean value of the estimation that gives the best estimate of y * (i.e., magnitude of E-field at x * ), and the κ * is the predictive variance of the testing input in terms of two standard deviations (95% confidence interval), which reflects the reliability of the predictive mean. These two equations are the basis in the forward GP modeling.

Predicting Depth of Hydrocarbon
The fitted GP model is denoted as Y f it (s, h), where h here includes all "tried" and "untried" depths. This means that, for every h (i.e., 200 m-1000 m with a gap of 20 m), the GP gives estimates of the EM profile, Y f it (s, h), at distances from 11,840.80 m to 19,004.97 m. Assume that the observational dataset Y obs (s, ) is observed at N = 73 distance points and the hydrocarbon resides underneath the seabed but with "unknown" depth. To find out the depth of hydrocarbon of the observational dataset Y obs , the least square criterion is formed between the Y f it and Y obs . Mean squared error (MSE) is used as the loss function in the optimization algorithm. In general, MSE measures the average amount of squared differences between Y f it and Y obs . The loss function L h j is expressed in Equation (6) as follows: where L h j j = 1, 2, . . . , M and M = 41. Equation (6) is a function of h. The loss function is minimized over the h with the minimizer serving as an estimate of the hydrocarbon depth. By using Equation (6), the L h j between the Y f it and Y obs are calculated. Here, the least L h j is then identified. As for now, the minimizer of the least L h j serves as the best estimate of hydrocarbon depth for Y obs . Thus, in this work, the respective minimizer is treated as the initial guess of depth. Then, the iteration index j is set as zero and the initial guess is denoted as h 0 .
Next, for the optimization using the GD method, the first-partial derivative of the L h j with respect to h j is required. The partial derivative of the loss function is defined in the equation below. The derivation of Equation (7) is provided in Appendix B (in Equation (A1)): where Y f it is the derivative of the GP using the SE covariance function with respect to h j . It must be noted here that the forward GP model is defined by estimations composed by the predictive mean as shown in Equation (2). To find out the Y f it , the predictive mean m * is differentiated with respect to h j as follows: where K (X, x * ) is the derivative of covariance function between X and x * with respect to h j . As mentioned earlier, each element of the covariance matrix is computed using the SE covariance function (Equation (4)). Thus, the derivative of SE between X and x * with respect to h j is defined in order to find the K (X, x * ), as expressed in Equation (9). Its derivation is provided in Appendix B (in Equation (A2)): Next, the following step is finding the correction or update of h j . It is denoted as h j+1 . By referring to the GD algorithm, the equation of h j+1 is expressed as follows: where α is the learning rate, which governs the step size of the optimization. Normally, as a start, a geometric sequence with a common ratio of 10 is used to indicate which minimizer can give a significant optimization. After getting the h j+1 , the EM profile at h j+1 , denoted as Y f it s i , h j+1 , is then estimated using the GP regression by defining the h j+1 as its testing input and the GP modeling steps are repeated. Next, the new loss function for h j+1 between Y f it and Y obs is calculated using Equation (11) as follows: The tolerance threshold ε, which acts as the stopping criterion for the iterations, is defined for verifying the L h j+1 before deciding the optimal hydrocarbon depth for the Y obs . In the context of this work, it was decided that the changes between each subsequent minimizer or depth had to be less than 0.05. As presented in Section 4, the specific limit leads to acceptable results. Equation (12) shows the equation of relative percentage change (RPC) used in this work: The next conditions used in this work are the following: the maximum number of iterations is 100 and the loss function (i.e., MSE) must converge. In other words, L h j+1 must be always lower than L h j . Each iteration must satisfy these three threshold conditions. The optimal depth of hydrocarbon for Y obs is believed to be successfully predicted if the process does not meet the conditions. Otherwise, it is returned to the step of evaluating the L h j onwards by defining a new iteration index, j = j + 1. All the flows are repeated until the ε is neglected, then the iterations stop. The methodological flow presented in this section is shown in Figure 3.

Estimate and GP Model Validation
Error measurement is very important to determine the reliability of estimations. In this work, for estimate validation, the percentage error between the computer input (i.e., hydrocarbon depth used in the computer simulation) and the depth predicted by the GP-based inverse modeling was calculated using Equation (13). Meanwhile, for GP model validation, three error validators were used in this work, namely mean absolute deviation (MAD), mean absolute percentage error (MAPE), and root mean squared error (RMSE). The errors of EM profiles provided by the GP model were calculated by assuming the EM responses modeled by the FE (through the CST software) as the true EM values. All equations of the MAD, MAPE, and RMSE are defined in Equations (14)- (16), respectively, as follows: where h true is the CST depth input, h estimate is the depth estimated by the GP-based inverse modeling, y FE is the EM response modeled by FE method, y GP is the EM response evaluated by the forward GP modeling, and N is the total number of observations.

Results and Discussion
All prior electromagnetic (EM) responses were generated using the computer simulation technology (CST) software. Nine seabed logging (SBL) models were parameterized with the same input properties, except for its depth of hydrocarbon, which were 200 m-1000 m with a gap of 100 m each. It must be noted that these EM responses were evaluated using the finite element (FE) method, which was used exclusively by the CST software for the low-frequency applications. In this work, only magnitude of electric (E)-field was considered as the corresponding output. Figure 4 shows the magnitude versus offset plot for all nine SBL datasets. Based on Figure 4, the x-axis represents the offset distances from approximately 11,840.80 m to 19,004.97 m and the y-axis represents the magnitude of E-field with a logarithmic scale. This figure reveals that the CST was capable of simulating the SBL application in the deep marine environment. The generated EM responses behave similarly to the real behavior of EM signals in media with various conductivities. The EM signals travel through a high-resistive body with less attenuation. Thus, a model with higher hydrocarbon depth generates a much smaller magnitude of E-field due to the factor of far EM signal propagation. Next, these EM responses were fully utilized by the GP regression in order to provide more EM profiles at "untried" depths. The contour plot of the twodimensional (2D) GP model is depicted in Figure 5. Based on Figure 5, the x-axis is the offset distance and the y-axis is the 41 depths of hydrocarbon from 200 m to 1000 m, while the magnitudes of the E-field are represented in color codes. This forward GP model is a combination of EM profiles at "tried" and "untried" hydrocarbon depths. The purpose of this GP modeling was to evaluate as many EM responses as possible at various depths of hydrocarbon for the inversion scheme, for which generating all these EM responses using the CST software could take high computational time. Simulating and generating all these 41 datasets using the CST could require around 3 h and 30 min, whereas the GP modeling only took 45 min for the prior knowledge acquisition and less than 1 min for developing the forward GP model and evaluating 41 profiles.
In this work, three random observational datasets were generated separately using the CST software where the depths of hydrocarbon were assumed to be "unknown", i.e., 350, 650, and 950 m. Here, the goal was to find the best estimates of hydrocarbon depths of the observational datasets with the smallest MSE values. In other words, a small MSE value indicates that the EM profile evaluated by the forward GP model gives the most likely values for the observational dataset, thereby giving the best depth information. As discussed before, the step size of the optimization was determined by the learning rate. For this work, three different learning rates were used for every inversion scheme. The first scheme (1st observational dataset) used 0.051 as its learning rate, whereas the second (2nd observational dataset) and third (3rd observational dataset) schemes used 0.087 and 0.081, respectively. Details of how this work determined the suitable learning rates are provided in Appendix C. Plots of MSE versus learning rate and iteration number versus learning rate for every observational dataset are depicted in Figures 2-6 and A1. After identifying the learning rate, optimization using the gradient descent was executed. Thus, three MSE versus hydrocarbon depth plots are plotted as shown in Figures 6-8.
By referring to Figures 6-8, the x-axis is the depth of hydrocarbon and the y-axis is the misfit, which is the MSE. The blue points are the 41 MSE values between the GP model from depths 200 m to 1000 m and the observational data. The red points indicate the optimized MSEs with respect to the updated minimizers. Every optimization must satisfy the convergence condition. Thus, in order to monitor the convergence status, this work illustrates MSE versus number of iterations plots for every observational dataset, which are depicted in Figures 9-11. In these figures, the x-axis and y-axis represent the number of iterations and MSE, respectively. Based on these figures, even though the first inversion scheme needed a higher number of iterations than the second and third ones, all the iterations were smaller than the defined maximum number, which was 100. Moreover, all the MSEs converged where the optimizations stopped at the lowest MSE values. To make it easy to interpret, the main results are presented in Table 2.       were also very small. This means that the deviations between its true EM values and the predicted EM profiles were very low. All the final RPCs were also lesser than 0.05 as compared to what was defined in the methodology. In addition, in the context of this work, the reliability of the GP regression in evaluating the EM profiles was investigated as well by considering data modeled by the FE method as the true EM values. Thus, MADs, MAPEs, and RMSEs between the EM responses modeled by FE and EM evaluated by the GP (at depths 200-1000 m with a gap of 100 m each) were calculated. The results are presented in Table 3. Furthermore, to find out the accuracy and performance of the GP-based inverse modeling in predicting the hydrocarbon depth, percentage errors between the true depths (i.e., computer input) and the estimates were calculated as well. The percentage errors are indicated in Table 4.   Table 3, all the MADs, MAPEs, and RMSEs for depths 200 m to 1000 m with 100 m gap each were very small, which were 1.49 × 10 −9 (on average), 0.114% (on average), and 3.49 × 10 −9 (on average), respectively. As seen in Table 4, the percentage errors were very low as well, which was 0.0260% on average. All the results revealed that GP regression was capable of evaluating the EM profiles at various depths of hydrocarbon accurately and with low time consumption. The main goal was also achieved in that the GP-based inversion methodology can predict the depth of hydrocarbon in the SBL application efficiently.

Conclusions
This work proposed Gaussian process (GP)-based inverse modeling to predict the depth of hydrocarbon in the seabed logging (SBL) application. The capability of computational modeling and inversion of the SBL data with good accuracy and low computational effort is one of the main concerns in electromagnetic (EM) problems. Thus, this work used machine learning-based approaches, namely GP, as the forward modeling of EM profiles and the gradient descent (GD) method as the optimizer for finding the input parameter of the hydrocarbon (e.g., depth) in the SBL data processing. This work utilized three datasets generated through the computer simulation technology (CST) software with random "unknown" hydrocarbon depths as the observational responses. Based on the results, GP regression can provide accurate EM profiles at various depths of hydrocarbon with low computational time (45 min for the prior data acquisition) and less than 2 min for the GP-based inversion methodology. We believed that if the forward modeling performed by the GP could be done rapidly for the repeated use, the inverse modeling (by the GD) would not significantly affect the total computational time of the inversion methodology. Next, the resulting mean absolute deviations (MADs), mean absolute percentage errors (MAPEs), and root mean squared errors (RMSEs) between the GP model and the data modeled by finite element (FE) method were very small. Furthermore, the GP-based inverse modeling could also predict the depth of hydrocarbon for each observation dataset accurately. All the calculated percentage errors between the estimates (by the GP-based inversion methodology) and the true depth values (computer input) were very small as well. Therefore, it can be concluded that the GP-based inverse modeling was successfully developed to predict the depth of hydrocarbon in the SBL. This attempt could help marine EM data processing, especially in order to de-risk hydrocarbon exploration in the deep marine environment.

Acknowledgments:
We would like to express our deepest appreciation to Universiti Teknologi PETRONAS (UTP) for the tremendous support in completing this work.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Data Acquisition Using CST Software
Typical synthetic seabed logging (SBL) models were replicated using the computer simulation technology (CST) software. The input parameters can be seen in Table 1. The first layer was air, followed by seawater and sediment. Based on Figure A1, the isotropic hydrocarbon layer was positioned in the sediment layer at 200 m from the seafloor (for demonstration purposes, only one model is shown here). The source was placed at the center of the model and elevated at 30 m from the seafloor (coordinate: x = 10,000 m, y = 10,000 m, z = 1270 m), and inline receivers were located on the seafloor (i.e., static source-receiver separation distances). All the generated electromagnetic (EM) responses were symmetrical from the center to the left and right boundaries. Equation (A1) shows the derivation of the first-partial derivative of the mean squared error (MSE) with respect to h j .
Equation (A2) shows the derivation of the derivative of the squared exponential (SE) covariance function with respect to the testing point, x * .

Appendix C. Identifying Learning Rate (α) of Gradient Descent
Geometric sequences with a common ratio of 10 were tested first in order to indicate which learning rates would give the significant optimization. From here, other values of learning rate were examined as well to find the smallest mean squared error (MSE) that could be optimized. Based on Figure 2, Figure 4, and Figure 6 3.30 × 10 −7 , 3.90 × 10 −7 , and 5.10 × 10 −7 were the smallest optimized MSEs for all observational data. Thus, based on the MSEs, the size of the learning rates was narrowed down. Next, the number of iterations was then considered as plotted in Figure 3, Figure 5, and Figure 7. From here, learning rates that gave the smallest MSE with the lowest iteration number and followed the relative percentage change (RPC) condition were chosen for the Gaussian process (GP)-based inverse modeling.