Interplate rheological contrast revealed by asymmetric deformation after the 2023 Kahramanmaraş earthquakes

Domain: chemistry
Published: 2026-02-25
Editorial outcome: major_revision

Abstract
Postseismic surface deformation offers key insights into subsurface lithospheric rheology and earthquake fault behavior, yet geodetic data often fail to capture its full complexity. Here, we use satellite radar images to derive time series of three-dimensional surface displacements following the 2023 Kahramanmaraş (Türkiye) earthquakes. The results show strong spatial and temporal deformation asymmetry across the East Anatolian Fault with the stiffer Arabian lithosphere exhibiting slower deformation decay than the Anatolian side. The temporal asymmetry excludes elastic afterslip from being the main postseismic process and instead the postseismic response is dominated by asymmetric viscoelastic relaxation as well as poroelastic rebound. These findings underscore the importance of both temporal and three-dimensional spatial deformation data in advancing understanding of post-earthquake recovery processes and subsurface lithospheric structure.

Introduction
Postseismic processes following major earthquakes, including viscoelastic relaxation, poroelastic rebound, and afterslip, are of paramount importance for understanding lithospheric rheology, fault slip characteristics, and regional tectonic movements1–3. One of the most direct manifestations of these postseismic processes is surface deformation, which can be measured using geodetic techniques such as interferometric synthetic aperture radar (SAR, InSAR) and global navigation satellite system (GNSS). While different postseismic processes reflect the properties of the subsurface medium from distinct perspectives, the surface deformation measured results from a combination of all these different postseismic processes. Consequently, significant efforts have been made to distinguish between these processes. For instance, by leveraging prior information or simplified assumptions, studies often focus on a specific postseismic process while neglecting other secondary processes4–6. Alternatively, based on the spatio-temporal characteristics of the different postseismic processes, attempts have been made to sequentially separate these processes from the observed surface deformation7–9. However, accurately distinguishing between different postseismic processes remains highly challenging due to the often incomplete geodetic coverage of the affected area or the inability to reveal the full-dimensional (temporal and 3D spatial) surface deformation. This is particularly true for the highly complex 2023 Kahramanmaraş earthquakes10,11.

The Kahramanmaraş earthquakes sequence on February 6, 2023, comprised two major events with magnitude 7.8 and 7.6, occurring at the boundary between the Arabian and Anatolian plates, i.e., on the East Anatolian fault (EAF) and Sürgü fault, respectively12. These two earthquakes were separated by only 9h in time and approximately 90km in distance, with a total fault rupture length extending 500km (Fig. 1). The seismic activity released the accumulated stress along EAF due to the westward extrusion of the Anatolian Plate13,14. The long-term tectonic evolution of Anatolian Plate has been interpreted by the westward tectonic escape driven by northward movement of Arabian Plate15,16, with the EAF acting as a key transform boundary accommodating this motion17–20. Later geodynamic models further suggested that this escape is sustained by large-scale mantle flow coupling westward slab pull in the Aegean Sea21–23. Across the EAF, the Arabian lithosphere collides into Anatolia at the EAF, and is thought to be relatively rigid and homogeneous24. Conversely, the Anatolian Plate has formed through multiple subduction and collision events, and is still deforming under a variety of tectonic processes from east to west25,26. Additionally, the pronounced lateral heterogeneity across the EAF is highlighted by the significant differences in the thickness of subsurface structural layers27–29 in the region’s seismic structure. These complex geological settings and the associated seismogenic faults suggest that the postseismic deformation pattern is highly intricate.

In this study, we combined nearly two years of postseismic Sentinel-1 images with time-series InSAR30,31 and burst overlap interferometry (BOI)32,33 methods to capture the full-dimensional (temporal, east-west, north-south, and vertical) deformation following the earthquakes. The availability of such comprehensive surface deformation data enables discrimination of the distinct postseismic processes, i.e., viscoelastic relaxation in the lower crust and upper mantle, poroelastic rebound of the crust, and afterslip on and below the coseismic fault slip. This capacity stems from the dissimilar contributions of these mechanisms to different components of the surface deformation34, as well as from the notable temporal and spatial deformation asymmetry observed in this case.

Results
Using time-series InSAR and BOI methods, we obtained four relatively independent postseismic displacement observations with distinct geometric differences from both ascending and descending orbits in the line-of-sight (LOS) and azimuth directions (see “Methods”). To ensure complete coverage of the deformation area, both ascending and descending datasets include three adjacent orbital tracks (each with a swath width of 250km), with each track consisting of 4–5 frames (Supplementary Fig. 1). The InSAR LOS observations (Supplementary Fig. 2) provide spatially continuous displacement maps in the east-west (EW) and vertical directions. The BOI observations, on the other hand, are only available in burst overlap areas (Supplementary Fig. 2), and the resulting north-south (NS) displacement is therefore not entirely spatially continuous. Nevertheless, BOI complements InSAR by revealing the overall characteristics of postseismic deformation in the NS direction35,36. To remove secular interseismic tectonic signals from the postseismic displacement field, we also processed Sentinel-1 data from 2015 to 2023 using both time-series InSAR and BOI methods (Supplementary Fig. 3). Comparison with publicly available GNSS solutions37,38 exhibits root-mean-square (RMS) differences of 17, 11, and 11mm between the GNSS and postseismic InSAR cumulative displacements in the EW, NS, and vertical directions, respectively (Supplementary Fig. 4).

The postseismic three-dimensional (3D) displacements (Fig. 2a, c, e) indicate that the deformation characteristics align well with the expected postseismic response from the left-lateral strike-slip earthquakes. Horizontal deformation dominates the postseismic surface displacement, reaching magnitudes of up to 15cm within two years. Regarding the near-field horizontal deformation, most areas exhibit a smooth cross-fault transition (except at the northeastern ends of the two primary fault ruptures), suggesting a limited contribution from shallow afterslip. In terms of vertical deformation, a notable characteristic is the presence of extensive scale (~200km) uplift and subsidence (~5cm) respectively north and south of the EAF. Compared to the coseismic deformation39, the spatial extent of the postseismic deformation is significantly larger, suggesting a substantial contribution from processes deeper than the coseismic fault slip40.

A pronounced spatiotemporal asymmetry is observed between the Arabian and Anatolian plates in the postseismic EW displacement field (Fig. 2). The mechanically weaker Anatolian plate exhibits higher and more spatially concentrated deformation, whereas on the Arabian side lower deformation is distributed over a broader area. To fit the temporal evolution of EW deformation, we applied a logarithmic function (i.e., d=a log 1+t/T )41,42. The results reveal that the decay time T (Fig. 2a, j) on the Arabian plate is systematically larger than that on the Anatolian side, implying that, compared to the Anatolian Plate, the postseismic deformation of the Arabian Plate originates from a more prolonged subsurface process.

Postseismic deformation that is spatially asymmetric across a fault could potentially be caused by elastic afterslip if the lithospheric structure of the two fault sides is different. However, the observed temporal asymmetry cannot be explained by elastic afterslip, as such afterslip would yield the same deformation decay times on both sides. For that reason, we model the postseismic deformation as a combination of viscoelastic relaxation and poroelastic rebound using the time series of 3D postseismic displacements. For the viscoelastic modeling, both EW and NS displacements are used. Given that near-field deformation could be influenced by shallow fault slip, we conducted an independent test by searching for the optimal model parameters using only data beyond 50km (approximately twice the seismogenic depth43–45). The results indicate that the inclusion of near-field observations has an insignificant impact on the optimal model solution (Supplementary Figs. 6 and 7). In addition, observed vertical displacements were excluded from the viscoelastic modeling process, as they primarily stem from poroelastic rebound, a mechanism we subsequently address. While lower crustal relaxation can generate vertical deformation signals resembling those observed, prioritizing the fit to vertical deformation leads to a quite poor fit in the horizontal component (Supplementary Fig. 8).

We accounted for viscoelastic relaxation in both the lower crust and upper mantle using information about the subsurface layered structure (Fig. 1b) derived from available seismic velocity data46. The shear-wave velocity v s profile beneath the Arabian Plate shows anomalously high values (~4.6km/s) at 30–55km depths24, suggesting a mechanically strong layer. We assigned a high viscosity of 1×10²⁰ Pa·s to this layer, having tested alternative viscosities and examined scenarios where the Arabian and Anatolian Plates share identical layered structures (Supplementary Fig. 9). These tests confirmed higher rigidity of this layer compared to adjacent strata. Through grid search optimization (Supplementary Fig. 6), we determined viscosities for other layers using the PSGRN/PSCMP software (see “Methods”) (Fig. 1b), revealing that the Arabian Plate’s lower crust (η M=1×1020 Pa·s) is substantially more rigid than the Anatolian Plate’s (η M=2×1019 Pa·s, η K=6×1018 Pa·s). In the upper mantle, the steady-state viscosity η M is similar on both sides (η M=2×1018 Pa·s) except for the strong Arabian layer at 30–55km, but a lower transient viscosity (η K=2×1017 Pa·s) representing faster postseismic relaxation is found on the Anatolian side compared to the Arabian plate (η K=6×1017 Pa·s). This lateral heterogeneity in structure and viscosity correlates well with observed spatiotemporal asymmetry of EW displacements. The viscoelastic model predictions achieve an excellent fit to the EW and NS deformation fields, as well as their temporal evolution (Fig. 2 and Supplementary Fig. 10), confirming the validity of the viscoelastic model and suggesting a minimal contribution from afterslip (Supplementary Figs. 11 and 12).

The viscoelastic model, while providing a good fit to horizontal deformation, still exhibits significant residuals in the vertical component (Fig. 2g). This indicates that poroelastic rebound also plays an important role in the postseismic deformation. Based on the coseismic slip model (Fig. 3a), we isolated the contribution of poroelastic rebound by computing the difference in forward-modeled vertical surface deformation under undrained and drained Poisson’s ratios in the Green’s function relating slip to surface displacement (see "Methods"). Our analysis reveals that during the postseismic period, the highest correlation coefficient of 0.76 between model predictions and observations (Fig. 2h) is achieved when the difference in undrained and drained Poisson’s ratios reaches 0.1-0.13 (Fig. 3d). Although the vertical residuals from the poroelastic rebound model remain non-negligible (Supplementary Fig. 5c), the model successfully captures the first-order pattern of the observed deformation.

Absolute values of Poisson’s ratio under undrained and drained conditions are usually considered in poroelastic rebound modeling4,47,48. However, we argue that specifying absolute values is unnecessary as long as the chosen value falls within a reasonable range (e.g., 0.2–0.35). Poroelastic rebound-induced vertical deformation is more sensitive to the change in Poisson’s ratio between the undrained and drained states than to their absolute values (Fig. 3c). While different assumptions regarding the undrained Poisson’s ratio could theoretically influence the estimated coseismic slip, and consequently the poroelastic deformation, we find this effect to be negligible (Fig. 3b). Furthermore, our poroelastic rebound modeling indicates lateral heterogeneity across the EAF, with Poisson’s ratio changes of 0.13 on the Arabian side and 0.1 on the Anatolian side (Fig. 3d). This contrast results in slightly larger poroelastic vertical rebound on the Arabian plate. A Poisson’s ratio change Δν=ν u −ν d as large as 0.1–0.13 has rarely been reported in previous case studies, where typical variations are only about half that magnitude (~0.05)4,47. This suggests substantial coseismic and then postseismic change in crustal material properties in the region caused by these two large earthquakes, particularly within the stiffer Arabian plate.

Methods
To capture the postseismic deformation from the 2023 Kahramanmaraş earthquakes, we employed time-series InSAR and BOI (Burst Overlap Interferometry) techniques to derive surface displacements along both the satellite line-of-sight (LOS) and azimuth directions. The Sentinel-1 dataset used in this study includes three ascending and three descending tracks, each comprising 4–5 frames. We utilized both short-temporal-baseline interferograms and interferograms with longer temporal baselines (around 1 year) but short spatial baselines (e.g., less than 15m) (Supplementary Fig. 15). Pre-processing steps, including SAR image co-registration, interferogram generation, and phase unwrapping, were performed using GAMMA software.

To derive the LOS displacement from time-series InSAR data with improved signal-to-noise ratio, we first applied a multi-looking process with a 30×10 (range×azimuth) window to the interferograms. An adaptive filtering technique was then used, with a filter factor of 0.4 and a window size of 64. Phase unwrapping was subsequently performed using the minimum cost flow (MCF) method. After unwrapping, we employed the Small Baseline Subset (SBAS) technique to invert the time-series deformation for each satellite track. Due to differences in temporal sampling between tracks, the time reference of the time-series deformation for each track also varies. To unify the temporal reference across tracks and to account for possible period non-tectonics signals, we fitted each track’s time-series deformation using the empirical function d=a*log 1+t/T +b*sin 2πt/T +c (Supplementary Fig. 16). By removing the fitted sinusoidal component and the constant offset from the original displacement data, we isolate the postseismic displacement time series (blue markers), referenced to the moment immediately after the earthquake. Although early postseismic deformation tends to be large and difficult to fit accurately, our results show that the cumulative deformation from different tracks are generally consistent within overlapping regions and are also in good agreement with GNSS observations (Supplementary Fig. 4). This consistency supports the effectiveness of our processing strategy. To further derive east-west and vertical time-series deformation from ascending and descending LOS observations, we interpolated each track’s time-series data onto a unified temporal sampling grid (Supplementary Fig. 16). For each epoch, we then applied a simple least-squares method to estimate the east-west and vertical components by neglecting north-south deformation75. Previous studies have shown that while this approach yields accurate east-west deformation, the presence of true north-south deformation can introduce systematic errors into the derived vertical deformation76. To obtain more accurate vertical deformation fields, we corrected the vertical component using the modeled north-south deformation (Fig. 2d) from our viscoelastic forward model76.

The BOI processing is generally similar to the workflow of time-series InSAR. After generating Single Look Complex (SLC) images from Sentinel-1 data in the Burst Overlap (BO) region using GAMMA software, we performed differential interferometric processing. For both forward- and backward-looking interferograms in the BO region, the same topographic phase was applied for terrain correction. Prior to differencing the forward- and backward-looking interferograms to obtain the BOI interferogram, we applied adaptive filtering to enhance the signal-to-noise ratio of the resulting BOI interferogram. Following the generation of the BOI interferogram, we employed the SBAS method to invert the time-series deformation along the azimuth direction. To derive the north-south (NS) deformation in the BO region, we corrected for the east-west (EW) displacement (i.e., Fig. 2a) in the BOI azimuth displacement based on the formula of BOI=sin(α)×EW+cos(α)×NS with α the azimuth angle of SAR satellite.

The ascending and descending LOS and azimuth displacements derived from the above data processing contain both postseismic deformation and long-term secular deformation. To eliminate the influence of secular deformation on the analysis of the postseismic process, we obtained interseismic displacements by processing Sentinel-1 images spanning 2015 to 2023. The processing steps for the interseismic data were largely consistent with those used for the postseismic period. Notably, a Mw 6.8 earthquake occurred along the East Anatolian Fault (EAF) on January 24, 202077. To remove the coseismic deformation introduced by this event from the estimation of the interseismic deformation rate, we fitted the interseismic time series using a Heaviside step function (see Supplementary Fig. 17). In this analysis, we did not account for the postseismic deformation from this event. Given the moderate magnitude of the earthquake and the relatively long duration of our time series (approximately 8 years), the postseismic deformation contribution of this 2020 event to the estimated interseismic rate can be reasonably neglected. After obtaining the LOS and azimuth interseismic displacement, we subtracted the secular components from the postseismic time series. Based on these corrected postseismic LOS and azimuthal measurements, we then calculated the postseismic east-west, north-south, and vertical displacements.

We employed the PSGRN/PSCMP software78 to model postseismic viscoelastic relaxation. This software requires a coseismic slip model as the driving force for postseismic deformation and divides the subsurface structure into multiple layers with laterally homogeneous parameters for each layer. Consequently, conventional modeling using this software cannot account for the lateral heterogeneities in both subsurface layering and viscosity coefficients across our study region for this earthquake. Our strategy involves separately determining the viscosity coefficients for the Arabian Plate and Anatolian Plate on either side of the fault, while maintaining the software’s default homogeneous subsurface structure configuration. Through our experiments, we found that when determining the optimal viscosity coefficient for the Arabian Plate, we configured the software’s subsurface layering to match the Arabian Plate’s structure (Fig. 1b). This configuration yielded optimal fitting between the forward-modeled surface deformation field and observations in the Arabian Plate region (south of the EAF) (Supplementary Fig. 18e), although the fit was suboptimal in the Anatolian Plate region (north of the EAF). Similarly, when determining the optimal viscosity coefficient for the Anatolian Plate, we aligned the software’s subsurface layering with the Anatolian Plate’s structure (Fig. 1b), achieving optimal fitting between modeled and observed surface deformation in the Anatolian Plate region (Supplementary Fig. 18f), despite poorer fitting in the Arabian Plate region. Supplementary Fig. 18 demonstrates the comparison between observations and modeled values when optimal fitting is achieved for each plate, respectively. The results show distinct differences in fitting quality across the fault, while demonstrating relatively consistent fitting within the target plate.

To obtain a complete and optimal viscoelastic model-predicted deformation field, we spatially merged the deformation fields that achieved optimal fitting for the Arabian and Anatolian (Supplementary Fig. 18e, f). Specifically, we defined two normalized weight matrices, W₁ and W₂. For matrix W₁, the value is set to 1 at locations corresponding to the Arabian Plate (south of the EAF), decreases smoothly to 0 at locations north of the second seismic fault (i.e., Sürgü Fault), and transitions gradually between 1 and 0 in the intermediate region between the EAF and Sürgü Fault. The spatial characteristics of the W₁ matrix are illustrated in Supplementary Fig. 18b, and the W₂ matrix is simply defined as W₂=1−W₁ (Supplementary Fig. 18c). Let D₁ represent the optimal deformation field for the Arabian Plate region (Supplementary Fig. 18e) and D₂ the optimal deformation field for the Anatolian Plate (Supplementary Fig. 18f), the merged optimal deformation field D (Supplementary Fig. 18d) is then calculated through the weighted combination: D=W₁*D₁+W₂*D₂ (where * denotes element-wise multiplication). Although we simply set the weighting factors between the two major earthquake ruptures to vary from 0 to 1, resulting in relatively larger modeling residuals within the transition region compared with areas outside the two major fault ruptures (Supplementary Fig. 5), the final merged displacement field is still able to capture the first-order postseismic deformation signals (Fig. 2).

We employed a grid search approach to determine the optimal viscosity parameters for the lower crust and upper mantle beneath both the Arabian and Anatolian plates. The search ranges and optimal values for each parameter are presented in Supplementary Tables 1 and 2. During the parameter optimization for both plates, we conducted forward modeling calculations for all possible parameter combinations to ensure comprehensive coverage of the parameter space.

In this study, the use of laterally homogeneous viscoelastic models to approximate a heterogeneous lithospheric structure implicitly assumes that postseismic viscoelastic flow beneath one side of the fault has a negligible influence on surface deformation on the opposite side. While we could avoid this assumption by using 3D finite element models, such models are computationally expensive and not well-suited for extensive parameter-space exploration. However, to validate the reliability of this assumption, we performed a series of independent finite element simulations (using PyLith79 see Supplementary Discussion 1) and demonstrated that the cross-fault influence of postseismic viscoelastic flow is negligible.

We use the coseismic slip model (Fig. 3a) and elastic half-space dislocation theory to compute vertical displacements associated with coseismic fault slip under undrained and drained Poisson’s ratios (ν u and ν d). The difference between the two vertical displacement fields is taken as the forward model of poroelastic rebound. We fix ν u=0.3 and vary the Poisson’s ratio change Δν=ν u −ν d from 0.01 to 0.2 with an increment of 0.01. For each Δν, we calculate the corresponding surface deformation field and compare the predicted and observed vertical displacements to identify the Δν that best fits the observations.

In addition, lateral heterogeneity in the crust across the fault may lead to different poroelastic rebound responses on the two sides of the fault. Therefore, we adopt a strategy similar to that used in our viscoelastic modeling, estimating the optimal Poisson’s ratio change on each side separately. This strategy is justified because, according to our poroelastic modeling framework, the poroelastic rebound deformation at a surface point depends only on the slip model, the location of the surface point, and the Poisson’s ratio of the medium between the fault and that point. In other words, Poisson’s ratio on the opposite side of the fault has negligible effect on the rebound deformation on the side being evaluated. By independently comparing the predicted and observed rebound signals on both sides, we find that the optimal Poisson’s ratio change Δν is 0.13 for the Arabian plate and 0.1 for the Anatolian plate (Fig. 3d).

To quantitatively assess the sensitivity of Poisson’s ratio changes constrained by the current InSAR vertical deformation data, we perform a synthetic experiment. We get the vertical deformation signal (Supplementary Fig. 19a) of the poroelastic rebound by using the Poisson’s ratio change of 0.12, and add spatially correlated noise (Supplementary Fig. 19b), with a fractal dimension of 2.2 and an amplitude range of [−5cm, 5cm], to generate synthetic observations (Supplementary Fig. 19c). This type of noise is commonly used to simulate atmospheric artifacts in InSAR data80 and can reasonably reproduce the residual signals present in real observations (Supplementary Fig. 5c). The chosen amplitude range ensures that the noise level is comparable to the misfit between the real data and model predictions (Fig. 3d). Applying the same search procedure as above yields one estimated Poisson’s ratio change; repeating this process 10,000 times results in a variance of 0.015 for the estimated value (Supplementary Fig. 19d). We therefore conclude that the Poisson’s ratio change constrained by the current InSAR data is 0.121±0.015, corresponding to a relative uncertainty of 12%.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Review of “Interplate rheological contrast revealed by asymmetric deformation after the 2023 Kahramanmaraş earthquakes” by J Liu et al.\n\nThis is a nicely written paper. The science is well described and at the forefront of the field (particularly the interferometric observations), and results are quite interesting.\n\nThe only, pretty minor, concern I have is with the authors’ use of the verb “validates” in discussing their modeling strategy. They use a laterally homogeneous simulation tool to approximate the deformation resulting from laterally heterogeneous viscoelastic structures by, in essence, modeling the deformation on each side of the fault independently. The discussion in lines 336-357 accurately describes the modeling process, but I do not agree that this is a validation of the strategy (it is a validation that the viscoelastic structure is laterally heterogeneous).. In modeling each side with a laterally homogeneous model, you are neglecting the deformation that would result from the lateral coupling ratios between the two sides. Without doing the calculations, I agree with the authors implicit contention that those coupling ratios would be less important than the viscoelasticity structure at depth, but I feel that the only true validation of their modeling strategy would be to simulate postseismic deformation using a model that included lateral viscoelasticity contrasts, and then attempt to recover those input viscoelasticity values using their two-step estimation process. I’m not sure what I would suggest as the correct word here, perhaps just remove this final sentence (“This validates the feasibility of our strategy ...”) and let the description stand on its own. I would also recommend that you include a sentence in this paragraph (either following the first or second sentence) that states that the need to follow this strategy is due to the fact that the computational cost of fully 3D finite element models would be too prohibitive to fully explore the model space.", "reviewer_2": "The manuscript by Liu et al. presents an analysis of post-seismic deformation following the 2023 Kahramanmaraş earthquakes by combining InSAR observations with viscoelastic and poroelastic modeling. The authors propose that the observed asymmetric deformation is primarily due to rheological contrasts across the East Anatolian Fault, specifically suggesting that the Arabian plate is mechanically stiffer than the Anatolian plate.\n\nOverall, the manuscript is well-structured, clearly written, and represents a substantial effort. However, I believe that several major aspects require clarification, further discussion, and revision before the manuscript is suitable for publication.\n\nMajor Comments:\n1. Introduction – Missing Key References and Context.\nThe authors should expand the introduction to better contextualize their findings within the framework of previous work. The concept of asymmetric deformation and rheological contrast across the East Anatolian Fault been previously proposed in the context of escape tectonics, where the Arabian plate acts as a rigid indenter. This has been discussed in studies such as Şengör et al. (1985), Dewey et al. (1986), Faccenna et al. (2013), among others. These works should be cited and discussed to properly situate the current study within the existing literature.\n\n2. Viscoelastic Modeling – Use of Homogeneous Models for Heterogeneous Structures.\nA significant concern lies in the use of a modeling approach that assumes homogeneous viscoelastic media to represent an heterogeneous lithospheric structure. This simplification warrants further justification and discussion. Validation steps, such as synthetic or checkerboard tests, should be included to demonstrate the robustness of the results. Moreover, the use of weighting factors (e.g., w = 0 or 1) to represent lateral variations in rheology requires clarification. How is this binary weighting justified physically? Why does the rheology of one plate seemingly not influence the deformation pattern of the adjacent plate? These modeling assumptions should be more transparently discussed and quantitatively tested.\n\n3. Selection of Data and Model Sensitivity.\nThe decision to use only east-west (EW) displacements at distances greater than 50 km is questionable. This approach oversimplifies the inherently three-dimensional nature of the deformation signal and neglects near-field viscous flow, where the signal is strongest. As a result, the model lacks sensitivity to key features of both the upper and lower crust, as well as the upper mantle. Notably, the model assumes that the upper 20 km of the lithosphere behave as purely elastic, thereby excluding them from contributing to viscous flow. Below this depth and down to 55 km (or 42 km beneath the Anatolian plate) a pure (or nearly pure) Maxwell rheology is applied, with a relaxation time of approximately (!)100 years (or ~20 years for Anatolia). Therefore, the deformation at these depths has minor affect on the surface displacements over the relatively short, ~3-year observational period, presented in this manuscript. These simplifying assumptions limit the model’s ability to accurately constrain viscoelastic properties within the crust and upper mantle. A more comprehensive treatment that accounts for depth-dependent rheology and includes the full displacement field is warranted.\n\n4. Poroelastic Modeling – Lack of Methodological Details and Justification.\nThe manuscript does not sufficiently describe the poroelastic modeling approach in the Methods section. It is also unclear why poroelastic and elastic properties are assumed to be homogeneous. Given the heterogeneity of the region, this assumption may significantly bias the results. The authors should discuss how heterogeneity in poroelastic parameters (i.e. drained and undrained Poisson ratio) could influence the predicted deformation and whether any sensitivity tests were conducted.\n\nAdditional Comments:\n5. The BOI method has already been applied to the 2023 Kahramanmaraş earthquakes in Magen et al. (2024). This reference should be acknowledged in the introduction.\n6. Line 77 and the Supplementary Material, please re-check the expression for alpha. It appears that you may have used μ_k instead of μ_m, which would lead to α = 0 for a Maxwell body. This correction is essential, and should follow the convention used in Wang’s PSGRN/PSCMP code.\n7. Lines 109–110 use potentially confusing terminology. Rather than \"stronger\" and \"weaker\", it would be clearer to refer to \"higher” or “greater\" and \"lower\" deformation. For example, the sentence could be rephrased to clarify that higher deformation is observed in the mechanically weaker plate.\n8. In Line 112, the use of a logarithmic decay function to fit post-seismic deformation should be supported by appropriate references. Studies such as Freed et al. (2010) and many others have used similar approaches and should be cited accordingly.\n9. Figure 2. Please consider adding corrected EW and NS displacement panels to Figure 2 (as was done for panel g), which include the poroelastic response. This would allow for a more complete assessment of the model-data fit.\n10. Figure 2. The label “13.1” should be corrected to “-13.1” to accurately reflect the data."}, "author_response": "Reviewer #1\nThis is a nicely written paper. The science is well described and at the forefront of the field\n(particularly the interferometric observations), and results are quite interesting.\nThe only, pretty minor, concern I have is with the authors’ use of the verb “validates” in\ndiscussing their modeling strategy. They use a laterally homogeneous simulation tool to\napproximate the deformation resulting from laterally heterogeneous viscoelastic structures by,\nin essence, modeling the deformation on each side of the fault independently. The discussion\nin lines 336-357 accurately describes the modeling process, but I do not agree that this is a\nvalidation of the strategy (it is a validation that the viscoelastic structure is laterally\nheterogeneous).. In modeling each side with a laterally homogeneous model, you are\nneglecting the deformation that would result from the lateral coupling ratios between the two\nsides. Without doing the calculations, I agree with the authors implicit contention that those\ncoupling ratios would be less important than the viscoelasticity structure at depth, but I feel\nthat the only true validation of their modeling strategy would be to simulate\npostseismic deformation using a model that included lateral viscoelasticity contrasts, and then attempt to\nrecover those input viscoelasticity values using their two-step estimation process. I’m not sure\nwhat I would suggest as the correct word here, perhaps just remove this final sentence (“This\nvalidates the feasibility of our strategy ...”) and let the description stand on its own. I would\nalso recommend that you include a sentence in this paragraph (either following the first or\nsecond sentence) that states that the need to follow this strategy is due to the fact that the\ncomputational cost of fully 3D finite element models would be too prohibitive to fully explore\nthe model space.\n>> Thank you for your recommendation. Following your suggestion, we removed the sentence\ncontaining the verb “validates” from the revised manuscript. In addition, we incorporated a set\nof finite element method (FEM) validation experiments, which demonstrate that postseismic\nviscoelastic flow beneath one side of the fault exerts only a very weak influence on surface\ndeformation on the opposite side, thereby supporting the effectiveness of our modeling strategy.\nAccordingly, we added the following statement to the main text and included the corresponding\nFEM results in the Supplementary Information. Detailed modifications can be found in our\nresponse to the comment #2 of Reviewer 2.\n“In this study, the use of laterally homogeneous viscoelastic models to approximate a\nheterogeneous lithospheric structure implicitly assumes that postseismic viscoelastic flow\nbeneath one side of the fault has a negligible influence on surface deformation on the opposite\nside. While we could avoid this assumption by using 3D finite element models, such models\nare computationally expensive and not well suited for extensive parameter-space exploration.\nHowever, to validate the reliability of this assumption, we performed a series of independent\nfinite element simulations (using PyLith82, see Supplementary Discussion 1) and demonstrate\nthat the cross-fault influence of postseismic viscoelastic flow is negligible.”\n\nReviewer #2\nThe manuscript by Liu et al. presents an analysis of post-seismic deformation following the\nporoelastic modeling. The authors propose that the observed asymmetric deformation is\nprimarily due to rheological contrasts across the East Anatolian Fault, specifically suggesting\nthat the Arabian plate is mechanically stiffer than the Anatolian plate.\nOverall, the manuscript is well-structured, clearly written, and represents a substantial effort.\nHowever, I believe that several major aspects require clarification, further discussion, and\nrevision before the manuscript is suitable for publication.\nMajor Comments:\n1. Introduction – Missing Key References and Context.\nThe authors should expand the introduction to better contextualize their findings within the\nframework of previous work. The concept of asymmetric deformation and rheological contrast\nacross the East Anatolian Fault been previously proposed in the context of escape tectonics,\nwhere the Arabian plate acts as a rigid indenter. This has been discussed in studies such as\nŞengör et al. (1985), Dewey et al. (1986), Faccenna et al. (2013), among others. These works\nshould be cited and discussed to properly situate the current study within the existing literature.\n>>We have added the following sentences and references in the revised manuscript.\n“The long-term tectonic evolution of Anatolian Plate has been interpreted by the westward\ntectonic escape driven by northward movement of Arabian Plate15,16, with the EAF acting as a\nkey transform boundary accommodating this motion17. Later geodynamic models further\nsuggested that this escape is sustained by large-scale mantle flow coupling westward slab pull\nin the Aegean Sea18-20.”\n15. Tapponnier P, Molnar P. Slip-Line Field-Theory and Large-Scale Continental Tectonics. Nature 264,\n319-324 (1976).\n16. Şengör A, Görür N, Şaroğlu F. Strike-slip faulting and related basin formation in zones of tectonic escape:\nTurkey as a case study. Strike-Slip Faulting and Basin Formation, edited by K. T. Biddle and N.\nChristie-Blick, Spec. Publ. SEPM Soc. Sediment. Geol. 37, 227-264 (1985).\n17. Dewey JF, Hempton MR, Kidd WSF, Saroglu F, Sengör AMC. Shortening of continental lithosphere:\nthe neotectonics of Eastern Anatolia - a young collision zone. Geol Soc Spec Publ 19, 3-36 (1986).\n18. Faccenna C, Becker TW, Jolivet L, Keskin M. Mantle convection in the Middle East: Reconciling Afar\nupwelling, Arabia indentation and Aegean trench rollback. Earth Planet. Sci. Lett. 375, 254-269\n(2013).\n19. Faccenna C, et al. Mantle dynamics in the Mediterranean. Rev. Geophys. 52, 283-332 (2014).\n20. Roche V, Jolivet L, Papanikolaou D, Bozkurt E, Menant A, Rimmelé G. Slab fragmentation beneath the\nAegean/Anatolia transition zone: Insights from the tectonic and metamorphic evolution of the\nEastern Aegean region. Tectonophysics 754, 101-129 (2019).\n2. Viscoelastic Modeling – Use of Homogeneous Models for Heterogeneous Structures.\nA significant concern lies in the use of a modeling approach that assumes homogeneous\nviscoelastic media to represent an heterogeneous lithospheric structure. This simplification\nwarrants further justification and discussion. Validation steps, such as synthetic or\ncheckerboard tests, should be included to demonstrate the robustness of the results. Moreover,\nthe use of weighting factors (e.g., w = 0 or 1) to represent lateral variations in rheology requires\nclarification. How is this binary weighting justified physically? Why does the rheology of one\nplate seemingly not influence the deformation pattern of the adjacent plate? These modeling\nassumptions should be more transparently discussed and quantitatively tested.\n>> We acknowledge the reviewer’s concern regarding the use of homogeneous viscoelastic\nmodels to approximate an intrinsically heterogeneous lithospheric structure. In this study, we\nadopt a modeling strategy in which the two sides of the fault are modeled separately using\nlaterally homogeneous viscoelastic structures. The implicit assumption underlying this\napproach is that postseismic viscoelastic flow occurring beneath one side of the fault exerts\nonly a minor influence on the surface deformation observed on the opposite side, and that such\ninfluence can be neglected to first order. To clarify this point, we added the following\nSupplementary Discussion:\n“Supplementary Discussion 1. Validation of the Homogeneous Modeling Strategy Using\nFinite Element Simulations\nTo evaluate the validity of approximating a laterally heterogeneous lithosphere using laterally\nhomogeneous viscoelastic models on each side of the fault, we performed a series of three-\ndimensional finite element forward simulations using PyLith. The finite element model (FEM)\nadopts the same layered lithospheric structure and rheological parameters as the preferred\nmodel shown in Fig. 1b. The computational domain extends 1800 km × 1800 km laterally and\n600 km in depth to minimize boundary effects on surface deformation in the study area. The\nfault geometry and coseismic slip distribution are identical to those used in the\nPSGRN/PSCMP-based modeling. A schematic illustration of the FEM mesh is shown in\nSupplementary Fig. 20.\nWe designed four sets of sensitivity experiments to quantify the influence of viscoelastic\nrheology on one side of the fault on the surface deformation observed on the opposite side.\nSpecifically, we varied the Kelvin viscosity 𝜂 of the lower crust and upper mantle beneath\n𝐾\neither the Anatolian or the Arabian plate, while keeping all other parameters fixed. The tested\n𝑜𝑝𝑡\n𝜂 values span one order of magnitude centered on the optimal value 𝜂 , with logarithmic\n𝐾 𝐾\nincrements of 0.25, i.e.,\n𝑜𝑝𝑡 𝑜𝑝𝑡\n𝜂 = 10\n[𝑙𝑜𝑔10(𝜂\n𝐾\n)−0.5:0.25:𝑙𝑜𝑔10(𝜂\n𝐾\n)+0.5]\n𝐾\nThe FEM-predicted surface deformation fields reproduce the main spatial characteristics of the\nobserved InSAR deformation (Supplementary Fig. 21). Sensitivity analyses at representative\npoints on both sides of the fault show that variations in 𝜂 beneath one plate result in only\n𝐾\nminor deformation changes on the opposite plate (Supplementary Fig. 22). Even in the most\nsensitive case, varying the lower-crustal 𝜂 beneath the Anatolian plate (Supplementary Fig.\n𝐾\n22a), the induced deformation change on the Arabian side is within about 2 cm. This\ndeformation amplitude is comparable to the uncertainty level of the InSAR time-series\nmeasurements used in this study. As the discussion of the grid search sensitivity in\nSupplementary Fig. 6, this sensitivity limitation does not affect the general conclusions of the\nstudy. Therefore, within the resolution of the available observations, the assumption that\npostseismic viscoelastic flow beneath one side of the fault has a negligible influence on surface\ndeformation on the opposite side is justified.”\nSupplementary Fig. 20. A schematic illustration of the finite element model (FEM) mesh used\nto validate our PSGRN/PSCMP modeling strategy.\nSupplementary Fig. 21. Comparison between the observed cumulative postseismic horizontal\ndisplacements and predicted displacements by the three-dimensional finite element model. The\nmesh of the 3D finite element model is shown in Supplementary Fig. 20, and the layered\nstructure and viscosity values are in Fig. 1b of the main text. Overall, the model predictions\nshow good agreement with the observations. Points B and B’ are two representative locations\nsituated on the Arabian and Anatolian plates, respectively. The time series at these two points\nunder different viscosity scenarios are shown in Supplementary Fig. 22. Selecting\ncharacteristic points on opposite sides of the fault allows us to investigate how variations in\nviscosity on one side of the fault influence surface deformation on the other side.\nSupplementary Fig. 22. Modeled post-seismic transient displacements at two representative\npoints, B and B’, located on the Arabian and Anatolian plates, respectively, for different\nsubsurface rheological parameters on one side of the fault. a Time series at points B’ and B for\ndifferent viscosity values assigned to the Kelvin body in the lower crust beneath the Anatolian\nplate. The results show that modifying subsurface parameters on the Anatolian side\nsignificantly affects surface deformation at the Anatolian-side point B’, while the impact on\nthe Arabian-side point B is negligible. b-d Time series at the two locations B and B’ resulting\nfrom variations in viscosity within the Anatolian upper mantle, the Arabian lower crust, and\nthe Arabian upper mantle, respectively. Overall, these results indicate that variations in\nsubsurface viscosity on one side of the fault—within one order of magnitude—have only a\nminor influence on the surface deformation on the opposite side of the fault. Different colored\ncurves represent different viscosity values: the yellow curve denotes the optimal viscosity, the\nblue curve indicates a viscosity 0.5 orders of magnitude lower than the optimal value, and the\ngreen curve represents a viscosity 0.5 orders of magnitude higher.\nWe also added the following sentences to the main text of the revised manuscript\n“In this study, the use of laterally homogeneous viscoelastic models to approximate a\nheterogeneous lithospheric structure implicitly assumes that postseismic viscoelastic flow\nbeneath one side of the fault has a negligible influence on surface deformation on the opposite\nside. While we could avoid this assumption by using 3D finite element models, such models\nare computationally expensive and not well suited for extensive parameter-space exploration.\nHowever, to validate the reliability of this assumption, we performed a series of independent\nfinite element simulations (using PyLith82, see Supplementary Discussion 1) and demonstrate\nthat the cross-fault influence of postseismic viscoelastic flow is negligible.”\nRegarding the simplified weighting scheme used to merge the independently and\nhomogeneously modeled results, we assign binary weights (0 or 1) to regions outside the two\nmajor fault ruptures. Within the region between the two major ruptures, the weighting values\nare prescribed to vary smoothly from 0 to 1. We acknowledge that this simplification\nintroduces modeling residuals; however, based on both above finite-element analyses and the\npresent modeling results (Fig. 2), the adopted strategy is capable of capturing the first-order\ncharacteristics of the postseismic observations. Although relatively larger residuals are\nobserved within the transition region, they do not affect our overall modeling conclusions. We\nhave added the following sentences to clarify this point.\n“Although we simply set the weighting factors between the two major earthquake ruptures to\nvary from 0 to 1, resulting in relatively larger modeling residuals within the transition region\ncompared with areas outside the two major fault ruptures (Supplementary Fig. 5), the final\nmerged displacement field is still able to capture the first-order postseismic deformation signals\n(Fig. 2).”\n3. Selection of Data and Model Sensitivity.\nThe decision to use only east-west (EW) displacements at distances greater than 50 km is\nquestionable. This approach oversimplifies the inherently three-dimensional nature of the\ndeformation signal and neglects near-field viscous flow, where the signal is strongest. As a\nresult, the model lacks sensitivity to key features of both the upper and lower crust, as well as\nthe upper mantle. Notably, the model assumes that the upper 20 km of the lithosphere behave\nas purely elastic, thereby excluding them from contributing to viscous flow. Below this depth\nand down to 55 km (or 42 km beneath the Anatolian plate) a pure (or nearly pure) Maxwell\nrheology is applied, with a relaxation time of approximately (!)100 years (or ~20 years for\nAnatolia). Therefore, the deformation at these depths has minor affect on the surface\ndisplacements over the relatively short, ~3-year observational period, presented in this\nmanuscript. These simplifying assumptions limit the model’s ability to accurately constrain\nviscoelastic properties within the crust and upper mantle. A more comprehensive treatment that\naccounts for depth-dependent rheology and includes the full displacement field is warranted.\n>>Thank you for your suggestion. We re-searched the optimal model parameters by comparing\nthe modeled and observed deformation fields in both the east–west and north–south\ncomponents across the entire study area. The comparison shows that the parameters derived\nfrom the full displacement field are consistent with those obtained using only the east–west\ncomponent beyond 50 km, as originally adopted in the manuscript (see the following figures).\nTherefore, this modification does not affect the overall conclusions of the study. Accordingly,\nwe have made the following revisions in the manuscript:\n“For the viscoelastic modeling, both EW and NS displacements are used. Given that near-field\ndeformation could be influenced by shallow fault slip, we conducted an independent test by\nsearching for the optimal model parameters using only data beyond 50 km (approximately\ntwice the seismogenic depth44-46). The results indicate that the inclusion of near-field\nobservations has an insignificant impact on the optimal model solution (Supplementary Figs.\n6-7).”\nRegarding the viscoelastic properties of different layers, our results show that the viscosity of\nthe upper mantle can be reasonably constrained with the available data. In contrast, the model\nis less sensitive to the viscoelastic parameters of the lower crust. This is primarily because the\nviscosity of the lower crust is higher than that of the upper mantle, resulting in a smaller\ncontribution to the observed surface deformation. Consequently, given the same signal-to-\nnoise ratio of observations, the sensitivity to the lower-crustal viscosity is relatively low.\nNevertheless, the grid-search results (the following Supplementary Fig. 6) reveal that, for the\nArabian Plate, the model–observation misfit decreases progressively as the lower crust\napproaches a Maxwell body behavior (η > 1×10²⁰ Pa·s). For the lower crust of Anatolian Plate,\na local minimum in misfit is found at η = 2×10¹⁹ Pa·s, although the misfit increases only slightly\nwhen η exceeds 2×10¹⁹ Pa·s. These results align well with the mechanically stiffer property of\nArabian Plate compared with the Anatolian Plate.\nFor the upper crust (0–20 km), the contribution from viscoelastic relaxation is expected to be\nminimal; therefore, viscous flow in this layer is not considered under the current observational\nconstraints.\nBased on these findings, we further refined the RMSE figure of the parameter search and\nperformed an additional sensitivity analysis in the caption of Supplementary Fig. 6 to illustrate\nhow the current observations constrain the model parameters.\nSupplementary Fig. 6. Root mean square (RMS) values for different parameters in the\nviscoelastic model. We determine the optimal viscoelastic parameters (Fig. 1b) separately for\na, b the Arabian plate and c, d the Anatolian plate. To better illustrate the sensitivity of the\nRMS values to different viscoelastic parameters, we show the RMS variations with respect to\nthe parameters of a, c the lower crust (LC) and b, d the upper mantle (UM), respectively. For\ninstance, in a, the RMS values for different lower crustal parameters beneath Arabia are shown\nwhile keeping the upper mantle parameters fixed at its optimal value. It should be noted that\nduring the grid search process, we explored a wide range of parameter combinations\n(Supplementary Tables 1-2) for both the lower crust and upper mantle, and performed forward\nmodeling for all possible combinations to compare with the observations and identify the\noptimal parameters (see Materials and Methods for details). Limited combinations are shown\nsolely for illustrative purposes. 𝜂 and 𝜂 are the steady-state viscosity (Maxwell body) and\n𝑀 𝐾\ntransient viscosity (Kelvin body) in a Burgers model, respectively, 𝜇 and 𝜇 are the\n𝑀 𝐾\ncorresponding shear modulus, and 𝛼 is the ratio between the effective and the unrelaxed shear\nmodulus (i.e., 𝛼 = 𝜇 /(𝜇 +𝜇 ) ). From the figure, it is evident that the available\n𝐾 𝐾 𝑀\nobservations are more sensitive to the viscoelastic parameters of the upper mantle, as indicated\nby the clear local RMS minima in panels b and d. In contrast, the model shows lower sensitivity\nto the viscoelastic properties of the lower crust. As shown in panel c, when the viscosity\nexceeds the optimal value, the RMS increases only slightly; similarly, in panel a, when the\nlower crust approaches a Maxwell body or even an elastic regime, the RMS value\ncorresponding to the optimal parameters does not show a pronounced advantage. This reduced\nsensitivity primarily arises because the viscosity of the lower crust is higher than that of the\nupper mantle, resulting in a weaker contribution of lower-crustal viscoelastic flow to the\nsurface deformation signal. Consequently, under the same signal-to-noise ratio conditions, the\ncurrent observations are less sensitive to the lower-crustal viscosity. Nevertheless, our data and\nmodeling results consistently indicate that the layers of the Arabian Plate are somewhat\nstiffer—i.e., characterized by a higher viscosity—than those of the Anatolian Plate.\nSupplementary Fig. 7. Same as Supplementary Fig. 6, except showing the misfit between the\nmodel and the displacement at distances greater than 50 km from the fault. The optimal\nviscosities obtained in this analysis are very close to those in Supplementary Fig. 6, indicating\nthat including near-fault data in the determination of the optimal viscosities does not\nsignificantly bias the results.\n4. Poroelastic Modeling – Lack of Methodological Details and Justification.\nThe manuscript does not sufficiently describe the poroelastic modeling approach in the\nMethods section. It is also unclear why poroelastic and elastic properties are assumed to be\nhomogeneous. Given the heterogeneity of the region, this assumption may significantly bias\nthe results. The authors should discuss how heterogeneity in poroelastic parameters (i.e.\ndrained and undrained Poisson ratio) could influence the predicted deformation and whether\nany sensitivity tests were conducted.\nThank you for your suggestion. In the revised manuscript, we have added detailed descriptions\ncovering (1) the poroelastic rebound modeling procedure, (2) the way we account for the\nheterogeneity of poroelastic rebound between the two sides of the fault, and (3) the sensitivity\nanalysis of poroelastic rebound parameters based on numerical experiments.\nIn the Results section:\n“Furthermore, our poroelastic rebound modeling indicates lateral heterogeneity across the EAF,\nwith Poisson’s ratio changes of 0.13 on the Arabian side and 0.1 on the Anatolian side. This\ncontrast results in slightly larger poroelastic vertical rebound on the Arabian plate. A Poisson’s\nratio change ∆ν=ν -ν as large as 0.1-0.13 has rarely been reported in previous case studies,\nu d\nwhere typical variations are only about half that magnitude (~0.05)4,41. This suggests\nsubstantial coseismic and then postseismic change in crustal material properties in the region\ncaused by these two large earthquakes, particularly within the stiffer Arabian plate.”\nIn the Methods section:\n“We use the coseismic slip model (Fig. 3a) and elastic half-space dislocation theory to compute\nvertical displacements associated with coseismic fault slip under undrained and drained\nPoisson’s ratios (𝑣 and 𝑣 ). The difference between the two vertical displacement fields is\nu d\ntaken as the forward model of poroelastic rebound. We fix 𝑣 =0.3 and vary the Poisson’s ratio\nu\nchange ∆𝜈 = 𝑣 −𝑣 from 0.01 to 0.2 with an increment of 0.01. For each ∆𝜈, we calculate\nu d\nthe corresponding surface deformation field and compare the predicted and observed vertical\ndisplacements to identify the ∆𝜈 that best fits the observations.\nIn addition, lateral heterogeneity in the crust across the fault may lead to different poroelastic\nrebound responses on the two sides of the fault. Therefore, we adopt a strategy similar to that\nused in our viscoelastic modeling, estimating the optimal Poisson’s ratio change on each side\nseparately. This strategy is justified because, according to our poroelastic modeling framework,\nthe poroelastic rebound deformation at a surface point depends only on the slip model, the\nlocation of the surface point, and the Poisson’s ratio of the medium between the fault and that\npoint. In other words, the Poisson’s ratio on the opposite side of the fault has negligible effect\non the rebound deformation on the side being evaluated. By independently comparing the\npredicted and observed rebound signals on both sides, we find that the optimal Poisson’s ratio\nchange ∆𝜈 is 0.13 for the Arabian plate and 0.1 for the Anatolian plate (Fig. 3d).\nTo quantitatively assess the sensitivity of Poisson’s ratio changes constrained by the current\nInSAR vertical deformation data, we perform a synthetic experiment. We get the vertical\ndeformation signal (Supplementary Fig. 19a) of the poroelastic rebound by using the Poisson’s\nratio change of 0.12, and add spatially correlated noise (Supplementary Fig. 19b), with a fractal\ndimension of 2.2 and an amplitude range of [-5 cm, 5 cm], to generate synthetic observations\n(Supplementary Fig. 19c). This type of noise is commonly used to simulate atmospheric\nartifacts in InSAR data83 and can reasonably reproduce the residual signals present in real\nobservations (Supplementary Fig. 5c). The chosen amplitude range ensures that the noise level\nis comparable to the misfit between the real data and model predictions (Fig. 3d). Applying the\nsame search procedure as above yields one estimated Poisson’s ratio change; repeating this\nprocess 10,000 times results in a variance of 0.015 for the estimated value (Supplementary Fig.\n19d). We therefore conclude that the Poisson’s ratio change constrained by the current InSAR\ndata is 0.121±0.015, corresponding to a relative uncertainty of 12%.”\nSupplementary Fig. 19. An example of the simulated vertical observations used to test the\nsensitivity of Poisson’s ratio changes constrained by the current InSAR vertical deformation\ndata. a Simulated vertical displacement due to poroelastic rebound with a Poisson’s ratio\nchange of 0.12. b Simulated noise designed to approximate the vertical residual signals\n(Supplementary Fig. 5c) in the real observations. c Simulated observations used to determine\nthe optimal value of the Poisson’s ratio change in the numerical experiments. a-c illustrate one\nexample of how simulated observations are generated for a single estimation of the Poisson’s\nratio change. We repeated the simulation 10,000 times, yielding 10,000 optimal estimates of\nthe Poisson’s ratio change, as shown in d. In these simulations, the poroelastic vertical\ndisplacement shown in a is kept identical, whereas the noise realization in b and the resulting\nsimulated observations in c vary for each run.\nAdditional Comments:\n5. The BOI method has already been applied to the 2023 Kahramanmaraş earthquakes in\nMagen et al. (2024). This reference should be acknowledged in the introduction.\nWe have updated the BOI references in the revised manuscript.\n36. Magen Y, et al. Fault Coalescence, Slip Distribution, and Stress Drop of the February 2023 Southeast\nTurkiye Earthquakes from Joint Inversion of SAR, GNSS, and Burst Overlap Interferometry.\nSeismol. Res. Lett. 95, 680-696 (2024).\n37. Nergizci M, et al. Refining 3D Displacement Fields and Coseismic Slip Models of the 2023\nKahramanmaraş Earthquakes Using Subswath and Burst Overlap Interferometry (SBOI). Journal of\nGeophysical Research-Solid Earth 130, e2025JB031421 (2025).\n6. Line 77 and the Supplementary Material, please re-check the expression for alpha. It appears\nthat you may have used μ_k instead of μ_m, which would lead to α = 0 for a Maxwell body. This correction is essential, and should follow the convention used in Wang’s PSGRN/PSCMP code.\nWe have double-checked Wang’s PSGRN/PSCMP code and its documentation. We confirm\nthat our mathematical definition of 𝛼 = 𝜇 /(𝜇 +𝜇 ) and the statement that “α = 1.0 means\n𝐾 𝐾 𝑀\na pure Maxwell body” are correct. Here is the screenshot of the software documentation.\nα = 1 means that the shear modulus μ of the Kelvin body approaches infinity. Although such\nK\na value is not physically realistic, in the PSGRN/PSCMP implementation, μ can be taken as\nK\nan effective shear modulus used to regulate the strength of the short-term Kelvin relaxation\nrather than the true long-term shear modulus. In this case, because the spring and viscous\nelements in the Kelvin body are arranged in parallel, an infinite spring stiffness causes the\nKelvin body to contribute no postseismic deformation. Then, the whole Burgers rheology is\nequivalent to a Maxwell body.\n7. Lines 109–110 use potentially confusing terminology. Rather than \"stronger\" and \"weaker\",\nit would be clearer to refer to \"higher” or “greater"}, {"round": 2, "reviewer_comments": {"reviewer_2": "I have no further comments on the revised version of the manuscript. Thank you.\nYariv Hamiel"}, "author_response": ""}]

Decision letter