Unraveling the role played by a buried mud diapir: alternative model for 2016 M w 6.4 MeiNong earthquake in southwestern Taiwan

Southwestern Taiwan exhibits multiple fold‑thrust systems as a consequence of the interaction between the Philip‑ pine Sea Plate and the Eurasian Plate. A prominent geological feature of this region is the extensive layer of GuT‑ ingKeng mudstone, with a thickness of approximately 4 km, which serves as a source material for the formation of mud or shale diapirs. The 2016 Mw 6.4 Meinong earthquake, striking southwestern Taiwan at a depth of 15–20 km and inducing approximately 100 mm of uplift, has prompted investigations into the potential involvement of shallow structures (< 4 km) in this uplift. Recent studies have proposed that such shallow structures may have contributed significantly to the observed uplift during the earthquake. This study aims to elucidate the role of buried mud diapirs in the context of coseismic deformation. Here, we present a modeling approach that utilizes sill‑like dislocations to simulate the deformation at the upper tip of the diapir. Our results indicate a vertical opening of approximately 60 mm at a depth of 1.4 km, which closely aligns with the spatial distribution of tomographic and gravity anomalies. We also examine how the coseismic stress changes induced by the Meinong earthquake can lead to a dilatational strain of about 1.2 microstrain within the shallow depth range of 0–4 km, resulting in extension within our modeled region. In contrast, the dilatational strain diminishes from 0.2 to − 1.2 microstrain at greater depths (4–8 km), implying compression in the subsurface beneath the diapir’s top. This study discusses the potential mechanisms how fluid‑rich and high‑pressure mudstone may be deformed through coseismic process and how mud diapirs may contribute to additional deformation within the seismic cycle.


Introduction
Identifying seismic sources within active mountain belts is important to understand the mechanism of orogenic earthquakes, which is also the key for seismic hazard assessments (e.g., Bilham et al. 2001;Burchfiel et al. 2008;Avouac et al. 2015) and long-term topographic evolution of orogens (e.g., Avouac 2007;Elliott et al. 2016;Whipple et al. 2016).As one of major active mountain belts in the world, the Taiwan mountain belt is a spotlight for orogenic studies and has experienced several earthquakes in recent years.The Taiwan mountain belt is the product of arc-continent collision between the Eurasian and Philippine Sea plates (Teng 1990) with southward propagation of collision (Suppe 1981).However, while the mechanism of the mountain building is thought to be dominated by the critical-taper wedge mechanics and thin-skinned model (e.g., Suppe 1987), the interpretation of the deformation in southwestern Taiwan is nonunique, and several geological and geophysical observations also suggest the existence of other local processes as diapirism or fault-related folding (Ching et al. 2016;Tsukahara and Takada 2018;Mai et al. 2021;Wang et al. 2022).
Since 2010, the SW Taiwan had experienced several moderate earthquakes, such as the M w 6.2 Jiashian earthquake in 2010 (Ching et al. 2011;Hsu et al. 2011;Rau et al. 2012;Lee et al. 2013;Lin et al. 2016), the M w 5.9 Wutai earthquake in 2012 (Chiang et al. 2016), and the M w 6.4 MeiNong earthquake in 2016, in which the MeiNong earthquake (6 February 2016, 03:57 UTC + 8:00) was the most devastating seismic event since the M L 6.1 Hsinhua earthquake in 1946 (Chang et al. 1947;Bonilla 1975).Its strong ground motion (Lee et al. 2016) and associated soil liquefaction (Lu et al. 2017) resulted in many damaged structures.In total, 117 people passed away mainly due to the collapse of residential high-rise buildings.
Comprehensive studies have been conducted and various features were addressed regarding the MeiNong earthquake in 2016.Here two major points related to our study are summarized.First, the fault ruptured on a previously unknown plane buried deep in the middle crust.The centroid moment tensor (CMT) reported by the Central Weather Bureau (CWB) of Taiwan and Global CMT (GCMT; Dziewonski et al., 1981;Ekström et al., 2012) both suggested a near-west striking and gentle north-dipping fault plane, with a strike, dip, and rake of 274.81°,41.74°, and 17.02° and 279°, 22°, and 21° according to CWB and GCMT data (IRIS 2016), respectively.The fault geometry is a gentle plane located at 15-20 km depth with maximum slip 1.2 m.In the map view, the slip appears to be nearly westward motion.(Lee et al. 2016).Second, high-resolution satellite geodetic observations indicated that several shallow tectonic structures were triggered by the earthquake (e.g., Huang et al. 2016;Le Béon et al. 2017;Tsukahara and Takada 2018;Rau et al. 2022).
This study, on the other hand, employed a different perspective of mud diapirism that has been considered as a potential mechanism for interseismic deformation (Ching et al. 2016;Tsukahara & Takada 2018).In SW Taiwan, offshore mud diapirism occurs when a compressional tectonic environment meets rapidly deposited, over-pressured, gas-bearing, thick sedimentary (Kopf 2002;Chen et al. 2014;Doo et al. 2015).This sedimentary layer is known as late Plio-Pleistocene GuTingKeng formation, a thick (> 4 km) mudstone interbedded with thin sandstones (Hsieh 1972;Huang et al. 2004;Sun et al. 2010;Le Béon et al., 2017).The compressional environment originates from the ongoing subduction-continent collision between the Philippine Sea Plate and Eurasian Plate (Mouthereau et al. 2001;Rau et al. 2012), with a convergent rate of approximately 8.2 cm/yr (Lin et al. 2010).
The concept of "mud diapirism" has been renew in recent studies, it has been associated with the term "mobile shale".According to Soto et al. (2021b)'s compilation, the shales represent all fine-grained mud-rich units, which can also be related to the GuTingKeng mudstone in our study area.In advance, Soto et al. (2021a) proposed a novel mechanism to explain how critical state "mobile shale" may act as penetrative (visco-) plastic deformation.This theory has broadened our understanding of mud/ shale diapirism, suggesting that within a detachment fold (or other folding structures), as shear strain accumulates, certain fine-grained materials can reach a critical state, enabling them to exhibit flow behavior.This phenomenon can manifest in various geometries in seismic imaging, including shale sheets, shale diapirs, and shale-cored folds.The dominant type of diapiric structure in our study area remains unclear, and we aim to address this question in the present work.
The relationship between offshore and onshore mud diapirism in SW Taiwan has been the subject of ongoing debate, with several unresolved questions.Key inquiries include whether offshore diapirs are directly connected to inland anticlines, whether inland diapirs exhibit behavior similar to their offshore counterparts, and how inland diapirs interact with folding and faulting processes.Some of these issues have been investigated in multidisciplinary studies (e.g., Hsieh 1970Hsieh , 1972;;Sun and Liu 1993;Huang 1995;Liu et al. 1997;Chen and Liu 2000;Chuang 2006;Yang et al. 2006;Chen et al. 2014;Doo et al. 2015;Giletycz et al. 2019;Wang et al. 2022;Lo et al. 2023).
Notably, investigations such as those by Doo et al. (2015) and Lo et al. (2023) have revealed that both offshore and inland diapirs exhibit high residual gravity anomalies compared to the surrounding rock formations, indicating that the high-density mudstone (0.1-0.2 g/cm 3 higher) may be a product of diapiric processes.In the case study of Liuchiu Island, an offshore diapir island, Wang et al. (2022) successfully modeled the geological uplift rate using a numerical model based on fault-related folding.They concluded that mobile shale or mud diapirism should be considered a secondary process in this context.On the other hand, Ching et al. (2016) reported the existence of an inland diapiric anticline with a rapid interseismic uplift rate of approximately 18 mm/yr, posing difficulties for modeling through a 2-D fault slip model, which suggests a potential contribution from diapirism.For coseismic scenario, Rau et al. (2022) utilized high-rate GNSS and strong motion seismograms to model the seismic source of 2016 Meinong earthquake, in their study, after compared the model prediction to coseismic displacement field derived from GNSS daily solutions, significant residuals (~ 80 mm) lead to an interpretation of triggered anelastic and aseismic deformation related to mud diapir nearby Gutingken anticline.
These findings underscore the need for further research to resolve the complex and multifaceted issues about mud diapirism in southwestern Taiwan.We analyze the residual deformation patterns and carry out a 2-step strategy to model the deformation caused by diapiric source triggered by Meinong earthquake.Through stress analysis and comparing with geophysical observation, we unravel how a fluid-rich flat-top diapir can be deformed by coseismic stress changes which may amplified coseismic deformation in SW Taiwan.Especially, where several important infrastructures were announced and claimed as urgently needed, such as candidate sites for MeiNong Dam and Industrial Waste Landfill in Fig. 1.

Data acquisition and processing
Four Copernicus Sentinel-1A SAR acquisitions were obtained through Alaska Satellite Facility (S1A, European Space Agency satellite; Copernicus Sentinel data 2016.Retrieved from ASF DAAC 2016, processed by ESA.).These acquisitions contained two pairs of SAR images, ascending (2-14 February 2016) and descending (4-16 February 2016) pairs, respectively.The differential interferometric SAR (DInSAR) results were generated with commercial software GAMMA (Wegmüller and Werner 1997).The digital elevation model used to subtract topographic phase in the interferograms was the 1 arcsecond resolution (~ 30 m) Shuttle Radar Topography Mission version 3 data (Farr et al. 2007).The interferograms were obtained with a resolution of approximately 30 m 2 and were unwrapped with Snaphu open-source software (Chen and Zebker 2002).The major feature of the InSAR result is a strong line-of-sight (LOS) range decrease (~ 120 mm) concentrated in an oval-shaped region ~ 20 km east of the epicenter, which stops with sharp edges that are east bounded by the Lungchuan ridge and west bounded by a north-south striking structure near Guanmiao (Fig. 1 and Additional file 1: Fig. S4).
Continuous global positioning system (GPS) coseismic displacement field data were also collected to compare with the InSAR results (Additional file 1: Figs.S1, S2 and S4).The coseismic GPS displacement field was the same as that used by Lee et al. (2016), compiled by CWB.Coseismic displacement is defined as the difference between 7 days' average positions before and after 6 February 2016.All GPS solutions were processed using GAMIT (Herring et al. 2015).Because only one GPS station is located within our uplifted area of interest (Additional file 1: Fig. S2b), however, GPS data were not used in the following source modeling.

Modeling for coseismic triggered diapirism
The main fault slip model of 2016 Meinong earthquake Three MeiNong fault slip models were proposed without including coseismic DInSAR data (Chuang et al. 2016;Lee et al. 2016;Rau et al. 2022), and here we compare them and choose the one with its predicted coseismic displacements closest to the DInSAR observations.Theoretical coseismic displacements were calculated by assuming the fault dislocation buried in an elastic halfspace (Okada 1992).Our chosen model is Lee et al. (2016; see comparison in Additional file 1: Fig. S3) that was based on a joint inversion of GPS displacements and seismic waveforms.Their fault slip model was obtained with strike 281° and dip 24° plane buried at 15-20 km depth.The major asperity with 1.2 m westward slip was located at 20 km.The LOS residuals after this model fitting are up to 50 mm, which implies that half of the coseismic deformation may not be associated with the main fault rupture.Moreover, the areas of the high LOS residuals shown in the S1A ascending and descending scenes (Fig. 2) suggest that these motions were mainly vertical (uplift), and residual GPS displacements also confirmed this perspective (Additional file 1: Figs.S2, S11).
Coincidently, the location and pattern of this residual LOS displacement are very similar to those of a rapid interseismic (2007)(2008)(2009)(2010)(2011) deformation observed by ALOS-2 InSAR (Fig. 2), with LOS rates up to 37 mm/yr, which was suggested mainly driven by subsurface mud diapirism (Tsukahara and Takada 2018).Accordingly, we treated this residual deformation as the coseismic motion induced by a subsurface mud/shale diapirism as the secondary source.

Model assumptions of diapiric source
On the geometry and mechanism of mud/shale diapir, Hudec and Soto (2021) provide a catalogue by reviewing diapirs in seismic reflection.Considering the regional tectonics of SW Taiwan, the passive, active and thrust piercement could be candidates for our study area.In their study, most of mobile shale mechanisms form a flat top, and the piercing mechanism favors both extensional and compressional stress environments.
From the perspective of geodetic modeling, an uplift dominated deformation generally prefers flat sources, such as tensile dislocation (Okada 1992) or penny-shape crack (Fialko et al. 2001) in an elastic half-space.Because these types of flat geometry generate less horizontal motion other than spherical sources (e.g., Mogi 1958;Yang et al. 1988).Therefore, we adopted a horizontal, sill-like tensile dislocation (Okada 1992) to model the opening on the top of a buried diapir.The dislocation model was often used for analyzing volcanic sources; for example, Wicks et al. (2006) and Chang et al. (2007)  We recognized that inelastic rheology would be considered for a mud-rich area such as SW Taiwan.In the case of the MeiNong earthquake, however, high-rate GPS data suggested that most of the coseismic deformation was formed in a short period, perhaps in minutes or less (Huang et al. 2016).Rau et al. (2022) conducted an exhaustive analysis of 1 Hz High-Rate GNSS time series, revealing the occurrence of two distinct stages of uplift.The onset of the second stage is identified to start at 7 s from the origin time then reached stable within a minute.Such short-term strain loading led a plausible assumption that the coseismic deformation was mainly elastic.Together with the evidence of low plasticity (average 14 ~ 18%) for mudstones in SW Taiwan (Lee et al. 2007), we proposed that tensile dislocation (sill-like dislocation in following text) in an elastic half-space can serve as a plausible working model for our study.

Two-step modeling strategy for diapiric source
To perform the source modeling, we introduced a twostep approach that includes nonlinear optimization of dislocation geometry and nonnegative least-square inversion for sub-patched openings.In the first step a MATLAB-based open-source software called Geodetic Bayesian Inversion Software (GBIS; Bagnardi and Hooper, 2018) was used to determine the optimal location, length, width, depth, strike, and uniform opening of a sill-like dislocation.The GBIS implemented a Metropolis-Hastings algorithm (Metropolis et al. 1953;Hastings 1970) based Markov-chain Monte Carlo (MCMC) method with adaptive steps to rapidly calculate the posterior probability distribution of model parameters.
The single-patch sill-like model is mechanically different from the class of pressure (or stress) driven crack model often used for fluid or gas filled sources (e.g., Mogi 1958;McTigue 1987;Yang et al. 1988;Fialko et al. 2001).While the top of a mud diapir was often addressed as an overpressure fluid-rich reservoir (see also later discussions), one may argue if the above single-patch sill-like dislocation can be adequate for our case.Davis (1983) demonstrated that, however, a rectangular sill-like dislocation can provide a computationally efficient method to predict similar deformation produced by a pressurized 2-D elastic crack model, despite the fact that real cracks tend to be elliptical in shape with non-uniform opening.
To reconcile these different model design, we adopted the above optimal single-patch sill-like dislocation model depth but subsequently evaluate its opening distribution by subdivided the model geometry (the second step).
After data preparation (Additional file 1: Text S1), we ran the MCMC sampling algorithm within plausible ranges of model parameters that were set fairly wide to explore the posterior probability (Additional file 1: Table S1).We also performed MCMC sampling several times to ensure result stability.Approximately 200,000 samples are required to accomplish the burn-in stage, and a million samples generally result in convergence for each parameter.In some extreme scenarios, if the deformation pattern is discontinuous with distant local maximums, a single-patch sill-like model may be exhausted to find best approximation.In our case, the residual deformation is rather concentrated in a small region, Additional file 1: Fig. S9 shows our optimal model succeed to fit the major deformation pattern on both ascending and descending LOS data.
Our best MCMC results are shown in Additional file 1: Table S1 and Additional file 1: Figs.S7-S9, which indicate a single-patch sill-like dislocation with uniform opening (48 ± 2 mm) at a depth of ~ 1.4 ( 1.4 +0.17 −0.14 ) km.The uncertainty from single-patch model assumption was detected through MCMC sampling process, it's about 0.1 km in depth and 2 mm for opening.A depth uncertainty of 0.1 km is negligible since it only alters the predicted uplift by less than 0.5 mm in the second step.
The second step, with optima depth at 1.4 km, we then subdivided the single-patch model, into 1 km × 1 km patches (Fig. 3) and performed the nonnegative leastsquares to invert the LOS residual deformation for the best-fit opening on each patch.This inversion problem is formulated as where d is the observed residual data (Fig. 2c and h), G is the Green's function formed with the dislocation solution, α is the factor controlling the smoothness incorporated with second-order Tikhonov regularization matrix L, and m is model matrix for opening.The smooth factor α was obtained by grid searching with evaluation of the trade-off between variance reduction (VR) and α.The definition of VR follows that of Huang et al. (2016), and is expressed as where s i is model prediction for ith data d i .The preferred model (Fig. 3) is satisfied when α is 0.6 and VR reaches 48.6% (Additional file 1: Fig. S10).

Model result
The InSAR residuals were significantly reduced by our diapir source model from < 50 mm to < 20 mm (Figs. 3, 4a and Additional file 1: Fig. S4), comparable with the shallow thrust model proposed by Huang et al. (2016), which suggests that both models are plausible for interpreting the secondary coseismic deformation of the MeiNong earthquake.The optimal sub-patched openings are up to 66 mm and clearly bounded by a sharp edge in the east (Fig. 3).Here we refer this feature to a shallow fault developed at the top of the diapir (Schultz-Ela et al. 1993;Dooley et al. 2009).According to the location and orientation, this fault may occur along the previously proposed Lungchuan backthrust (Le Béon et al. 2017).

The relationship between diapiric source and geophysical data
Seismic tomography is widely used to detect the fluid content in lithosphere.Huang et al. (2016) reanalyzed V p /V s ratio in previous published tomography (Huang et al. 2014) and pointed out a high V p /V s (~ 2.0) volume, which could indicate a fluid-rich region above the main fault plane of the MeiNong earthquake.However, their study included limited seismic stations, which may lead to insufficient resolution for analyzing the physical property of the rock volume.Another local seismic tomography recently conducted by Kuo-chen et al. (2017), with a deployment of 36 temporary stations for one month following the MeiNong earthquake, otherwise provided a high-resolution shallow velocity structure with ambient Rayleigh wave.Remarkably, a low V s anomaly (~ − 20% with respect to average value of the tomography) was found at depths of 0-2 km right under the InSAR uplift area and is consistent with the location of our modeled sill-like diapir source (Figs. 2, 3).Considering the evidence that the mudstone has generally low permeability (1.0 ~ 25 × 10 −8 cm/s) and much higher porosity (20-60%) at depths above 2 km (Lee et al. 2007;Mondol et al. 2007), the low V s anomaly could be interpreted as an isolated fluid-rich and potentially over-pressured mud reservoir (e.g., Bonini 2012;Mazzini and Etiope 2017;Ijiri et al. 2018) within the shallow GuTingKeng mudstone formation.
Another evidence from residual Bouguer gravity anomalies analysis made by Lo et al. (2023) clearly revealed a series of high gravity residual from 0-5 km significantly correlated with anticlines in SW Taiwan.Their work shows the gravity residual associated with Gutingkeng anticline (GTKA in Fig. 1c) is slightly lower than other anticlines, which is interpreted as "water and gas containing mudstone" or relative buried.
The strengths of mudstone are notably influenced by fluid contents.For air-dried specimens from SW Taiwan, triaxial compression tests revealed that their brittle failures were greatly suppressed with confining pressures higher than 40 MPa (Lee et al. 2007), which correspond to depths of > ~ 1.6 km if the density of 2600 kg/m 3 was taken (Doo et al. 2015).Lee et al. (2007) also found that the specimen strength can be further reduced by increasing water contents (WC), for example the uniaxial compressional strength drops from ~ 10 MPa for dried (WC = 0.1%) samples to ~ 4 MPa for wet (WC = 5%) ones.Background seismicity also support the above evidence where most of the pre-and post-MeiNong earthquakes in our study area occurred deeper than ~ 6-8 km, or the thickness of the GuTingKeng formation (Tsai et al. 2017).These observations all suggest that our coseismically expanded diapir, with the top depth of ~ 1.4 km, is likely under an environment where mudstone deforms mainly as a ductile material and may also relate to the rapid interseismic deformation in the vicinity (Ching et al. 2016;Tsukahara and Takada 2018).2017) and Gutingkeng anticline (Wang et al. 1998;Chen 2013).Where sudden changes of final residuals occurred are denoted as Lungchuan, east-dipping fault and west-dipping backthrust merged at surface.b Cross-section after Mai et al. (2021).Mud diapirs (MD) are outlined with dash lines.Here we re-interpreted these diapirs with mobile shales in the center of folds, denoted as gray area.The opening region from sill model proposed by this study is denoted as light green region, its location overlapped with MD-3 diapir from Mai et al. (2021) and extend to west.The main fault plane (thick gray line) and hypocenter (red star) of 2016 Meinong earthquake located at 15 km in depth

Coseismic strain and stress analysis for diapiric source region
We also carried out a coseismic stress change analysis to evaluate how the main fault slip can mechanically affect the shallow diapir opening.Bonini et al. (2016) proposed a feeder dike concept to analyze the role of coseismic normal stress changes in mud diapirs or mud volcanoes, which assumes a vertical planar dike as a pipe to supply muddy material or as a root of the mud diapir.Mai et al. (2021) applied this concept to investigate the coseismic normal stress changes along a series of N-S trending vertical planes.Their analysis led to the conclusion that a mud diapir beneath the Gutingken anticline (MD-3 in Fig. 4) could be triggered.Additionally, Mai et al. (2021) computed Coulomb stress changes on surrounding shallow faults, such as the Lungchuan fault and Lungchuan backthrust.However, the observation of negative Coulomb stress changes suggested that none of these faults were likely to be triggered.We followed the same concept to evaluate the static normal stress and dilatation changes induced by the MeiNong mainshock on a series of vertical planes beneath the coseismic uplift area.Coulomb 3.3 (Lin and Stein 2004;Toda et al. 2005) was used to calculate the coseismic stress change in a homogeneous elastic half-space induced by the fault slip model of Lee et al. (2016).
Figure 5a-c shows the dilatational strain changes at different depths, in the sill model region (black box in Fig. 6) volumetric increases to 1.2 microstrain (10 -6 ) at shallow (< 2 km) while decreases at deep depths are revealed (0.2, − 0.5, − 1.2 microstrain at 4, 6, 8 km depth).Changes of normal stress, moreover, show a consistent positive normal stress pattern for north-striking vertical receiver fault planes (Figs.5c, d, 6d-f ), unclamping to 0.1 MPa in depths of 0-4 km.Then switch to negative normal stress (clamping, − 0.2 MPa) in deeper depths.Similar results were found when we altered the strikes of vertical receiver fault planes to east-west (Fig. 6g-i).The overall figures of these coseismic stress analysis can be summarized as follows: vertical cracks experienced a lateral expansion right beneath the uplifted area, which changed to compression at depths exceeding 4 km.Considering pre-existing subsurface mud diapirs, the scenario is similar to squeezing toothpaste while loosening the cap (Dooley et al. 2009;Mai et al. 2021).

Reconciling the role of coseismic folding and diapirism
Coseismic folding is another mechanism often considered when uplifted/tilted regions coincide with existing anticlines; for example, the 1983 Coalinga earthquake (Stein and King 1984), the 1987 Whittier Narrows earthquake (Lin and Stein 1989) and the 2015 Pishan earthquake (Li et al. 2016;Wu et al. 2019).When a coseismic fault slip occurs along a fold-related thrust, the deformation can contribute to folding (Stein and Yeats 1989).However, the fault slip of the MeiNong earthquake is buried deep at 15-20 km, relatively far from the existing regional shallow folds (Fig. 4b).Moreover, the Coulomb stress change analysis shows that EW extension is favored in the uplifted area (Figs. 5,6), which would discourage folding.In additional, Rau et al. (2022)'s strain analysis from GNSS residual deformation also revealed that a strong E-W extension (2.0-3.8 μstrain) right on Gutingken anticline.Thus, we suggest that shallow coseismic folding is unlikely a secondary source responsible for the residual uplift, whereas coseismic mud diapirism should be more plausible for the study area.
The evolution of mud diapirism is often described in geological time scale, but a case study was found in which the mud diapir was surprisingly sensitive to shortterm earthquake (Abbasi et al. 2016).In the Makran Coastal Belt in SW Pakistan, an offshore small island named Zalzala Jazeera (earthquake island in the Urdu language) emerged above sea level after the 2013 Mw 7.7 Balochistan earthquake.Farther than 380 km from the epicenter, this mud island was raised a total height of 18-21 m and was mostly composed of sheared mudstone.Similar events along the offshore Makran coast were also found back to 1945 (Delisle 2004).

Conclusions
We discover how mid-crust moderate earthquake could trigger coseismic mud diapirism through induced coseismic stress changes and produce additional amount of coseismic deformation.While the mechanism of mud diapirism in SW Taiwan is relatively unclear, we demonstrate an effective way to model diapirism through sill-like tensile dislocation and provide crucial information for advance study.Our model is also consistent with interseismic deformation pattern, which reveals the potential of crustal deformation produced by mud diapirism in different stages of earthquake cycle.Furthermore, our findings prompt a reassessment of hazards risk of mud diapirism where major infrastructures are proposed.

Fig. 1
Fig. 1 InSAR coseismic deformation fields with geological structures.a Ascending acquisitions (2-14 February 2016) and b descending acquisitions (4-16 February 2016) obtained from Sentinel-1A.Both satellite orbiting directions of ascending and descending path are denoted with white arrows aside (a and b).The look angles of each acquisition are denoted with black arrows respect to vertical line are put in the corner in a and b, for ascending path the look angle is 33° and for descending is 39°.The warmer the color, the greater the movement toward the satellite, in other word, upward motion of the ground.The thin gray rectangle indicates the main fault model of the MeiNong earthquake from Lee et al. (2016) with dark red line denoting its top edge.The small bold black rectangle is the position of the sill-like dislocation model used in this study.The red lines and black lines indicate active fault traces and axes of anticlines (Chen 2010) from the Central Geological Survey (CGS) of Taiwan, respectively.Brown shades denote the distribution of mud diapirs compiled by Chen et al. (2014).Yellow triangles indicate the location of mud volcanoes.Circles with colors indicate projected line-of-sight (LOS) displacements from continuous GPS station.The blue square denotes the announced location of the south industrial waste landfill, and the salmon square in a and b is the suspended site of the MeiNong dam, which is also close to the epicenter of the earthquake (yellow star).The straight black line denotes the location of profile A in Fig. 4. c and d Are zoom in for area of interest.Names of anticlines noted in c, TNA Tainan anticline, CCA Chungchou anticline, DGSA Dagangshan anticline, BPSA Banpingshan anticline, GTKA Gutingkeng anticline

Fig. 2
Fig. 2 Model fitting for the InSAR observations.The upper row (a-e) exhibits Sentinel-1A ascending InSAR observations and model fitting from the main fault (Lee et al. 2016) and the sub-patched sill-like dislocation model (black rectangular in d and i).The lower row (f-j) indicates the same data for the descending observations.The star is the MeiNong earthquake's epicenter.The black diamond indicates the location of the Chungliao tunnel, where interseismic uplift reaches 20 mm/yr (Ching et al. 2016).The black-lined polygon in the subfigures denotes the interseismic uplift area (37 mm/yr) observed from ALOS-2 InSAR data (Tsukahara and Takada 2018).The red-lined polygon denotes the low Vs anomaly (< − 10%) at depths of 0-2 km (Kuo-Chen et al. 2017).Yellow triangles indicate the location of mud volcanoes.Black line indicates Gutingkeng anticline (GTKA in Fig. 1)

Fig. 3
Fig. 3 Sill-like dislocation model in 1.4 km depth (this study) and main fault model (Lee et al. 2016; obtained from http:// tnem.earth.sinica.edu.tw) spacing relationship.a Is 3D view of these two models.b Top view of our sill-like dislocation.The red line indicates the top edge of the main fault.The yellow star denotes the GCMT hypocenter of the MeiNong earthquake.The red-lined polygon denotes the low Vs anomaly (< − 10%) at depths of 0-2 km (Kuo-Chen et al. 2017)

Fig. 4
Fig. 4 Schematic interpretation of the coseismic triggered mud diapirism along profile A in Fig. 1. a Projected profile of ascending InSAR data and model fitting.Descending data are presented in Additional file 1: Fig. S4.Transparent white squares with error bar denotes the continuous GPS coseismic displacement projected to LOS direction.The gray areas indicate the Pitou and Naplin anticlines from Le Béon et al. (2017) and Gutingkeng anticline(Wang et al. 1998;Chen 2013).Where sudden changes of final residuals occurred are denoted as Lungchuan, east-dipping fault and west-dipping backthrust merged at surface.b Cross-section afterMai et al. (2021).Mud diapirs (MD) are outlined with dash lines.Here we re-interpreted these diapirs with mobile shales in the center of folds, denoted as gray area.The opening region from sill model proposed by this study is denoted as light green region, its location overlapped with MD-3 diapir fromMai et al. (2021) and extend to west.The main fault plane (thick gray line) and hypocenter (red star) of 2016 Meinong earthquake located at 15 km in depth

Fig. 5 a
Fig. 5 a Coseismic dilatational strain at depth of 1.4 km.An unclamping region (positive value) matched our sill model (black rectangle).The gray rectangle with red edge is the main fault plane from Lee et al. (2016).b The locations of north-south and east-west planes above the main fault.c and d Are different view angles for static normal stress changes on these planes

Fig. 6
Fig. 6 Coseismic strain and stress change fields induced by main fault slip of 2016 MeiNong earthquake.Dilatational strain (a-c), normal stress changes on E-W striking vertical planes (d-f) and on N-S striking planes (g-i), respectively.Warm colors are expansion in dilatational strain and unclamping in normal stress change, cold colors are the opposite