Tracking surface and subsurface deformation associated with groundwater dynamics following the 2019 Mirpur earthquake

Abstract

The Mirpur Mw 5.8 earthquake on September 24, 2019, produced extensive liquefaction-induced surface deformation (LISD) in the surrounding villages. Due to the complexity of seismic hazards and the occurrence of their effects on a large spatial scale, the resulting surface, and subsurface deformation are often poorly resolved. To cover spatially extended LISD, the PSInSAR technique provided subsidence and uplift rate values ranging from −110 to +145 mm/yr consistent with the spatial distribution of the mapped liquefaction features. The most prominent surface change occurred in Abdupur and Sang villages. GPR measurements were conducted to map the near-surface cracks produced by transported liquified sand into the shallow subsurface layers and other liquefaction features (elevated groundwater table, conductive clay pockets, fractures, sand dikes, and water-enriched zones). Thus, the GPR survey assisted in the reconstruction of these structural and hydrogeological features on the near surface. In addition, the highly vulnerable zones were identified and mapped using space- and ground-based remote sensing measurements supported by the field observations. The results highlight the effectiveness of the proposed novel approach for detailed assessment of the coseismic liquefaction-induced deformation on- and near-ground surfaces by identifying areas prone to failure during earthquakes and thereby can help with hazard mitigation.

Keywords:

• Earthquake
• coseismic effects
• liquefaction
• ground penetrating radar (GPR)
• geological survey

Introduction

The Mirpur area is located in the floodplain of the Jhelum River. Alluvial systems have been recognized worldwide for several decades as liquefaction-prone settings (Khan et al. 2021a). During a seismic event, a granular deposit transforms from a solid into a liquefied state due to the increased pore-water pressure. When subjected to seismic shaking, the water-saturated sandy soil behaves as a viscous liquid when it loses its shear strength, and the phenomenon is called liquefaction (Youd 1977; Jefferies and Been 2015). Earthquakes are considered one of the most important geohazards causing liquefaction, and surface manifestations in the form of uplift, slope failure, subsidence, lateral spreading, ground fracturing, and sand blows (Youd 1995). The coseismic liquefaction hazard results in extensive damage to residential buildings, tens of meters long ground fractures, disruption of lifeline facilities, and partial to complete collapse of transportation infrastructure (e.g. bridges and roads) (Galli et al. 2012; Orense et al. 2012). In this context, the Mirpur earthquake is not an exception and documents the latest example in Pakistan, leaving its footprints in the flood plain of the Jhelum River (Khan et al. 2021b). In the event of such a disaster, irrespective of the transmitter’s location and receiver’s location, remote sensing is considered a valuable source of information for mapping earthquake-induced deformation compared to traditional destructive techniques (e.g. drilling and trenches). Broadly, two types of radar remote sensing tools (ground penetrating radar (GPR) and Permanent Scatterer Interferometric Synthetic Aperture Radar (PSInSAR)) have proven their potential in mapping earthquake-induced deformation around the world.

Within the domain of natural hazards, the use of space-borne remote sensing has become popular and efficient (Muhammad and Tian 2020; Muhammad et al. 2021). The evolution of geospatial technologies and advanced data processing has the ability to provide improved quality data (Muhammad and Thapa, 2020; 2021) with repeated spatiotemporal coverage covering extensive areas in rough geomorphological and geological conditions (Gul et al., 2020; Hassan et al. 2021). The post-disaster remote sensing image visualization helps in emergency interventions and decision-making processes during the disaster (Tian et al. 2017). Coseismic effects such as ground fractures of around 100 meters in length accompanied by extensive lateral spreading and thousands of landslides were observed in the 2019 Mirpur and 2005 Kashmir earthquakes, respectively (Owen et al. 2008; Khan et al. 2021a). In shallow aquifer systems, changes in surface water bodies, groundwater levels and geochemical variations following an earthquake have been studied in Asia (Lapenna et al. 2004; He et al. 2017; He and Singh 2019; Jing et al. 2022). The Mirpur earthquake produced extensive liquefaction effects in the surrounding villages of Mirpur city. Seismic liquefaction is a cascade type of coseismic hazard that carries implications for agriculture, water resource management, the safety of residential, and key infrastructure. Despite its significance for various sectors relevant to food security, societal health, and public safety, it is difficult to identify and precisely map coseismic deformation mainly due to its large spatial extent. The detection and adequate characterization of liquefaction-induced features require the application of surface and subsurface techniques. Within this context, ground-based remote sensing techniques, such as GPR is one of the non-destructive geophysical methods that provide high-resolution subsurface images to detect structural, stratigraphic, and hydrogeological features. However, GPR cannot provide all the required information for hazard mitigation due to certain limitations (e.g. less spatial coverage, rough topography, labour-intensive, time-consuming, and high cost compared to the space-borne imagery) even on a small-scale like Mirpur district. To meet the above requirements, freely available space-borne data sets have successfully been tested as possible candidates for the assessment of surface coseismic effects in seismically active regions (Chini et al. 2008, 2011). Remote sensing data from a number of satellites and sensors are commonly used for hazard mapping. For instance, Landsat satellite imagery has been used to study the dynamic relationship between observed seismicity and lineament density (Nath et al. 2021). The Synthetic Aperture Radar (SAR) is generally accepted as the best candidate for mapping coseismic deformation. Interferometric SAR (InSAR) is one of the powerful remote sensing techniques of the SAR family, used to detect surface deformations over large areas with high accuracy (Kiseleva et al. 2014). InSAR-based methods allow low-speed surface deformations to be detected over vast areas with centimeter to millimeter precision (Bianchini et al. 2015; Sara et al. 2015; Huang Lin et al. 2019). From the 1980s (Zebker and Goldstein 1986; Gabriel et al. 1989) to date (Lanari et al. 2007; De Luca et al. 2017; Zinno et al. 2018; Boixart et al. 2020; Poreh and Pirasteh 2020), several studies have successfully employed the Differential SAR Interferometry (DInSAR) technique to detect surface deformations. However, several factors (atmospheric and orbital errors, spatial and temporal correlation distortion) limit the applicability of the conventional DInSAR method in studies related to surface deformations (Muhammad et al. 2021). The permanent Scatterers InSAR (PSInSAR) method (Ferretti et al. 2001), using a series of SAR images recorded at different dates, was developed to overcome these limitations (Abdikan et al. 2014; Yazici and Tunc Gormus 2020). The PSInSAR is generally accepted as the best candidate for mapping coseismic deformation in the field conditions as well as movements of urban infrastructure (Wang et al. 2017; Fiorentini et al. 2020, Saralioglu 2022). The main difference in PSInSAR-based methods arises in the mathematical approaches to obtaining PSs. Since pixels with ‘stable behavior’ are used when determining PSs, they are called Persistent or Permanent Scatter (PS). Hooper et al. (2010) developed the StaMPS/MTI (Stanford Method for Persistent Scatterers/Multi-Temporal InSAR) method by considering the PS interferometry technique with a different approach. In this method, as in PSI, more stable PS points are determined by using both amplitude and phase instead of only amplitude in the selection of PS points. Thus, PS points used to represent earth features in the StaMPS/MTI method are resistant to different atmospheric and climatic conditions and small changes in satellite angle (Hooper et al. 2010). In addition, although PS points usually occur in structures such as concrete, metal structures, and rock fragments, they can also occur in other land types. Due to this advantage, the determination of surface deformations in this study was carried out using the StaMPS method. On the contrary, space-borne data cannot provide complete information on the near-surface features produced as a consequence of an earthquake with the potential to damage the built environment severely. Hence, our proposed methodology utilize the state-of-art GPR technology with the PSInSAR method to reduce interpretational ambiguities common in a single technique and develop an improved understanding of the coseismic liquefaction hazard.

Like other non-destructive near surface geophysical methods (e.g. ERT), the ground penetrating radar (GPR) method has been applied to shallow subsurface investigations due to its high-resolution, time and cost-effective nature (Lapenna et al. 2003; Khan et al. 2019, 2020, 2021c). Recently, GPR has gained popularity in studies related to the detection of faults and fracture networks (Schwarz and Krawczyk 2020), slope instabilities (Khan et al. 2021c), and landslides (Hu and Shan 2016). The geophysical methods, for instance, GPR is one of the reliably accurate mapping tools to study a single site and/or imaging of a localized subsurface deformation but difficult to perform such surveys over an extensive earthquake/Karst depressions/landslide affected area to detect the near-surface target features. Among the aforementioned geohazards, few studies focusing on coseismic liquefaction and related ground failure have been conducted using field GPR measurements (Liu and Li 2001; Baradello and Accaino 2016).

To the best of our knowledge, GPR technology has not been tested to image the coseismic features as a consequence of the soil liquefaction phenomenon through combined results derived from the traditional geological survey in the field, PSInSAR, and high-resolution GPR measurements. Thus, our work presents a novel approach to study coseismic liquefaction effects (surface and subsurface deformation, hydrogeological changes) following the Mirpur earthquake.

On September 24, 2019, an earthquake of Mw 5.8 with ∼10 km of depth struck Mirpur city, causing dozens of fatalities and leaving scores of people injured with severe damage to civil infrastructure (www.ndma.gov.pk). The maximum ground shake intensity of VI resulted from the 2019 Mirpur earthquake (Figure 1b) (https://earthquake.usgs.gov/). The induced liquefaction and long fractures (>100 m) damaged all types of infrastructure, including buildings, walls, bridges, agricultural fields, water supply networks, and roads (Figure 2). In addition, the coseismic effects are observed in the nearby villages (Abdupur, Sang, Manda) located within 6 km from the epicenter (Figures 1–3).

Figure 1. (a) Tectonic map showing study area (modified after Jayangondaperumal and Thakur 2008). The red, black, and blue stars on a, b, and c, respectively, denote the epicenter of mainshock associated with the earthquake and its focal plane mechanism solution (data source: https://www.usgs.gov/, last accessed: 30-09-2019), heavy black line denotes the faults in the region, blue diamond denotes major cities in the vicinity of epicenter; (b) Macroseismic intensity map; (c) Liquefaction probability map; (d) surficial geological map of the study area showing locations of the ground deformation (after Arif et al. 2000; Khan et al. 2021a).

Figure 2. Field photographs showing ground deformation at various locations (1-8, see Figure 8 for their locations). surface fractures in the field (1-3,6); damaged infrastructure (4,8); small-scale landslide (5); and ground subsidence/uplift (7).

Figure 3. Map showing location of the study sites of field observations where GPR measurements were taken. Inset shows location of the study area at the regional scale. The rose diagram and trends (line tracing on the map) of the fractures at representative sites (S1, S3, S4, and S7), and sand blows as green ellipses at S4 are also shown.

Hydrogeological and seismological setting

In the northwestern Himalayas, the Indian plate’s northward motion at 12–14 mm/yr with respect to the Eurasian plate makes the Himalayas one of the most seismically active regions across the globe (Sreejith et al. 2021 and references therein). This collision process has produced devastating earthquakes in the whole region, including Pakistan, for instance, the 2005 Kashmir earthquake of Mw 7.6, which is the deadliest one in human history, killing more than 70 thousand people and triggering several thousand landslides (Owen et al. 2008, Yousuf et al. 2020). There are several thrust and strike-slip faults in and around the study area. Among the major seismic events associated with the fault system, very few ruptured the surface of the Himalayan Frontal Thrust (HFT) (Mencin et al. 2016; Bilham 2019). The 2019 Mirpur earthquake occurred in a region characterized by faults of different natures, such as Salt Range Thrust (SRT), Himalayan Frontal Thrust (HFT), Domeli thrust, and Jhelum strike-slip fault (Figure 1a).

The 2019 Mirpur earthquake caused 39 fatalities, injured more than 700, and caused damage to civil infrastructure (NDMA, 2019; USGS, 2019). It occurred in the Jhelum flood plain resulting in a maximum ground shake intensity of VI (Figure 1b,d). Despite a seismic event of small magnitude, the Mirpur earthquake caused severe surface deformation patterns which are comparable to seismic events of larger magnitude (Lunina and Gladkov 2016; Rodríguez-Pascua et al. 2016), including the Mw 7.6 Kashmir earthquake in the same region (Jayangondaperumal and Thakur 2008). The local hydrological setting and anthropogenic features such as the Jhelum River, Mangla dam, Upper Jhelum Canal (UJC), rice fields, and shallow aquifer system) may have significantly contributed to the extensive liquefaction effects. Moreover, the geology is dominated by the loose alluvial deposits overlying a liquefied sandy soil that also played a key role during the seismic hazard. The liquefaction-induced fractures of more than 100 m in length affected the diverse nature of infrastructure, including mosques, residential buildings, schools, the UJC, and bridges. The ground failures produced severe damage manifestations which were observed on the surface and within the near-surface in the epicentral area, including the most common ones within 50 m from the UJC (Figures 1b–d, location 4 on Figure 2).

Materials and methods

To study surface and subsurface deformation features resulted from the Mirpur earthquake (Figure 4), this work presents a novel methodology using space-borne and ground-based radar techniques and field geological survey (Figure 5).

Figure 4. (a, b) Schematic representation of before and after the Mirpur earthquake liquefaction hazard and related coseismic surface and subsurface deformation features.

Figure 5. A flowchart of ground and space-borne radar-based approach for coseismic liquefaction.

Space-borne RADAR/Permanent Scatterer Interferometric Synthetic Aperture Radar (PSInSAR)

Sentinel-1 SAR images were used in the space-based determination of surface deformations. The Sentinel-1 carries a C-band SAR sensor that enables the acquisition of data regardless of weather conditions in day and night. Sentinel-1 consists of two satellites including Sentinel 1 A and Sentinel 1B. The temporal resolution of each satellite is 12 days, the repeat cycle can be reduced to 6 days using the data from both satellites. Currently, Sentinel 1 B satellite data is not available due to an anomaly in the satellite’s power system. This study used 30 images from Sentinel 1 A satellite in descending mode between April 02, 2019 and March 27, 2020. The Single Look Complex (SLC) required for interferometric analyses was selected as the data type, and the Interferometric Wide Swath Mode (IW) was selected as the Sensor mode. An area of 250 km is covered by the data obtained in IW mode. Images of the Sentinel program run by the Copernicus program of the European Space Agency (ESA) can be accessed freely through the Sentinel Scientific Data Hub (https://scihub.copernicus.eu/dhus/#/home). Processing of images was carried out with SNAP-StaMPS integration.

PSInSAR Preprocessing

A total of 30 sentinel-1 images were acquired for this study. Using the TOPSAR Split command, the relevant one from the IW1, IW2, and IW3 modes falling into the working area was selected. Precise Orbit Determination (POD) information is then downloaded for each image. Based on a single main image, a total of N + 1 SAR images are used to create N interferograms in the PSI technique. Such interferograms carry the phase differences between the dependent images and the selected master image. For the selection of the main image, priority is given to the temporal baseline length. It is ensured that the main image is in the middle of the other dependent images at the temporal base length. In this study, the image dated August 24, 2019, was chosen as the master image. Single master-slave image stacks are created in the S-1 Back Geocoding stage. TOPSAR-Deburst operator is used to joining all burst data into a single image. After this stage, the study area was subsetted. The subtract topographic phase function is activated during the interferogram phase. In addition, 1 arc second SRTM (Shuttle Radar Topography Mission) data was used for topographic phase effect removal. At this stage, the function that enable the creation of an elevation band and orthorectified latitude and longitude coordinates has been activated. The generated data was saved with the StaMPS export function. In the next step, the converted data are imported into StaMPS using a special script called mt_prep_snap. In this part, the amplitude dispersion was set to a threshold value of 0.42, which minimizes random amplitude variability and eliminates highly unrelated pixels in some areas covered with vegetation, farmland, or snow. The step-wise procedure for StaMPS/MTI 4.1b1 was followed as reported by Hooper et al. (2018).

Estimate phase noise and correction

The second phase includes the estimation of Reflective point candidates and estimated phase noises (γ) based on various preliminary data. This is an iterative step that estimates the phase noise value for each candidate pixel in each interferogram. Filtered_grid_size = 50 m was selected from the parameters controlling this process.

PS points are selected according to noise characteristics. This step also estimates the percentage of random (non-PS) pixels in a scene from which density per km² can be achieved. If the noisy pixels are more than the estimated threshold limit, the PS is subtracted from the set point. It was selected as select_method='DENSITY' from PS selection parameters. The very noisy ones from the selected PS points are cleaned.

Selected pixels are verified for mismatched viewpoints and correlation errors due to DEM are calculated and debugged. Phase analysis is performed with two different analysis methods called 3D and 3D_QUICK. In this study, 3D was chosen as the unwrap method.

The spatially related viewpoint error (SCLA) and atmospheric trajectory errors (AOE) are calculated. For more details on image processing and the stages of obtaining PS points, please refer to Saralioglu (2022).

PSInSAR survey results

After analyzing the InSAR data of the liquefaction-affected area, we have identified the severely damaged zones on the ground that were further validated through obvious surface manifestations during the geological field campaign. Thus, space-borne measurements helped us in locating important sites undergone severe ground failure that subsequently allowed us to effectively plan ground-based geophysical radar surveys to detect the near-surface liquefaction structures (Figure 6).

Figure 6. The velocity (mm/year) of PS points in the whole study area. The Blue and Red regions indicate locations with anomalous velocity values.

Figure 6 shows the PS points produced by the StaMPS method in the region within the study area and the annual displacement rates of these points. In total, 261544 (Figure 7) PS points were created. The displacements of the PS points in the LOS direction range from −110mm to +145mm. It is important to examine the PS points produced in Figure 6 for accuracy. In Figure 6, the red dots indicate subsidence, and the blue dots show the uplift. The fact that the red and blue points with high-velocity variation differences are not adjacent to each other and the clustering of the point groups belonging to the close velocity values in certain places shows that the PS points are compatible in terms of accuracy (Figures 6, 8, and 9). The baseline distance between the images should be short of achieving this harmony. In addition, if a point with a different velocity change value is detected between neighboring points during detailed examinations, this point is not used in drawing the time graph (Figure 10). The distributions of all PS points generated are shown in Figure 7. Most of the spots display velocity values between −14 and 18 mm/year. These points are shown in light yellow in Figure 6.

Figure 7. Distribution of PS points.

Figure 8. The velocity of PS points at selected locations covering the sites (S1-S8) where field observations were made.

Figure 9. The red dots represent the PS points at locations 1-11 representing high deformation in the study area. Time-series plots of these points are shown in Figure 10 (For locations, see Figure 8).

Figure 10. Selected points with high deformation rates (The date with the red box indicates the date of the earthquake) in all plots.

Other PS points show higher annual velocity changes. There are 53177 (Figure 7) points with −46 to −110 mm/year change rates. These points are shown in Figure 6 with the colors red and orange. The (-) direction change in this region means that the points move away from the RADAR sensor. The + directional PS point count is 35715 (Figure 7), with high change rates (49 to 145 mm/year)). The + direction change means that the points in this region get closer to the RADAR sensor. In short, PS points with positive signs indicate surface deformation movement upwards and eastwards, while negative values indicate surface deformation towards the bottom and west (Tofani et al. 2013). In general, Figure 6 shows us that most of the study area is subjected to high surface deformation. The red, orange, and blue dots show the places where the deformation is most severe.

Figure 10 shows the time-series graphs of the 10 points shown in Figure 9. The region in the red box in the figures includes the two images obtained before the earthquake (September 17, 2019) and after the earthquake (September 29, 2019) and is closest to the earthquake date. In Figure 10, the effect of the earthquake can be observed in selected places within the study area. For example, in Figure 10a, it is seen that the time series graph, which continues steadily, jumps with a difference of approximately 50 mm between the two pre and post earthquake images. The time-series graphs indicate an uplift in the post-earthquake scenario. In some of these points (2, 4, 6, 8, 9, 11). This amount of increase is relatively high. The anomaly revealed by the earthquake become evident in the time-series graphs.

The satellite imagery allowed us to identify stable areas with no to moderate deformation, and the ones showing extensive unstable areas with ‘intense deformation’ (Figure 8). Spaceborne radar remote sensing results also helped us in prioritizing and selecting the sites for ground-based radar remote sensing (GPR) measurements. A GPR survey was carried out for further investigation.

GPR survey

Ground Penetrating Radar (GPR) is an effective geophysical technique used for non-invasive investigations (e.g. soil profiles and fractures/cracks mapping) in the shallow subsurface (Piscitelli et al. 2007; Loperte et al. 2011; Khan et al. 2019, 2021c). A GPR survey was conducted to map subsurface liquefaction structures and detect hydrogeological features following the Mirpur earthquake. The GPR results were aimed to verify the PSInSAR results by mapping the structural features and estimating how deep these features can penetrate in the near surface during the liquefaction hazard.

To confirm the subsurface deformation, and reconstruct the subsurface structural and lithological setup, GPR data were collected on both sides of the UJC. We conducted field GPR experiments using different frequencies (250, 500, and 800 MHz) along ten survey profiles at four sites S1, S3, S4, and S7 in the vicinity of UJC (Figure 3), in the 2019 Mirpur earthquake-affected area. Literature (Allen 1986; Ambraseys 1988, Ko et al. 2017) shows that an earthquake of magnitude (Mw 5–6), for example, 5.8 Mw Mirpur seismic event, could induce coseismic liquefaction deformation within 5 m depth (Ambraseys 1988; Obermeier 1996; Owen and Moretti 2011; Zhao 2012; Shao et al. 2020, Khan et al. 2021a; and the references therein). Therefore, our geophysical results are primarily based on the 500 MHz antenna system due to its suitability in terms of resolution and its depth detection capability to reasonably resolve the most common liquefaction features in the target depth interval (0–6 m) (https://www.guidelinegeo.com).

Additionally, we employed 250 and 800 MHz at a few sites to map target features located at very shallow and relatively larger depths, respectively. For an improved understanding of the coseismic deformation following the Mirpur earthquake, GPR data were collected on both sides of the UJC to reconstruct the subsurface structural and lithological setup. The radar measurements were subjected to standard processing steps such as time-zero correction, band-pass filtering scheme, and gain control using RadExplorer software to extract important information about the targeted features in the context of coseismic liquefaction (Figures 11 and 12).

Figure 11. GPR profiles showing different coseismic structures in the Mirpur earthquake affected area. Locations of the GPR profiles are (a) along the bank of UJC (location 4 on Figure 8); (b) at site S4; (c) at site S3; (d) at site S7; and (e) at site S1; and (f) at site S4.

Figure 12. GPR profiles showing different coseismic structures in the Mirpur earthquake affected area. Locations of the GPR profiles are (a) at site S4; (b) at site S1; (c) along the bank of UJC (location 4 in Figure 8); and (d) at about 100 m distance from site S7.

The quality and penetration depth of the radar sections was highly variable due to the high moisture content in the irrigated area and the clay content of the uppermost soil. For instance, we were unable to track any significant radar signature on the GPR profile (Figure 12d) in the context of liquefaction. On the other hand, some 500 MHz GPR profiles displayed important coseismic features such as Figures 11b and 12a showed the radar signatures that maybe related to liquified sand. This documents one of the diagnostic structures for the liquefaction occurrence following the seismic event. The presence of the liquified upwelled material/sand dike is also supported by the sand blow on the ground surface as field evidence from the same area (Figures 3 and 4b).

Figures (11a and 12c) show inclined reflectors and deformed strata based on 500 and 250 MHz antenna systems, respectively, that may probably be associated with the subsidence phenomenon as a consequence of the earthquake. Similar deformational patterns related to ground subsidence in this area could also be seen from the space-borne radar measurements (Figures 6 and 8). Field observations (Figure 2) attest to the radar-derived results in the study area. The Mirpur earthquake-induced deformation in fractured near surface can be seen in Figures 11c,e, and 11d, where typical hyperbolic signatures represent the existence of fractures, and some movement along one of such fractures can be observed from the high resolution radargrams of 250 and 800 MHz, respectively. The disrupted strata and movement along fracture observed on GPR profiles further highlight the deformed character of the near-surface as a consequence of coseismic hydrogeological changes (Figures 11e and 12b).

The seismically perturbed groundwater table resulted in the fractured shallow subsurface aquifer system (Figures 11f and 12a). The elevated groundwater table shown by a strong reflector on 500 MHz radargram in the near-surface contributed as one of the main factors in producing widespread liquefaction hazard. The ground-based radar measurements showed different deformation patterns of coseismic liquefaction, consistent with the satellite-based radar results at some locations and surface manifestations observed in the subsequent phase of the field investigation.

Survey of ground deformation

To confirm the results of space-borne radar data analysis (i.e. PSInSAR) and to get a picture of the subsurface deformation using ground-based geophysical radar technique (i.e. GPR), a detailed field campaign was undertaken at sites showing sufficient surface deformation on InSAR results (Figure 6). Ground deformation and infrastructure damages were observed across the epicentral region (Figures 1–3). Notable surface deformation was confined to a 50 m wide zone along the banks of UJC and at eight different locations (S1-S8) across the floodplain deposits of the Jhelum River in the Mirpur area downstream of the Mangla reservoir (Figures 1 and 8). Sand blows/dikes were found at three sites (S2, S4, and S6), whereas surface fractures were documented at all eight sites (S1-S8). The flow pattern of the extruded sediments was symmetric along the fractures, whereas the geometry of the sand blows was elliptical (Figure 3g of Khan et al. (2021a)).

Fractures running along the canal banks were 0.5–160 cm wide and 6–187 m long and were characterized by vertical and lateral displacement and subsidence of coherent sediment masses toward the canal (location 4 on Figures 2 and 8; Khan et al. 2021a). These fractures are oriented parallel to the UJC (i.e. ∼ N60°W ± 15°). Predominantly, fractures found at a distance of >50 m from the canal bank (i.e. at sites S1-S8) displayed width and length of 0.5–50 cm and 4–134 m, respectively. However, at site S2, the fracture width, and length reached 180 cm and 280 m, respectively. Most of the fractures displayed lateral extension, with some showing vertical movement of 2 cm to 135 cm (locations 2 & 7 in both Figures 2 and 8). Modest left lateral strike-slip movement (0–15 cm) was also recorded along a few fractures at sites S1, S5, and S7 (Figure 3f of Khan et al. (2021a)). All the fractures are striking towards N25°W–N45°E.

Discussion

The Mirpur earthquake produced liquefaction-induced ground deformation in the surrounding villages of the epicentral region covering a large area (Figures 1–3). Satellite-based imagery helped us in analyzing spatially distributed deformational patterns in the epicentral area (Figures 6 and 8). Thus, PSInSAR measurements made possible the tracking of surface deformation at a large scale and assisted in prioritizing the sites for detailed ground-based investigations (geological and geophysical). Substantial deformations are observed in the northern and central parts of the study area. In the northern part, annual velocity variations have reached as high as −115 mm. In the middle part, annual velocity changes have reached +145 mm.

Significant ground deformation was observed at selected sites (S1-S8) and along the banks of UJC. It included the development of surface fractures, lateral spreading, damage to the infrastructure, sand blows, and ground subsidence. Generally, an increase in the amount of deformation is observed with decreasing distance from the UJC. It shows that UJC plays an important role in the coseismic liquefaction-induced ground deformation, which may have increased the water content of the unconsolidated sediments in its immediate surroundings. Also, the ground fractures observed at eight stations show no relationship with the northeast dipping causative fault (Tan et al. 2022). All this ground deformation is confined to the broad zone showing positive velocity values on PSInSAR results interpreted as the areas which have experienced uplift (Figure 6) and high to very high ground deformation (Figure 8). The localization of the ground deformation within the unconsolidated sediments is attributed to the uplifting. The same is also reported for the other earthquake events in the surrounding region, for instance, the 2005 Kashmir earthquake (Sayab and Khan 2010).

Coseismic groundwater dynamics fractured the shallow aquifer system in the study area. GPR sections showed the subsurface fractures disrupting the geoelectric layers with varied thicknesses and dielectric characteristics within the upper 5 m. The radar-stratigraphic section confirms the presence of an elevated groundwater table and water-filled fractures in the near surface. Other important liquefaction structures included sand dike, inclined reflectors, subsidence, and movement along the disrupted reflectors due to ground shaking and enhanced pore pressure during the liquefaction phenomenon. In the context of hydrogeological hazards, similar abnormal pore pressure is observed at larger depths within the subsurface formations/petroleum reservoirs. Tectonic stresses, among other factors, trigger such pressured zones that may result in catastrophic events (e.g. well blowouts) that directly threaten safe drilling operations (Khan et al. 2022).

The mapped fractures facilitated groundwater movement toward the surface and transported the liquified sandy material resulting in the sand plugs in shallow subsurface. These processes contributed to soil moisture variability revealed by geophysical results, thereby emphasizing the importance of improved understanding and monitoring of hydrologic and agricultural systems in tectonically complex and seismically active areas. The reason for such monitoring is critical in seismically active areas where irrigated lands are more vulnerable to such types of coseismic hazards and hence may directly affect agricultural food production. This may subsequently compromise the food security and economic situation in a densely populated country like Pakistan where the agriculture sector is considered a major source of food and contributes to its economy (Kishida et al. 2009).

Conclusions

We proposed a new methodological approach based on space-borne and ground-based radar techniques (GPR and the PSInSAR) to assess coseismic liquefaction-induced surface and subsurface deformation.

The ground displacement with a high spatial and temporal resolution in the epicentral area was derived using PSInSAR. The PSInSAR-derived deformation data allowed the identification of uplift and subsidence areas that would otherwise have been overlooked by conventional geological mapping due to their subtle geomorphic signatures or their buried character. These maps delineated the deformation zones on the ground that are triggered by the earthquake-induced liquefaction phenomenon. The results from the interferometric analysis were instrumental in the identification of the locations with severe deformation. The annual changes in velocity values highlight high deformation hotspots in Abdupur and Sang villages. Our result shows that key structures such as the Mangla dam should be made more sensitive to earthquakes in the future, and existing structures should be re-evaluated. The PSInSAR also assisted in selecting the potential sites for GPR experiments.

GPR results confirmed the presence of the most common coseismic liquefaction structures, such as sand dikes, fractures, and water-filled cavities in the uppermost 5 m of the subsurface. A very promising large hyperbolic shape-like structure at 1.6 m depth attests to well-developed fractures of considerable dimension (ranging between 0.5 and 180 cm wide and 4 and 280 m long) that may potentially damage the common residential structure. GPR sections reconstructed the shallow subsoil stratigraphy and captured the hydrogeological features (the rise in the groundwater table). Parallel to this, the PSInSAR analysis indicates that this section appears to be affected by very high to high deformation closely located in the Abdupur and Sang villages.

The geological field survey was conducted to map and confirm the surface manifestations of both space-borne (PSInSAR) and ground-based (GPR) radar techniques and map liquefaction structures, providing an indirect clue of the occurrence of the seismic hazard at shallow depths. All the ground deformation is restricted to the floodplain deposits of the Jhelum River highlighting the control of lithology (i.e. loose floodplain sediments) and moisture content (i.e. shallow water table) over the ground deformation in case of an intermediate magnitude seismic event (i.e. 2019 Mirpur Earthquake). For detailed documentation of the liquefaction phenomena and associated ground deformation, trenching at selected sites will provide further insights into the liquefaction hazard.

Acknowledgments

The authors thank the National Centre of Excellence in Geology (NCEG), University of Peshawar, Pakistan, and GIS and Space Applications in Geosciences (GSAG) Lab, National Centre of GIS and Space Applications (NCGSA) for providing research facilities and the Higher Education Commission (HEC), Pakistan, for financial support to conduct the research work. This work acknowledges the support of the International Centre for Integrated Mountain Development (ICIMOD), funded by core funds of the ICIMOD and contributed by the governments of Afghanistan, Australia, Austria, Bangladesh, Bhutan, China, India, Myanmar, Nepal, Norway, Pakistan, Sweden, and Switzerland. The views and interpretations in this publication are those of the authors, and they are not necessarily attributable to their organizations.

Disclosure statement

The authors have no potential conflict of interest.

Data availability statement

The data in this study are available from the first and corresponding authors upon reasonable request.

Additional information

Funding

This study was supported by Higher Education Commision, Pakistan.