RETRACTED: Imaging Enhancement of Gas Reservoir via Full-Waveform Inversion Imaging: Case Study, Gas Field Open Access

Successful exploration and production operations in the oil and gas industry fundamentally rely on robust seismic imaging techniques. Clear seismic imaging results enable precise analyses at various evaluation stages, including assessing the distribution of hydrocarbons, reservoir characterization, shallow hazard assessment for drilling, and volumetric assessment of oil and gas reserves. In general, effective seismic imaging relies on several key requirements, including the quality of data acquisition, the signal-to-noise ratio (S/N ratio), and robust processing procedures. Under various conditions, the ultimate goal is to maximize the preservation of both low and high frequencies within a reasonable S/N ratio. This broadens the bandwidth, allowing for the building of accurate velocity models essential for interpreting subsurface structures and properties [1].

However, there have been significant limitations in fully utilizing the acquired bandwidth within the conventional imaging process. Over many decades, various methods have been developed to obtain high-resolution images. However, these conventional methods have largely been confined to three main stages: preprocessing, velocity model building, and migration. Each of these stages progresses sequentially, with parameter testing and result comparisons to derive optimal outcomes. Typically, subjective evaluations by processors are involved in this process, and extensive testing procedures result in time consumption [2].

As the initial step in the conventional imaging process, the ultimate goal of preprocessing is to retain only the single scattered (i.e. primary only) wavefield from the entire raw wavefield. Since conventional migration operations consider only a single scattered wavefield from the raw wavefield, other wavefields, such as multiscattered wavefields (multiples, ghosts, transmitted waves, etc.), must be eliminated [3-5]. However, the preprocessing procedure faces the challenge of removing only unwanted noise while preserving primary reflections. Seismic processors must consider a trade-off between adequate noise removal and preservation of primary reflections when determining the final preprocessing parameter settings.

Migration is the process of relocating reflections to their accurate subsurface positions, thereby creating subsurface images [6]. A precise velocity model is fundamentally required to accurately relocate wave energy. However, in complex geological settings and with poor data quality, there are limitations in accurately implementing velocity model building through conventional velocity analysis. Furthermore, even if an accurate velocity model is established, conventional migration still utilizes only primary reflections. Thus, it fails to leverage significant information such as frequency, wave path, and illumination carried by other wavefields [7].

To overcome the limitations of conventional velocity model building and migration techniques, Tarantola [8] and Lailly [9] devised the full-waveform inversion (FWI) method, which utilizes the full wavefield in a data-fitting procedure. The term “full wavefield” encompasses primary reflections, diving waves, supercritical reflections, and multi-scattered waves such as multiples [10]. FWI involves an optimization scheme that iteratively minimizes the misfit between observed data and predicted data via the least-squares norm. Typically, FWI primarily utilizes a single parameter (P-wave velocity) for data fitting; thus, it is mainly employed for velocity model building. Routinely, the results of FWI serve as high-resolution velocity inputs for conventional migration. One of the well-known limitations of FWI is the occurrence of cycle skipping artifacts. If the travel time of the starting model deviates from the observed travel time by more than half of the period, the optimization problem falls into a local minimum (Figure 1) [10]. To prevent local minima, it is essential to determine relatively accurate initial velocities and low-frequency components that are not sensitive to cycle skipping. Recent advancements in computing power, diverse acquisition designs, and developments in FWI algorithms have contributed to solving these issues. Additionally, the construction of high-resolution velocity models through FWI has further contributed to the improvement of migration image quality [7].

Based on the discussions thus far, FWI is considered conventionally limited to the task of building the velocity model, which is then used as an input for migration. However, Kalinicheva et al. [11] and Zhang et al. [5] developed a technique to directly convert the FWI results into seismic images. This is achieved by taking directional derivatives of the FWI results, referred to as FWI imaging. By directly converting FWI results into reflectivity, the conventional three-step process of preprocessing, velocity model building, and migration can be condensed into a single step (Figure 2). This allows for the rapid derivation of seismic images. Moreover, the use of the full wavefield incorporates additional information lost in conventional migration, resulting in an improved S/N ratio and the suppression of migration artifacts. Additionally, through the least-squares fitting process of FWI, its inherent robustness can be further exploited. Essentially, FWI imaging enables the generation of subsurface images with better structural information and can reduce the required processing time [5].

Huang et al. [7] presented a case study enhancing imaging in subsalt regions with OBN and WATS data, as well as improving the imaging of subgas zones via source-over-spread NATS data using FWI imaging. Similar to the latter case, this article aims to apply FWI imaging to deepwater regions where subsurface imaging issues exist due to shallow gas. FWI up to 25 Hz and FWI imaging were conducted on streamer data, and the results were compared with conventional reverse time migration (RTM) results. An analysis was performed to assess how effectively and accurately FWI imaging can be used to restore the structure of shallow gas-bearing areas. Furthermore, the technical strengths and weaknesses of FWI imaging are discussed.

This section aims to elucidate the methodology concerning FWI, FWI imaging, and RTM. First, a comprehensive explanation of the optimization problem related to FWI is provided to enhance the understanding of FWI. Second, the theory behind FWI imaging is elaborated, focusing on the formulation of FWI imaging and highlighting its advantages. Finally, a discussion of the formulation of RTM for comparison with FWI imaging is presented.

The acoustic wave equation is a fundamental equation that describes the propagation of acoustic waves through a medium, typically within the Earth’s subsurface. FWI begins with forward modeling that simulates seismic data by perturbing specific parameters (in this case, the P-wave velocity while the density remains constant). The acoustic wave equation is used for the wave simulation. While the actual behavior of the Earth is elastic in nature, for the sake of simplification and computational convenience, it is common in geophysical applications to assume that the Earth is an acoustic medium and to utilize the acoustic wave equation. The acoustic wave equation in the time domain can be described as follows.

where $pu$ is the seismic wavefield, $vp$ is P-wave velocity, $fp$ is the source function, $X$⁠, $t$⁠, and $xs$ are the spatial coordinates, time, and spatial coordinates of the seismic sources, respectively.

For forward modeling, this article utilized the fourth-order staggered-grid finite-difference method [12-14].

FWI is a highly effective computational tool utilized to predict subsurface properties of the Earth based on full-waveform seismic data. To predict subsurface properties, FWI relies on an inversion process that updates the model to match the modeled seismic data to the observed seismic data. In geophysical applications, the relationship between the model and data is generally nonlinear. Furthermore, owing to the ill-posed nature of geophysical problems, which fail to yield unique solutions, FWI involves a process of iterative local optimization. Essentially, the misfit function or objective function is defined as the residual between the modeled seismic data and the observed seismic data, and iterative modeling processes are performed to minimize this function (Figure 3) [10, 15]. Typically, the misfit function used in FWI is defined using the $l$2-norm (least-squares norm) for computational efficiency [8].

where $pu$ is the modeled seismic wavefield, $pd$ is the observed seismic wavefield, and $Xj$⁠, $Xi$⁠, and $t$ are the receiver location, source location, and time, respectively.

To obtain minimization, it is necessary to take partial derivatives of the misfit function with respect to the $k$-th model parameter. Based on this gradient of the misfit function, the model can be updated in the direction where the misfit function reaches a smaller value. The gradient of the misfit function is as follows:

where $m$ is the model parameter vector and $mk$ is the $k$-th model parameter.

Based on the gradient of the misfit function obtained in this manner, the model parameters can be iteratively updated according to inversion theory to minimize the misfit function; through this process, the FWI method derives the final model. Based on the steepest descent method, one of the simplest optimization techniques, the updated model parameters can be defined as follows:

where $iter$ is the number of iterative processes and $α$ is the step length.

The gradient direction of the model update equation may vary depending on the optimization technique used.

According to equation 3 above, the term $∂puxj,t;xi∂mk$ within the gradient of the misfit function is referred to as the sensitivity, Fréchet derivative matrix, or Jacobian matrix [10]. The dimension of this Fréchet derivative matrix corresponds to the number of model parameters and data points. Typically, in seismic data analysis, there are many model parameters and data points, resulting in a high dimension for the Fréchet derivative matrix. Consequently, calculating this matrix becomes computationally burdensome.

Therefore, to avoid directly computing the Fréchet derivative matrix, we can utilize the adjoint-state method [16] to modify the gradient of the misfit function with respect to the $k$-th $vp$ as shown below (refer to appendix A):

where $λ$ is the adjoint wavefield, $xk$ is the $k$th grid point, and $T$ is the total duration of seismic data recorded.

When discussing the physical interpretation of the gradient, according to equation 3, the gradient represents the zero-lag cross-correlation between the partial derivative wavefield and the data residual. However, through the adjoint state method, it can be regarded as the zero-lag cross-correlation between the incident wavefield emitted from the source and the back-propagated residual wavefields (equation 3; [10]).

There are various methods for finding solutions in optimization problems. One of the most fundamental algorithms is the steepest descent method. As previously mentioned, the gradient of the misfit function is calculated, and the model parameters are updated in the direction in which the misfit function decreases (equation 4). While the steepest descent method offers computational simplicity, its convergence toward a solution is slow, and it typically updates only the shallow region of the model, leading to less robust results.

In contrast, the full-Newton method updates the solution as follows:

where $H$ is the Hessian matrix, the second derivative of the misfit function.

Unfortunately, computing the Hessian matrix is extremely intensive, particularly when dealing with seismic field data. Consequently, considerable efforts have been made to approximate the Hessian matrix.

One prominent example is the Gauss‒Newton method, which uses an approximate Hessian by neglecting the second term of the Hessian matrix.

where $J$ denotes Jacobian matrix and $†$ denotes the conjugate transpose.

Finally, to reduce computational complexity, one can utilize a further approximation of the approximate Hessian matrix known as the pseudo-Hessian matrix. This approach is termed the pseudo-steepest descent method. Through use of the pseudo-Hessian matrix, it is possible to achieve computational efficiency and correct geometric spreading, resulting in a balanced amplitude distribution between deep and shallow regions [17].

where $S$ denotes modeling operator matrix.

In this article, a preconditioned gradient using the pseudo-Hessian matrix is employed for computational convenience.

Seismic reflectivity is a property of interfaces between two different subsurface layers with contrasting properties, such as density and velocity. Reflectivity is often quantified via the reflection coefficient. In turn, the reflection coefficient is the normalized contrast in the acoustic impedance between the two layers. A typical representation of impedance contrast is as follows:

where the $I$ is the acoustic impedance, which is calculated using density and velocity, $I=ρv$⁠, and $θ$ and $φ$ are dip angle and azimuth angle of the normal vector to the subsurface reflectors. In simplifying the equation, assuming density to be constant or a smooth function leads to the following expression [5].

Through the above equation, the relationship between the reflectivity and velocity can be determined. Ultimately, by taking the directional derivative of the high-frequency velocity obtained through FWI, one can derive the desired reflectivity image. This technique of creating reflectivity images is referred to as FWI imaging, and the images produced through this method can be termed FWI images [5].

Since its inception in the 1980s [18-20], RTM has become widely used as a technique for creating detailed images of subsurface structures of the Earth as computing power has rapidly increased. The basic algorithm of RTM shares an identical structure with FWI. Tarantola [8] and Lailly [9] indicated that RTM corresponds to the initial iteration of FWI using the gradient method. As previously noted, in FWI, the gradient is computed via zero-lag cross-correlation between the observed data and the backpropagated residual wavefields. In contrast, RTM derives the migrated section through zero-lag cross-correlation between the forward propagation of the source wavefield and the backward propagation (reverse time) of the receiver wavefield (Figure 4). This process replaces the residuals in the FWI algorithm with the receiver wavefield [21].

Additionally, RTM assumes that the source wavefield consists of downward-propagating waves and that the receiver wavefield consists of upward-propagating waves. To satisfy these conditions, it would be ideal to use a one-way wave equation in RTM. However, owing to limitations in the propagation angle and the inability to account for complex structures or severe horizontal velocity variations, a two-way wave equation is generally used. The wavefield propagated from the source is computed by performing modeling with an assumed source in the model, and the receiver wavefield is obtained by propagating the recorded waves backward in time through the modeling area via reverse time modeling. This can be expressed by the following equation.

where $us$ is the source wavefield, $ur$ is the receiver wavefield, $Tmax$ is the maximum recording time, $ns$ is the total number of sources, and $np$ is the total number of grid points.

In addition, because the algorithm structures of RTM and FWI are identical, RTM can be obtained via zero-lag cross-correlation of the partial derivatives of the wavefield with respect to the model parameter and the measured data [21].

where $us$ is the source wavefield, $Tmax$ is the maximum recording time, $ds$ is the recorded seismic data, $ns$ is the total number of sources, and $np$ is the total number of grid points. $∂ust∂mi$ is the partial derivative of the wavefield $ust$ with respect to the $i$-th model parameter $mi$⁠. This term represents how changes in the subsurface model affect the wavefield.

However, using the partial derivative wavefield requires a significant amount of computation time. Therefore, RTM can be performed indirectly by using the adjoint property of the wave equation matrix. To achieve this, we can express equation 14 in the frequency domain.

where $ω$ is the angular frequency, $u~sω$ is the source wavefield in the frequency domain, $∂u~sω∂mi$ is the partial derivative of the frequency domain source wavefield with respect to the model parameter $mi$⁠, and $d~s*ω$ is the complex conjugate of the recorded seismic data in the frequency domain.

To obtain the partial derivative wavefield of equation 15, we consider the wave equation in the frequency domain in matrix form.

where $S$ is the complex impedance matrix, $f~s$ is the source function, $K$ is the stiffness matrix, $C$ is the damping matrix, and $M$ is the mass matrix.

When differentiating both sides of equation 16 with respect to the parameter $mi$⁠,

and when expressing equation 18 in terms of partial derivative wavefields, it can be written as follows.

$-∂S∂mius~$ is defined as the virtual source, and for convenience, it is denoted as $fvs$⁠.

When substituting equation 21 into equation 15, the following equation is obtained.

Assuming matrix $S$ is symmetric, the reciprocity property holds, leading to $S-1T=S-1$⁠. Therefore, equation 22 can be summarized as follows:

In conclusion, RTM images can be achieved by multiplying the backpropagation of field data and the virtual source rather than directly using the partial derivative wavefield.

RTM, which is based on the two-way wave equation, demonstrates high imaging accuracy. RTM ensures accurate propagation amplitudes, can relocate any type of multiple reflections to their correct positions, and is capable of accurately interpreting steep dips. It is particularly effective in handling multipath arrivals and has no dip limitations, allowing it to image even overturned energy. Additionally, this technique has the advantage of preserving the true amplitude of the elastic wavefield. Moreover, it generates reflected, refracted, diffracted, and multiple waves in a manner that is consistent with actual wave propagation patterns.

The area selected for the case study is a gas field located in a deepwater region. This area suffers from amplitude attenuation due to shallow gas, so it is difficult to restore the structures in the deeper section. In particular, since conventional methods can use only primary waves, conventional methods have limited ability to improve the imaging quality of subgas zones. Accordingly, this site is selected to validate the extent to which FWI imaging, which uses the full wavefield to derive seismic images, can improve the image quality beyond the limitations of conventional methods.

This case study analyzes the effectiveness of FWI imaging in comparison with conventional RTM. The initial velocity is updated by applying FWI up to 25 Hz. Then, using this same velocity model, both RTM and FWI imaging are performed to derive seismic images, and the results are compared.

This study area is located in the deepwater area of the Bengal Basin, which is notable for its extensive clastic sediments. A significant portion of the sediments in the Bengal Basin is sourced from the Himalayas. Owing to the intense continent‒continent collision between the Indian and Eurasian plates, the Himalayas formed during the middle Eocene. This tectonic activity resulted in an enormous volume of sediment being transported into the Bengal Basin [22].

The northward to northeast oblique subduction of the Indian plate beneath the Burman plate subsequently formed the N–S-trending fold belt known as the Indo-Burman Range. As a result, during the early Miocene, sediment was supplied from the suture zone between the Himalayas and the Indo-Burman Range.

These compressional tectonic forces continued, eventually causing the uplift of accretionary sediments and the formation of thrusts and folds. Consequently, during the late Miocene, sediment was still mainly supplied from the Himalayas, but there was also a significant supply of sediment from this fold belt. As a result, the Bengal Fan extends over 3000 km and reaches a maximum thickness of approximately 16 km, making it the largest submarine fan system in the world [23, 24]. Moreover, tectonic activity associated with the Indo-Burman Range has formed various anticlines and faults in the Bengal Basin, leading to the discovery of gas fields associated with structural and stratigraphic traps. This tectonic history has allowed the area to become a favorable environment for hydrocarbon accumulation through its rich sediments and structural complexity.

The 3D survey utilized in this study is based on towed-streamer data acquired in 2006. The shooting configuration was a dual-source air gun with eight streamers and 400 channels. The shot spacing was 18.75 m flip-flop with a source depth of 6 m below sea level, the streamer separation was 100 m, the streamer length was 5000 m, and the streamer depth was 7 m below sea level (Table 1).

The process of creating the initial velocity for FWI is as follows: smoothing was applied to the legacy migration velocity obtained from a previous commercial company, and FWI was performed using diving wave energy up to 8 Hz. Subsequently, FWI-guided Q tomography was employed to update the Q model to mitigate the amplitude attenuation effects in the shallow gas area, and the resulting outcome was used as the initial velocity.

In terms of modeling dimensions, a total of 1600 × 642 grids were utilized for the velocity model with a grid spacing interval of 12.5 m; similarly, a grid spacing interval of 12.5 m was employed along the depth axis for a total inversion depth of 5 km. Consistent with the characteristics of the FWI imaging method, no specific preprocessing steps were undertaken. To prevent cycle skipping, bandpass filtering was applied to the raw data, and multistage inversion was implemented. Each frequency band underwent sequential iterations: four iterations in the range of 2–15 Hz and five iterations in the range of 2–25 Hz. Given that the initial velocity was derived using only diving waves up to 8 Hz through FWI, this study aimed to utilize both diving and reflection waves up to 25 Hz to construct the impedance structure. Source estimation was performed by matching the near-offset direct wave between the observed and modeled data, utilizing a single average source for the entire modeling process. Additionally, to expedite the computation, only every third source from the shot was utilized during modeling. The offset range used for modeling was set with a minimum offset of 100 m and a maximum offset of 5000 m. Consequently, the majority of the raw data, excluding ultranear-offset data, was employed. To minimize the effects of amplitude differences, the objective function employed was the trace-by-trace normalized least-squares norm as follows:

Ultimately, by running FWI up to 25 Hz, a FWI image was derived by applying spatial derivatives along the depth direction to the velocity model.

In contrast, RTM was performed using the velocity model of 25 Hz, and preprocessing was applied for RTM. The workflow of preprocessing can be seen in Figure 5.

As shown in Figure 6, the updated velocity model obtained through 25 Hz FWI (Figure 6(b)) reveals a higher resolution of structural information than the initial velocity model obtained through 8 Hz FWI (Figure 6(a)). Not only does it delineate the detailed sediment structure in very shallow areas near the water bottom (blue ellipses in Figures 6(a) and 6(b)) but it also effectively provides the geological features of complex geometry in slope areas (blue arrows in Figures 6(a) and 6(b)). The seismic section in Figure 6 illustrates this improvement, highlighting the enhanced resolution and detail achieved through 25 Hz FWI compared with 8 Hz FWI.

Moreover, the 25 Hz FWI velocity not only captures structural information in shallow areas but also effectively delineates structural features in deeper regions. Figure 7 illustrates this feature, with blue arrows pointing to structures at depths greater than approximately 3 km. The initial velocity was derived solely from diving waves via 8 Hz FWI. Owing to the limited reach of diving waves, velocity updates were constrained to depths of approximately 2.5 km, corresponding to half of the streamer length of 5 km. Thus, constraints exist in representing structural information beyond a depth of 2.5 km (blue arrow in Figure 7(a)). However, the 25 Hz FWI velocity incorporates information from reflection waves in addition to diving waves, enabling the representation of velocity details in deeper areas (blue arrow in Figure 7(b)). This finding suggests that improvements in seismic imaging of deeper regions can be attained through updates in higher frequency bands of FWI (Figure 7(c)).

Figure 8 compares the RTM image (Figure 8(b)) and the FWI image (Figure 8(c)), both implemented using the same 25 Hz FWI velocity. In the RTM image, despite the use of the 25 Hz FWI velocity, migration swings and noise—common limitations of conventional migration—are prevalent throughout the section. This noise makes it difficult to achieve continuous strata representation (blue arrow in Figure 8(b)). In contrast, the FWI image exhibits an effective reduction in the noise and significant improvement in the reflector continuity (blue arrow in Figure 8(c)). When the zoomed-in displays (Figures 8(d) and 8(e)) are examined, the marked improvement in lateral resolution and noise suppression in the FWI image is evident.

This finding demonstrates that the results can vary depending on the type of input wavefield used in the imaging process, even when the same velocity model is used. Conventional migration relies solely on primary reflection, whereas FWI imaging utilizes full wavefield data. By incorporating the information from diving waves and multi-scattered wavefields, which are typically discarded in conventional migration, the quality of seismic imaging can be significantly improved. Diving waves, in particular, dominantly move in the horizontal direction, aiding in enhancing the lateral resolution of the imaging process. As a result, compared with RTM, FWI imaging results in improved structural continuity (blue arrows in Figures 8(b) and 8(c)). Additionally, the FWI imaging algorithm, derived through an iterative least square data fitting process, produces more robust results than RTM. This enables FWI imaging to control migration artifacts and achieve higher-resolution images.

Figure 9 allows for a comparison between the RTM image and the FWI image in deeper regions at approximately 3 km depth. Like in Figure 8, by utilizing the full wavefield, FWI imaging results in a continuous structure (arrows in Figure 9), a higher S/N ratio (rectangles in Figure 9), and balanced amplitude characteristics compared with RTM.

Since FWI imaging uses a full wavefield, it provides sufficient illumination from various raypaths. This helps mitigate the acquisition footprints found in conventional migration (Figure 10). Compared with depth slices from RTM and FWI imaging, the RTM image suffers from severe acquisition footprints (Figure 10(a)). In contrast, the FWI image in Figure 10(b), benefiting from the additional illumination information, almost completely eliminates these acquisition footprints and provides more accurate information about the shallow geometry.

FWI imaging is also effective in enhancing imaging in shallow gas-bearing areas. This study area has an abundant distribution of shallow gas, as confirmed by the velocity model in Figure 11(a) (white circle). In the gas-bearing region, RTM results in significant amplitude attenuation and severe migration artifacts due to energy absorption by shallow gas (rectangles in Figures 11(b) and 11(d)). However, FWI imaging yields excellent results in improving imaging and noise suppression (Figures 11(c) and 11(e)). Conventional migration relies solely on primary reflections, making it challenging to acquire structural information from the gas-bearing area, as the primary reflections pass through the gas cloud twice. In contrast, FWI imaging can utilize information from wavefields other than primary reflections, allowing for the acquisition of structural information through various raypaths. Through use of the full wavefield with waves that either avoid or pass through shallow gas only once, FWI imaging can enhance images and identify structures in the subgas area [7]. Consequently, FWI imaging can serve as an alternative to conventional migration for data afflicted with significant seismic imaging issues caused by shallow gas.

FWI imaging significantly outperforms RTM by utilizing the full wavefield, leading to enhanced structural continuity, improved resolution, and superior noise suppression. This allows FWI imaging to deliver clearer subsurface images, particularly in challenging areas like subgas zones and complex subsurface environments. Furthermore, FWI imaging effectively reduces acquisition footprints and migration artifacts, offering more accurate and robust imaging results. Additionally, it bypasses the time-consuming preprocessing and migration steps required in conventional methods, allowing for faster and more efficient image generation while maintaining superior accuracy.

While FWI imaging offers substantial improvements in seismic imaging, FWI imaging has its limitations. FWI images in this study show the enhanced detail and robustness of imaging up to approximately 2.5 km, where diving waves can penetrate. However, deeper sections still present structural ambiguity or resolution limitations. Addressing these challenges requires careful consideration in acquisition design, including the use of longer streamers to capture diverse and broadband full wavefields. Furthermore, inherent limitations in FWI algorithms, stemming from simplified physics, underscore the need for future advancements in computational power to accommodate various parameters or develop more accurate inversion algorithms.

Nevertheless, even with these considerations, FWI imaging consistently outperforms conventional migration methods when the same velocity model is employed. FWI imaging offers distinct advantages, such as a continuous structure, a high S/N ratio, and the removal of acquisition footprints. Additionally, its ability to generate rapid results by bypassing preprocessing and migration steps positions FWI imaging as a promising alternative to conventional migration methods, with the potential to significantly enhance seismic imaging practices in the future.

The data used in this research are confidential and the author does not have permission to share or release to the public.

The authors declare that there is no conflict of interest regarding the publication of the paper.

This work was supported by a grant from the Human Resources Development program (No. 20204010600250) of the Korea Institute of Energy Technology Evaluation and Planning (KETEP), funded by the Ministry of Trade, Industry, and Energy of the Korean Government. J.L. work was supported by the Korea Institute of Energy Technology Evaluation and Planning (KETEP) and the Ministry of Trade, Industry & Energy (MOTIE) of the Republic of Korea (No. 20224000000120). Y.C. work was supported by the Basic Research Project (24, 3809-위탁1) of the Korea Institute of Geoscience and Mineral Resources (KIGAM) funded by the Ministry of Science, ICT of Korea.

We appreciate POSCO International’s permission to use the dataset and publish this work.