Abstract
The assumption of a homogeneous elastic half-space model is widely used to model the earth’s deformation. However, the homogeneous assumption would not accurately reflect the complexity of the shallow crust. We performed a 3D coseismic deformation model using the finite element method and referred to the 2010 Mentawai earthquake. The 2010 tsunami earthquake was located at the Mentawai segment, which is a part of the accretionary wedge in the Sumatra subduction zone. This active accretionary wedge is identified as the most complicated structure on earth and lies along the Sumatra subduction zone, at which most destructive earthquakes happen in this region. We examined the impact of the accretionary wedge geometry and material properties by considering the wedge as a single different property separated from the continental plate. Various geometrical features, such as topography and wedge dimension, as well as physical properties, were simulated. Those features are then observed for their responses on the surface deformation. The topography affected the magnitude of the horizontal deformation up to 10% but only the pattern of the vertical deformation. The wedge dimension seems to have an insignificant influence on the surface deformation compared to the topography. Different physical properties of the accretionary wedge affect not only the magnitude of the horizontal deformation up to 40% but also the orientation. The direction of the lateral movement is seemingly affected by the material under the GPS station and by the source. On the other hand, the variations in the physical properties resulted in discrepancies of 0.5 meters in the vertical deformation near the source. These results indicated that regional physical property information and geometrical features are critical in estimating coseismic deformation, leading to more accurate slip inversion and earthquake and tsunami hazard prediction, particularly in regions with significant inhomogeneity.
1. Introduction
Studies of the earth’s deformation using elastic dislocation have been carried out for over 5 decades [1, 2]. Okada’s [2] formula is widely used to estimate the earth surface deformation by assuming a homogeneous half-space medium model with a flat surface. This assumption may be used when insignificant medium contrast in the vicinity of the source is considered. However, the shallow part of the crust is more complex than a simple homogeneous model. Using a homogeneous elastic half-space model may lead to less accurate prediction in the coseismic simulation [3] and may not be sufficient to obtain the Green function especially for a region where the existence of inhomogeneity is large [4]. Zhao et al. [5] show that the effect of layered rigidity contributes to at least 20% of the surface deformation. However, incorporating a lateral rigidity variation on surface deformation is still rarely done. On the other hand, accumulated sediment called an accretionary wedge exists at some part of the subduction zone around the world. These sediments predominantly come from the oceanic plate deposits, have a wedge-shaped form, and are characterized as the most complex tectonic structure in the world [6]. The accretionary wedge is divided into two parts. The inner wedge consists of older accreted sediment and is located more landward. This inner part is less compressive and more stable than the outer wedge. The outer wedge is near the trench and is more deformable. The base of the accretionary wedge is the location of the seismogenic zone at which most devastating earthquakes happen [7]. Okamura et al. [8] analyzed the source of the giant tsunami Ryukyu 1771 earthquake which was estimated due to the collapsed accretionary wedge. On the other hand, some geometrical features are suggested contributing to the coseismic deformation such as topography and wedge dimension. The effect of topography yields different impacts in some areas and depends on its gradient [9]. A large topography changes likely have a significant effect on predicted coseismic deformation in the case of 2005 Nias earthquake [10] but minimal at the region where significant topographic change is absent for the case of 2009 L’Aquila earthquake [11]. The topography may have a major impact, especially due to slip at shallow depth, on coseismic deformation (e.g., [12]). Moreover, other geometrical features such as wedge dimension are rarely examined. Lotto et al. [13] found that larger compliant prisms enhance shallow slip and larger tsunami. In their simulation, only the outer part of the wedge is considered.
On Sumatra Island, Indonesia, the Sumatra subduction zone runs from the north-west to the south-east between the Indo-Australian plate and the Eurasian continental plate (Figure 1). With a convergence rate ranging from 60 millimeters per year in the southern Sumatra to 52 millimeters per year in the northern Sumatra, the Indo-Australian plate subducts beneath the Eurasian plate in an oblique direction [14]. This oblique subduction induces a strain partitioning into thrust motion along the plate contact and strike-slip motion at the Great Sumatra fault [15]. Unique nonvolcanic forearc ridges that are exposed above sea level can be seen in the Sumatra subduction zone [15–17]. Between the trench and the mainland (the Sumatra Island), this forearc ridge causes the creation of islands (e.g., Simeulue Island at Aceh, Nias Island at North Sumatra, Siberut and Pagai Islands at West Sumatra, and Enggano Island at Bengkulu). The Mentawai segment which is located at the seaward side of West Sumatra consists of Siberut Island, Sipora Island, North Pagai Island, and South Pagai Island. Based on the geometry, the accreted Mentawai segment is classified as compressional wedge [18]. This segment is thought to have been created by mélange complexes that were once a part of an uplifted accretionary wedge [19]. The ultramafic (ophiolitic) basement rocks composed of serpentinite, pyroxenite, and serpentinised dunite are the source of the rock formation and hydrothermal alteration of the oceanic crust and mantle that makes up the Mentawai segment’s oldest rocks. It was indicated that the mélange complexes were formed since Late Eocene. The tectonic activity in the Mentawai segment started in Early Oligocene-Early Miocene and formed the Sipora and Pagai Islands at the first period due to collision between the Indo-Australian plate and the Eurasian plate. The Siberut Island was then formed in the Early-Middle Miocene at the second period.
The Mentawai segment’s earthquakes typically occur every 200 years [20]. However, other sequences (for example, 1350/1388, 1658/1703, and then 1797/1833) have shorter periods every 40 years. In addition, a few earthquakes occurred in a period of fewer than 40 years, between 1500 and 1600. This historical data suggests that the next earthquake might occur no later than 2047. In the years 1797 and 1833, there were two significant earthquakes with magnitudes of Mw 8.5-8.7 and Mw 8.6-8.9, respectively [21]. The 1797 earthquakes produced an uplift of 0.8 meters and a tsunami run-up height of 5 meters. The earthquake of 1833 caused a 2 m uplift and a 3–4-meter tsunami run-up height. A recent large earthquake of Mw 8.4 occurred in 2007 and was followed by Mw 7.9 within a 12-hour period in the same region [22]. Only a small fraction of the 1833 ruptured area was ruptured by these earthquakes. The updip portion of the 2007 area was unexpectedly ruptured on October 25, 2010, by an earthquake with an Mw of 7.8. According to the Global CMT [23, 24], this earthquake was located at -3.71S and 99.32E, with a strike angle of 322°, dip 7°, and rake 98°. The earthquake rupture lasted for ~124 s and had a magnitude . This moderate earthquake, which is categorized as a tsunami earthquake, caused a tsunami that was higher than anticipated and had a run-up height of up to 16.9 meters.
We chose the 2010 Mentawai earthquake as our study area which is located at the Mentawai segment as a part of the accretionary wedge. Over 10-20 GPS stations which are part of Sumatran GPS Array network are installed around this segment and allow us to observe different behaviors of the surface deformation at the accretionary wedge area and the continental plate. Current computational technique, such as finite element modeling, has the advantage of including more complexity of the earth’s structure to estimate the surface deformation such as lateral variation and geometrical features. Therefore, we performed 3D coseismic deformation modeling to investigate those features. Our objective is to observe the impact of topography, wedge dimension, and different rigidity of the accretionary wedge on the predicted surface deformation.
2. Materials and Methods
The surface deformation was calculated in a 3D finite element simulation to observe the effect of all features mentioned in the previous chapter. The computation was conducted using an open-source FEM software PyLith [25–27] from Computational Infrastructure for Geodynamics (CIG). PyLith is widely used to simulate crustal deformation across spatial and temporal scales, especially for quasistatic and dynamic modeling earthquake faulting (e.g., [10, 28, 29]). PyLith implements fault slip using a domain decomposition approach with Lagrange multipliers. At the beginning, a 3D model was meshed with the software Coreform Cubit 2020.2 [30]. The model’s dimensions were 1500 km long, 1700 km wide, and 200 km in depth, covering from 93.2°E to 108.4°E of longitude and 2.5°N down to 6.5°S of latitude (Figure 1). The model had over 800.000 tetrahedral meshes, with the finest size of about 2.5 km at the source and the coarser size up to 100 km at the edge of the model (Figure 2(a)). The earth’s structure was divided into five blocks: elastic oceanic plate, continental plate, oceanic mantle, continental mantle, and accretionary wedge (Figure 2(b)). The continental and oceanic plate width was 65 km and 75 km, respectively [29]. The simulations were then constructed into three categories: first, the effect of topography was investigated by applying flat topography (zero elevation) and topographic relief (bathymetry and topography included) to the top surface of the models. Both models were then given different physical properties: homogeneous model and layered model. Identical physical properties were applied for each material in the homogeneous model (HOM). The velocity profile proposed by Collings et al. [17] at the Mentawai segment was adopted in our model. Vp value ~6.52 km/s, Vs ~3.26 km/s, and density ~2700 kg/m3 were chosen for the homogeneous model. On the other hand, vertical variation of the physical properties or layered model (HET) is shown in Figures 3(a)–3(c). Thus, there are four models in the first category, and we named them HOM-0topography, HOM-topography, HET-0topography, and HET-topography. The homogeneous material was included in the simulation to observe whether the topography effect was consistent for both homogeneous and layered models. The bathymetry and topography data from GMRT with ~500 m grid size [31] and Slab 2.0 [32] were used for the subducting surface. The subducting surface was then interpolated up to the trench with ~6° of slab angle following the seismic profile by Singh et al. [33] and Hananto et al. [34]. In the second category, the dimension of the accretionary wedge was varied to evaluate the impact of different lengths and widths of the accretionary wedge on surface deformation. As the length and width of the accretionary wedge along Sumatra subduction are different, the accretionary wedge was evaluated with length variations of 175 km, 200 km, and 225 km from the trench and 16 km, 21 km, and 26 km of width from zero elevation. Identical layered physical parameter in the first category for all wedge variations (Figure 3) was applied. In the third category, the accretionary wedge was considered as a single material separated from the continental plate to observe its impact on surface deformation. All the materials except the accretionary wedge used layered material as in the previous simulation. The physical parameter of the accretionary wedge was applied as a single value to observe the impact of the lateral variation on the surface deformation. The compressional wave velocity (Vp) varied from 5.1 to 5.7 m/s and 1.6-2.2 for the Vp/Vs ratio, respectively, following the minimum and maximum values from Collings et al. [17]. The values of Vp and Vp/Vs ratio are associated with the shear wave velocity (Vs) and rigidity values, as shown in Figures 3 and 4. The second last and third last categories were simulated using the geometry model including the topography. The slip source solution from Hill et al. [35] was used in all categories. This coseismic slip distribution was obtained by performing standard linear inversion techniques and using the EDGRN/EDCMP code from Wang et al. [36]. The structure was assumed as a horizontally layered crustal model based on the CRUST 2.0 model. Using this slip source model, we compared our results to the real observational data and observed the impact of each parameter in estimating coseismic deformation on the surface. Twelve c-GPS stations from Sumatran GPS Array (SuGAr) [37] were used as reference to be compared with our results.
(a)
(b)
(c)
(a)
(b)
(c)
(d)
3. Results
3.1. First Category: The Impact of Topography
The impact of topography on both horizontal and vertical deformations is observed in the first category. The horizontal displacement due the topography has less magnitude compared to the flat model for both homogeneous and heterogeneous model. The deformation above the source from the homogeneous and flat model gives the largest horizontal movement among all the models (Figures 5(a), 5(b), 5(d), 5(e), 5(g), 5(h), 5(j), and 5(k)). Similar behavior is shown on the GPS offset (Figure 6(a)). The differences between the HOM-topography model compared to the HOM-0 topography ranged between 2 and 40 mm at BSAT, SLBU, PRKB, SMGY, and KTET stations for the horizontal offset. The horizontal offset difference at the rest of the stations is affected by less than 1 mm. Larger discrepancy is resulted by the heterogeneous model on the horizontal displacement. The differences around 20-108 mm at BSAT, SLBU, PRKB, SMGY, and KTET stations, around 20 mm at LNNG and MKMK stations, and less than 1 mm at the rest of the stations resulted from the HET-0topography model. The last model, HET-topography, yields 20-160 mm differences at BSAT, SLBU, PRKB, SMGY, and KTET stations, around 20 mm at LNNG and MKMK stations, and less than 1 mm for the rest of the stations. On the contrary, the vertical deformation magnitude due to a flat surface is slightly larger compared to the model with topography (Figures 5(c), 5(f), 5(i), and 5(l)). The largest vertical deformation above the source is given by the HOM-0topography model (Figure 5(c)). The flat model gives larger vertical movement on both homogeneous and heterogeneous models (Figures 5(c) and 5(i)). However, the heterogeneous model for flat and topography model yields larger vertical movement at the GPS stations (Figure 6(a)) with differences around 1-7 mm for all stations. However, the topography affects the uplift pattern due to the surface irregularity. The cross-section (Figure 7(a)) shows that the magnitude’s uplift is almost the same between both models but different in its pattern.
(a)
(b)
(c)
(d)
(e)
(f)
3.2. Second Category: The Impact of the Wedge Dimension Variation
Different dimensions of the wedge affect the magnitude of the surface displacement but not the direction of the offset. At the top of the accretionary wedge, larger horizontal displacement is given by a longer wedge Figure 6(b), as shown by the red arrow. The same effect is given by a thinner wedge with the same length as pointed by the black arrow (length: 200 km and width: 26 km) as the thinnest, followed by a yellow arrow (length: 200 km and width: 21 km), and the thickest gray arrow (length: 200 km and width: 16 km). The difference given by the widest of the wedge compared to the least (26 km to 16 km) is around 5-20 mm for the horizontal displacement and 1-8 mm for the vertical displacement at BSAT, SLBU, PRKB, SMGY, and KTET stations. The differences less than 1 mm for the horizontal deformation and less than 0.1 mm for the vertical deformation are obtained for the rest of the stations. On the other hand, the largest length of the wedge compared to the smaller shows 1-11 mm of horizontal displacement discrepancies at BSAT, SLBU, PRKB, and SMGY stations and 1-5 mm for the vertical displacement. The differences at the other stations are less than 1 mm for the horizontal offset and less than 0.2 mm for the vertical offset. The vertical displacement at the cross-section (Figure 7(b)) shows differences for all wedge variation around 5 mm in maximum. Based on this simulation, the wedge dimension variation gives small effect on surface deformation at the mainland but quite large effect on the island or close to the source.
3.3. Third Category: The Impact of the Rigidity Variation for the Accretionary Wedge
Since the accretionary wedge is considered a single material, its impact on the surface deformation can be observed, especially at the GPS station. The Vp value variation impacts only the horizontal offset’s magnitude (Figures 6(c) and 6(d)). Higher Vp contributes to larger horizontal deformation and smaller deformation due to lower Vp value. The Vp variations change the magnitude without changing the direction of the horizontal deformation. In comparison to the model with Vp/Vs values of 1.6 and 2.2 (Figures 6(e) and 6(f)), the variation of the Vp/Vs value also affects the magnitude of the deformation at the island and also at the mainland where higher Vp/Vs value yields larger horizontal deformation at the GPS stations. Interestingly, even though the material of the continental plate does not change, the accretionary wedge’s different material properties have an impact on the deformation of the surface on the mainland. It should be noted that the mainland is a part of the continental plate, while the island is a part of the accretionary wedge. The differences of the surface deformation at the GPS stations obtained from the rigidity variation are 8-35 mm (Vp/Vs 1.6, Vp 5.1-5.7), 3-23 mm (Vp/Vs 2.2, Vp 5.1-5.7), 8-23 mm (Vp 5.1, Vp/Vs 1.6-2.2), and 12-35 mm (Vp 5.7, Vp/Vs 1.6-22) for the horizontal displacement, respectively. The differences of the vertical displacement are 1-5 mm for the Vp variation and 5-50 mm for the Vp/Vs ratio variation. The higher Vp/Vs almost fits the GPS offset at the mainland, as shown in Figure 6(d). Meanwhile, the uplift at the GPS station seems to agree with the model from Vp/Vs ratio equal to 2.2. The lower Vp/Vs ratio leads to an uplift while the higher Vp/Vs ratio leads to subsidence. The difference on the vertical deformation at the cross-section (Figures 7(c)–7(f)) due to the Vp variation is very small (~0.006 m) compared to the Vp/Vs variation (0.6 m). The variation of Vp/Vs with a constant Vp value shows the rotating of the horizontal displacement at the GPS station. Higher Vp/Vs value leads to the horizontal deformation reaching the GPS offset, as shown by the yellow arrow in Figures 6(e) and 6(f). But at some GPS stations (BSAT, PRKB), the offset direction is not changing. The impact of different rigidity actually occurred at the first category simulation. The model with heterogeneity gives smaller horizontal displacement compared to the homogeneous model regardless the flat and topography model. The heterogeneity seems fail to fit the GPS data on the island but fit better in the mainland (Figures 6(c)–6(f)). On the contrary, even though the homogeneous model fit better on the island, it gives much larger offset or overshoot compared to the GPS data (LNNG and MKMK).
4. Discussions
The irregularity of the surface from this simulation affects its magnitude and pattern. The ground deformation is larger on the horizontal direction for the homogeneous model. A larger magnitude of the horizontal displacement is obtained from the flat model compared to the topography model in both homogeneous and heterogeneous models (Figure 6(a)). [9] simulated a 3D coseismic deformation using the spectral-infinite element method for the 2010 Maule-Chile earthquake with a magnitude of ~8.8. The flat mesh model generates larger displacement, especially on land. The difference between the topographic and flat models is up to 30%. They argued that the Green function obtained from the flat model, used in a finite fault inversion, might lead to incorrect displacement shape on the slip distribution. [10] obtained different behavior which previously observed the same feature in the northern part of our study area. Their results show that the topographic model produced larger horizontal displacement. However, the differences given by the topography model are around 2-40 mm (equal to ~10%) (HOM-topography) and 20-160 mm (equal to ~40%) by the topography model including heterogeneity (HET-topography) at the stations in the vicinity of the source. Including the topography to predict the coseismic offset especially for the earthquake happened at shallow depth should be considered. Large discrepancies of the GPS offset near the source (on the island) and far from the source (mainland) show that the topography must be included into the calculation when using the data nearby the source. In addition, in our case, the flat surface gives higher uplift and lower subsidence (Figure 7(a)) even though the uplift amplitude is slightly larger at the GPS stations for the heterogeneous material for both flat and topography models. Since the topography shows the realistic geometry with a deeper surface near the trench and higher as it goes landward, the topography still contributes to more accurate vertical deformation. Thus, the topography needs to be considered to calculate the deformation especially related to the tsunami prediction.
Different responses to the offset on the mainland and the islands may provide us with some information regarding the effects of the wedge’s dimension. A wider wedge () results in a higher horizontal offset, and a narrower wedge causes a smaller offset on the island. While the remaining stations have variances of less than 10%, those closest to the source of the wedge length variation have discrepancies of 10% or more. However, the wedge’s width accounts for a differential of 30% at the station closest to the source and only less than 10% at the other stations. The vertical offset, on the other hand, responds differently depending on its location. The wedge topography and the station’s placement, which is related to the wedge’s edge, could be related. The vertical offset will behave for all wedge dimension variations based on its location, in contrast to the consistent horizontal offset behavior at the station point. In the meantime, the station on the mainland is showing a slight impact on surface deformation. Given that the difference between all the offsets is so minor except for the width variation, the wedge dimension appears to have a minimal impact on the surface deformation (Figure 6(b)).
As a single accretionary wedge is considered, it shows that lateral variation on the material affects the surface deformation and significantly impacts a specific area where contrast material exists. Higher Vp leads to higher rigidity at which a larger horizontal offset is obtained from the simulation (Figures 6(c), 6(e), and 6(f)). On the contrary, lower rigidity leads to higher vertical deformation. [38] applied low rigidity to the shallow sediment to increase the slip and match the tsunami data for the 2010 Mentawai earthquake. This result verified that low rigidity contributes to greater vertical deformation changes but smaller horizontal displacement. This result also confirms the cause of smaller horizontal deformation at the first category (topography-HOM vs. HET). The lateral offset at the station due to heterogeneity for both flat surface and topography models is smaller than the homogeneous model. The rigidity value for the homogeneous model is ~30 GPa; meanwhile, the rigidity value at the shallowest part of the heterogeneous model is 15.5 GPa. Thus, the regional information is necessary and may impact the surface displacement. Furthermore, the variation of the Vp value only affects the magnitude of the horizontal offset without changing its direction (Figures 6(c) and 6(d)). The rigidity value slightly changes due to different Vp values, but it is insignificant (up to 3-10 GPa). In addition, a constant Vp value with a different Vp/Vs ratio produces a significant change in rigidity value (up to 10-20 GPa). This contrast change of rigidity seems to affect the response at the surface. In this case, the horizontal offset caused by low rigidity Figure 6(d), Figures 6(e) and 6(f) (yellow arrow) is directed toward the source and reaches the GPS data (SLBU, SMGY, KTET, PRKT, and PPNJ). It should be noticed that for all the accretionary wedge rigidity variation models, none of the models can match the direction of the GPS data at station BSAT and PRKB.
Interestingly, different physical parameters applied to the accretionary wedge also affect the surface of the continental plate. The deformation response on the mainland for each variation is the same as on the island but with a smaller value. Lower rigidity of the accretionary wedge causes a smaller horizontal offset at the station, and higher rigidity leads to a larger offset even though the physical parameter of the continental plate is constant. Small topography gradient contributes only 5% difference on seismic potency which is smaller compared to heterogeneity effect (up to 50% on seismic potency) [39]. The deformation right at the top of the accretionary wedge is certainly impacted. But the fact that the offset on the mainland is also affected by this different material has to be pointed out. The overshoot at the mainland shown from the first category by the homogeneous model but fit offset by the heterogeneous model (Figure 6(a)) seemingly shows an existence of different rheology between the accretionary wedge and the continental plate.
Whether its deformation on the mainland is related to its response due to different rheology of the accretionary wedge may be considered. The direction of the offset at the GPS station is another interesting response. As the Vp/Vs ratio gets larger, leading to a lower rigidity, the offset reaches the GPS data. Is this direction related to the material’s rheology in the subsurface in the vicinity of the station but also the source? Will the deformation of the accretionary wedge give a different response from those at the island at some specific rheology? The elastic properties especially at the upper plate of subduction zone dominantly control large shallow slip and slow rupture which is related to a larger coseismic surface deformation. Incorporation realistic megathrust elastic properties and rigidity variation of the upper plate cannot be neglected and contribute to a better earthquake and tsunami hazard prediction [40].
5. Conclusions
We have conducted a 3D crustal deformation modeling for the 2010 Mentawai earthquake using the finite element method. The Mentawai segment is a part of the accretionary wedge at which the 2010 Mentawai earthquake was located. The accretionary wedge is mentioned as the most complex tectonic structure on earth. By considering the accretionary wedge as a different material from the continental plate, we observed its impact on coseismic deformation. Including material heterogeneity especially on lateral direction to simulate the coseismic deformation is still rarely done. Some geometrical features were also added in the modeling such as topography and wedge dimension. It is shown that the topography significantly contributes to the horizontal displacement compared to the vertical deformation, and it should be considered for the slip inversion modeling. Even though the impact on vertical deformation is insignificant, the topography impacts the pattern of the deformation. Including topography in the simulation leads to a more accurate predicted coseismic uplift. The wedge dimension dominantly impacts the horizontal deformation at the wedge surface but has a very small effect on the continental plate. A larger accretionary wedge produces a larger horizontal offset and small scale on vertical deformation. The wedge’s dimension is considered to have minor impact except the wedge’s width. Lastly, the lateral material variation produces a more specific response depending on its location. Higher rigidity of the accretionary wedge leads to a larger horizontal offset at the wedge surface and continental plate surface but lower vertical displacement. The lateral variation impacts not only the magnitude of the offset but also its direction. For some areas, a specific rigidity value produces a predicted coseismic offset that matches the GPS data, but some areas cannot fit the GPS data. The Sumatra subduction zone is one of the accretionary wedges in the world, which are shown up above the surface. This gives us the advantage of knowing better the behavior of the accretionary wedge and its impact on any natural phenomena. These simple simulations can show different behaviors between materials in the real case of earth deformation. The detail of the accretionary wedge geometry and regional information is important to be provided and should be included to estimate the coseismic deformation. More complex geometry and materials can be incorporated using finite element modeling. These information can improve the predicted coseismic deformation and, thus, can be used to model more accurate slip inversion especially the region in which the existence of the inhomogeneity is significant. Future studies are still necessary, such as investigating source slip models which consider heterogeneity and geometrical structures for the 2010 Mentawai earthquake.
Data Availability
The slab geometry, topography, GPS coseismic offset, 2D velocity model, and slip source solution data supporting this study are from previously reported studies and datasets, which have been cited.
Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.
Acknowledgments
The academic financial support of the first author was funded by the Indonesian Endowment Fund for Education (LPDP) provided by Ministry of Finance, Indonesia (contract no: PRJ-1633/LPDP.3/2016). This research was also partially supported by the World Class University (WCU) Program-ITB 2019 and Research Innovation and Community Empowerment (PPMI)–Faculty of Mathematics and Natural Sciences ITB 2022, Indonesia. The authors would like to thank Laboratory of Modeling and Inversion and Laboratory of Physics of Porous Media and Fluid Dynamics, Physics of Earth and Complex System, ITB, for the permission to use its facilities in support of this study. The first author would like to give thanks to Farhan Javed for supporting the computational method.