Detectability of low-viscosity zone along lithosphere–asthenosphere boundary beneath the Nankai Trough, Japan, based on high-fidelity viscoelastic simulation

After the 2011 Tohoku-oki earthquake, a seafloor observation system observed rapid landward deformation. These seafloor observation data are potentially explicable via modeling of a low-viscosity zone under the subducting plate, called the lithosphere–asthenosphere boundary (LAB). However, the effect of a low-viscosity LAB has not been empirically examined in past earthquakes because of the absence of seafloor observation data, which are expected to have high sensitivity in the detection of a low-viscosity LAB. The Nankai Trough, southwest Japan, is one location where seafloor geodetic stations are in operation that could detect a low-viscosity LAB after a future earthquake. Therefore, we quantitatively model how the low-viscosity of the LAB is expected to affect postseismic crustal deformation following an anticipated future large earthquake under the Nankai Trough. Calculating viscoelastic deformation using a model that takes into account the realistic crustal structure in the southwestern Japanese Archipelago, we examine the effects of LAB viscosity on surface deformation patterns. The results show for very low LAB vis-cosities ( 2.5 × 10 17 Pa s), rapid landward displacement as large as 38 cm/yr occurs in the offshore area, with uplift and subsidence patterns highly dependent on the viscosity structure. Especially in the first year after the earth-quake, the difference in deformation between the cases with and without a low-viscosity LAB was much larger than the observational error level of the current seafloor observation system. For such low LAB viscosities, the probability of detection is therefore high. On the other hand, for moderately low LAB viscosities ( 2.5 × 10 18 Pa s), the deformation patterns in the cases with and without a low-viscosity LAB are similar. Indeed, the difference in displacement lies within the observational error level of the current observational network, so detection is expected to be difficult for a LAB of only moderately low-viscosity. However, such a LAB may be detected by improving observational accuracy in the Global Navigation Satellite System–acoustic observation technique by performing more frequent measurements. Additional detection potential lies in expected future accuracy improvements in vertical displacement estimates at the DONET stations.


Introduction
After a large earthquake occurs in a plate subduction zone, additional crustal deformation follows in several to tens of years.This postseismic crustal deformation is typically attributed to viscoelastic stress relaxation in the asthenosphere (Wang et al. 2012) and afterslip (Hsu et al. 2006) after the mainshock.For example, the 2011 M9 Tohoku-oki earthquake in the Japan Trench caused continuous crustal deformation for more than ten years afterward.Rapid landward deformation was observed immediately after the earthquake based on seafloor crustal deformation data obtained by the Global Navigation Satellite System-acoustic combination technique (GNSS-A) (Sun et al. 2014;Watanabe et al. 2014).Relocking of the Tohoku-oki source region cannot explain the data, because the magnitude of this crustal deformation far exceeds the long-term rate of plate convergence, even if the plate boundary relocked immediately after the earthquake.Nor can the afterslip effect explain the landward motion, since afterslip causes deformation in the opposite direction (Watanabe et al. 2014).
To fully explain the seafloor observation data, relaxation phenomena within the low-viscosity region under the subducting plate, called the lithosphere-asthenosphere boundary (LAB), seems to be required.Therefore, it has been proposed that the LAB be modeled, along with the viscosity structure of the mantle wedge and oceanic mantle (e.g., Sun et al. 2014;Freed et al. 2017).
The low-viscosity region has been ascribed to either concentrated melt zones (Kawakatsu et al. 2009;Rychert and Shearer 2009;Schmerr 2012) or combination of lower temperature plastic flow and dislocation creep in the lower part of a bending slab (Buffett and Becker 2012).Both hypotheses suggest that a low-viscosity LAB may be universally present at subduction zones.To prove the existence of a low-viscosity LAB at each subduction zone, it requires seafloor geodetic observation data as obtained along the Japan Trench.
The Nankai Trough is located off southwest Japan, where the Philippine Sea Plate is subducting under the southwest Japan Arc.It is a major plate subduction zone where large earthquakes are expected to occur and have occurred.The most recent ruptures in the region were the 1944 Tonankai earthquake (M7.9) and the 1946 Nankai earthquake (M8.0).Historical records reveal that megathrust earthquakes have occurred repeatedly along the Nankai Trough with recurrence intervals of 100-200 yr (Ando 1975).Therefore, a major earthquake could occur within the next few decades.Around the Nankai Trough, an advanced seafloor observation network has been developed, and interseismic crustal deformation data are available from the GNSS-A (Yokota et al. 2016) and ocean bottom pressure (OBP) sensors of the Dense Oceanfloor Network system for Earthquake and Tsunamis (DONET) (Kaneda et al. 2015).Thus, crustal deformation patterns associated with the low viscosity of the LAB may be observed in connection with a future ∼ M8 event in the Nankai Trough, as were observed in the 2011 Tohoku-oki case in the Japan Trench.Notably, estimations of afterslip and viscosity have been carried out for the 1944 Tonankai and 1946 Nankai earthquakes in order to provide a rational explanation for postseismic deformation data (Sherrill and Johnson 2021).However, land deformation data from triangulation and leveling surveys are the only data available for these earthquakes, and potential viscosity heterogeneities, including at the LAB, have not been examined.A study of the viscoelastic structure of this region, including the influence of the LAB, has been made by Suito (2017) for the 2004 M7.1 and M7.4 earthquakes off the Kii Peninsula.However, the misfit scores in his analysis do not resolve the viscosity structure well probably because he used only 10 stations in the vicinity of the epicenter.
Therefore, we quantitatively evaluate how the low viscosity of LAB affects postseismic crustal deformation in association with a hypothetical future large earthquake in the Nankai Trough.We implement the crustal deformation analysis using a highly detailed model that realistically reflects the crustal structure around the Nankai Trough.The computational cost for analysis using this highly detailed crustal structure model would be too high to conduct with standard analysis methods.However, such analysis has been made possible by the development of a crustal deformation analysis method that reduces computational cost (Fujita et al. 2017) and improvement in the computing capabilities of high-performance computing systems.Recently, several studies have performed viscoelastic analyses using highly detailed crustal structure models (e.g., Yamaguchi et al. 2018).In the present study, we also perform a highly detailed crustal deformation analysis using GPUs with newly improved computational performance to quantitatively show the dependence of amount of deformation on the viscosity of the LAB.
In the following sections, we first describe the crustal deformation analysis methodology and problem setting.We also explain the results of the analysis.Then we discuss the detectability of the LAB and the validity of the analysis.We summarize the paper's conclusions.

Problem setting
In this study, we evaluated the viscoelastic response of a M9-class earthquake scenario with a large slip along a shallow plate boundary (Hyodo et al. 2013) (Fig. 1).The direction of slip was set to N 125 • E (Miyazaki and Heki 2001).We computed the 4-year viscoelastic response which is sufficiently long to observe the displacements associated with rapid viscous relaxation within the low-viscosity LAB as observed in the 2011 Tohoku-oki case (Sun et al. 2014;Freed et al. 2017).We targeted a complex three-dimensional (3D) crustal structure model based on the Japan integrated velocity structure model (JIVSM) version 1 (Koketsu et al. 2009, 2012).In JIVSM, the plate boundary geometry and the elastic structure of the crust and mantle are considered, but the viscosity structure is not explicitly considered.Therefore, we created new layers for the upper plate, the subducting plate, the viscoelastic mantle and LAB according to the distribution of viscosity (Fig. 2 and Table 1).The upper plate and the subducting plate were assumed to be elastic based on Sherrill and Johnson (2021).The thickness of the upper plate was assumed to be 45 km based on the thickness of the cold wedge following Sherrill and Johnson (2021).Other studies have estimated the thickness of the lithosphere to be in the range of 30−60 km, which would explain the postseismic deformation data (e.g., Fukahata et al. 2004;Ito and Hashimoto 2004;Noda et al. 2018).The thickness of the subducting plate was assumed to be uniformly 30 km based on Sherrill and Johnson (2021).Although the low-viscosity LAB is supported by several studies, the possible thickness of low-viscosity LAB has not been discussed in detail.
The thickness of the LAB was set to 10 km based on Sun et al. (2014) and Suito (2017).The low-viscosity LAB was modeled up to the end of the subducting slab based on Suito (2017).Though past studies have applied several types of viscoelasticity to explain transient postseismic deformation (e.g., Sun et al. 2014), Maxwell (steady state) viscosity is the most essential in subduction dynamics.We thus focus on the Maxwellian viscosity based on Sherrill and Johnson (2021).
To evaluate the influence of (1) the presence of a low-viscosity LAB and (2) differences in viscosities of oceanic mantle and mantle wedge, we considered five cases, each with a different Maxwell viscosity (Table 2).Based on Suito (2017) and Sun et al. (2014), we considered the case of lower viscosity of the mantle wedge compared to that of the ocean mantle (e.g., case 2; Hirth and Kohlstedt 2003), the case of uniform mantle viscosity (case 1, 5), and the case of a low-viscosity LAB beneath the subducting plate (case 3, 4).In particular, the viscosity of the LAB was set to 2.5 × 10 17 Pa s in case 4 following Sun et al. (2014).
We then solved the system of equations in which the crust and mantle were each regarded as a linear viscoelastic body based on the Maxwell model: (1) (3) Here, σ and f represent the stress tensor and external force, respectively.(˙) , η , ǫ , and u are the first derivative in time, viscosity coefficient, strain, and displacement, respectively., µ , and δ are the first Lamé parameter, shear modulus, and Kronecker delta.The Einstein summation convention was used.Since the complex 3D geometry and heterogeneous material properties of crustal structures can have significant effects on crustal deformation, the analysis was performed by finite-element analysis based on unstructured elements, which can accurately consider stress-free boundary conditions at the ground surface as well as complex geometry and material heterogeneity.To evaluate stresses and strains with high accuracy, we used tetrahedral quadratic elements.Fault slip was evaluated using the split-node technique (Melosh and Raefsky 1981).The effect of gravity was evaluated by adding external forces corresponding to the ground surface displacement (Barbot and Fialko 2010).Dirichlet boundary conditions were imposed on the bottom and sides of the finite-element model.As in Ichimura et al. (2016), the evolution of the viscoelastic analysis was computed by using Algorithm 1.To perform crustal deformation analysis for 4 yr from the time of the earthquake, we set the time-step increment to dt = 648, 000 s (i.e., 7.5 days) and the number of time steps to N t = 200.and Table 3 (points 1-8) Algorithm 1 Algorithm for solving linear viscoelastic response of crust following Ichimura et al. (2016).Here, superscript () i is the variable in the i-th time step.dt is the time increment.B is the displacementstrain transformation matrix.A and D are matrix indicating matrix property.
,where β′ is the Jacobian matrix of β and α is the controlling parameter.Ω is the viscoelastic body.Detail of the algorithm is shown in Ichimura et al. (2016) To accurately evaluate the difference in crustal response due to the presence of an LAB with low viscosity and due to its particular viscosity, a large-scale 3D crustal deformation analysis was required considering a highly detailed, large, 3D rheological structure.Setting the origin at the zero-meter height of the Earth ellipsoid model GRS80 at ( 33.5 • N, 135 • E), we generated the finite-ele- ment model based on this crustal structure model for the region of −1248 km ≤ x ≤ 1248 km in the east-west direction, −1248 km ≤ y ≤ 1248 km in the north-south direction, and −1100 km ≤ z ≤ 1 km in the vertical direc- tion in a Cartesian coordinate system of mutually perpendicular axes.To ensure computational convergence around the fault plane, an initial mesh with 500 m-sized elements was used in the vicinity of the fault, while up to 2000 m-sized elements were used for domains far from the fault with elastic material properties.Figure 3 shows the generated finite-element model, which consists of 1.0 × 10 9 tetrahedral elements and has 4.2 × 10 9 degrees of freedom.
The computational cost of solving systems of linear equations K v δu i = f i in line 8 of Algorithm 1 at each time-step would be huge using a highly detailed largedomain model with 4.2 × 10 9 degrees of freedom; such models are typically considered difficult to analyze.In this study, we used the multi-grid solver algorithm in Ichimura et al. (2016), which can efficiently solve largescale linear systems of equations by the finite-element method in a CPU-based, large-scale computing environment.We implemented this multi-grid solver in GPU environments, since their performance has improved remarkably in recent years, making them useful for computer models demanding high computational cost.The single-precision arithmetic used in this multi-grid solver  to reduce the data transfer and computational cost is also effective in GPU environments, enabling high-speed analysis that exploits GPU performance (Yamaguchi et al. 2018).GPU computation was achieved by OpenACC, which enabled the porting of CPU codes to GPU environments by adding a small number of instructions called directives to the CPU code without using GPU-specific languages such as CUDA.Here, 8 compute-nodes, each equipped with a 2.80-GHz 24-core AMD EPYC 7402P CPU and four NVIDIA A100 GPUs, were used for computation (i.e., 32 GPUs in parallel with 32 MPI processes).By using appropriate computation algorithms implemented on GPUs, it was possible to reduce the computation time per case to about 4300 s.

Results
The horizontal coseismic displacement at the surface roughly coincided in direction and scale with the coseismic slip (Fig. 1), with large displacements occurring near the trough axis (Fig. 4a).Coseismic displacements typical for reverse-fault-type earthquakes occurred, with subsidence at the lower end of the slip zone and uplift at the upper end (Fig. 4b).
Figure 5 shows the postseismic horizontal displacement rate at the surface with respect to a reference point near Oki Island (36° N, 133° E), where the GNSS station farthest from the Trough axis is located.At sufficiently low LAB viscosities, large landward deformation was seen offshore despite small onshore deformation (case 4).At sufficiently low oceanic mantle viscosities, landward displacement on the offshore side was large, and onshore seaward deformation was larger than those of other cases (case 5).
Figure 6 shows the postseismic displacement rate in the vertical direction.As in the horizontal results, large deformations appeared at sufficiently low LAB viscosities (case 4) or at sufficiently low oceanic mantle viscosities (case 5).In the first 2 years after the earthquake, surface deformation near the trough axis varied greatly depending on the viscosity, with some subsidence observed in cases 3 and 4. On the other hand, in the third and fourth years, surface deformation near the trough axis became less viscosity-dependent, with vertical rises appearing regardless of the LAB viscosity.This indicates that the earlier state of postseismic stress in the LAB contributes to subsidence near the trough axis due to viscous relaxation.In addition, rapid stress relaxation in the lowviscosity LAB resulted in the completion of the displacement rate decay and no longer caused subsidence in the third and fourth years.On land, the displacement magnitude varied with time and case, but there was no significant trend in the uplift-subsidence pattern from north to south.
Figure 7 shows the profile along line AB (see Fig. 1) of horizontal displacement per year, with displacement toward the offshore side plotted as positive.Large, positive coseismic horizontal displacement (up to 20 m) typically occurred on the hanging wall side for reversefault-type earthquakes (Fig. 7a).After the earthquake, a displacement peak appeared around the d = 270 km point in all cases except case 4. Focusing on the movement of the hinge position with time, where the direction of displacement changes, case 1 showed almost no movement of the hinge position.In cases 2 and 3, the hinge moved about 20 km landward over 4 years.In case 5, the hinge moved about 10 km seaward over 4 years.When the viscosity is not uniform (cases 2-4), differential stress relaxation in the viscous layers and differential displacement rate decay at the ground surface result in a shift of the hinge position.On the other hand, the uniform viscosity causes the almost uniform damping of stresses and displacement rates, resulting in little shift of the hinge position (cases 1 and 5).The offshore displacement initially took on a very large negative value of −40 cm/ yr in the case of low viscosity (cases 4 and 5), but rapidly decayed over time.In case 3, due to the moderately low viscosity of the LAB, the landward displacement was 14 cm/yr at maximum, which was slightly larger than the 10 cm/yr of case 2 without a low-viscosity LAB.When the viscosity of the LAB layer is moderately low, the viscous relaxation proceeded rapidly, enhancing the offshore displacement rate due to stress relaxation in the early stages.The stress accumulation in the other layers and associated displacement rate decay would have resulted in a displacement rate close to case 2 after 4 years.In case 1, the deformation occurred at a constant rate over 4 yr with little change in the distribution of the displacement.
Figure 8 shows the profile along line AB (see Fig. 1) of vertical displacement per year.Unlike the horizontal displacement, the vertical coseismic displacement exhibited a short-wavelength component in the offshore area because the vertical component is sensitive not only to the slip distribution, but also to plate boundary surface shape (Fig. 8a).The displacement showed a longwavelength pattern that was symmetric with respect to the coastline ( d = −200 km) for the case of uniform viscosity (cases 1 and 5) but asymmetric for the case of non-uniform viscosity (cases 2-4).For the onshore area, the displacement pattern in all cases showed a peak with subsidence near the coastline and an inland hinge where the subsidence transitioned into an uplift.In cases 1 and 5, uplift was observed over a wide area ( −500 km ≤ d ≤ −280 km), while in cases 2-4 no uplift was observed within d ≤ −400 km.The peak of subsid- ence remained near the shoreline regardless of viscosity or time, while the hinges moved landward with time in all cases, especially in case 4 the hinge moved as much as 30 km (Fig. 8c).For the offshore area, in cases 1 and 5, the peak of uplift was observed around d = −100 km, and the amount of offshore uplift was almost the same as on the onshore side.On the other hand, in cases 2-4, the uplift and subsidence trend differed from those of the onshore area.For sufficiently low LAB viscosities, the trough axis and its seaward side ( d ≥ 0 km) tended to subside (cases 3 and 4).For very low LAB viscosities, offshore uplift exceeded onshore uplift (case 4).
Figure 9 shows the time variation in the horizontal ground surface displacement rate at the six observation points shown in Fig. 1.The presence of a layer with lowviscosity caused a large difference in displacement rate immediately after the earthquake that decayed rapidly over time (cases 4 and 5).In particular, in case 4, very large landward displacement rates of more than 50 cm/ yr occurred immediately after the earthquake at points 4, 5, and 6, but over the course of about 200 days, viscous relaxation in the LAB progressed rapidly, and the displacement rate decayed significantly to 30 cm/yr.Thus, onshore point 4, which was very close to the offshore area, and points 5 and 6 on the seafloor provided strong indications regarding the viscosity of the LAB.Comparing cases 2 and 3, the difference in displacement rate immediately after the earthquake at points 4, 5, and 6 was about 5 cm/yr at maximum, suggesting no significant effect of moderate viscosity in the LAB.In case 5, in which the viscosity of the whole oceanic mantle was low, large seaward displacement rates were observed, especially at onshore points 1, 2, and 3.In case 1, in which the overall mantle viscosity was high, little change over time was apparent in displacement rate.
Figure 10 shows the time history of vertical displacement rate at the observation points.As with the horizontal displacement rate, a large displacement rate immediately after the earthquake decaying over time, was observed in case 4, in which the viscosity of the LAB was very low.However, the displacement rate in the vertical direction was at most about 30 cm/yr, which is smaller than in the horizontal direction.On the other hand, the vertical component showed behavior that differed significantly between neighboring observation points 2 and 3 and between 5 and 6, where horizontal displacement direction and decay over time were almost identical.Specifically, in case 5, uplift occurred at observation point 3 but subsidence occurred at observation point 4. In addition, only a small uplift of 5 cm/yr occurred at observation point 5, while a relatively large uplift of 12 cm/yr occurred at observation point 6.In case 4, in which the viscosity of the LAB was very low, the displacement rate varied significantly with time at observation points 3, 5, and 6, leading to the detection of a low-viscosity LAB.In particular, at observation point 6, the direction of the displacement rate reversed over time, giving a strong indication of the low-viscosity LAB.
Figure 11 shows the horizontal and vertical displacement in the vertical section along line AB.Here, we explain the mechanism of surface deformation as shown in Figs. 5, 6, 7, 8, 9 and 10, with viscous relaxation and deformation patterns reflecting each of the internal structures.In all cases, a large downward flow occurred around d = −200 km with converging flow on both flanks, but for case 4, the downward flow was tilted in the landward direction.In case 2, the viscosity As a result, the displacement in the elastic slab above the LAB was large, and the relative displacement elsewhere was very small.Thus, large horizontal deformations on the offshore side were observed at the ground surface.In case 4, the displacement rate rapidly decayed significantly over time.In particular, the vertical component attenuated significantly due to the gravity, so that the displacements observed in the elastic slab were almost entirely horizontal in 2-4 years.In case 5, the displacement was very large due to the low-viscosity of the entire viscoelastic mantle resulting in quite large displacement on the onshore side, but the displacement pattern was similar to that in case 1.
Figure 12 shows dd-, zz-, and dz-components of the instantaneous deviatoric stress s (0) in the viscoelastic layer and the average viscoelastic stress glut tensor rate during the first year ṁ(0−1) in the vertical section along line AB.Note that color scales are different among cases.Here, we define stress glut tensor rate ṁ as: where τ = η/µ and s ij = σ ij − σ kk δ ij /3 are the relaxa- tion time and deviatoric stress, respectively.Superscripts (4) ṁ = −s/τ , ( (0) ) and ( (0−1) ) indicate the instantaneous state and aver- age during the first year.Reorganization of the Maxwell model (Eqs.( 2) and ( 3)) allows us to interpret viscoelastic relaxation as standard elastic deformation driven by the equivalent force rate as follows: where The strain rate ǫ and σ ′ ij are given by the linear elastic relation the body force rate ḟi is equivalent to the tensor derivative of stress glut tensor rate.The stress glut tensor rate is reciprocally proportional to the viscosity; thus, it increases when the viscosity is low (e.g., case 4).In case 4, the average stress glut tensor rate during the first year exceeds the instantaneous deviatoric stress.From Eq. 8, this suggests that σ ′ are mainly counter to the stress glut tensor rate and that the effective relaxation time is longer than the relaxation time (0.11 yr in the LAB).We should integrate the response from stress gluts over a wide region, to quantitatively evaluate displacement rates (Barbot and Fialko 2010).However, when the stress source is locally concentrated, a qualitative explanation is available.Based on the stress glut tensor distribution (Fig. 12), the displacement rate described above can be explained as follows: (5)  12b-1).Dilation sources are also observed in the deeper Oceanic mantle that is not indicated in the figures.In cases 1, 2 and 5, vertical integration of m zz results in column dila- tion, which led to an uplift at the surface directly above.On the other hand, when the viscosity of LAB was low (cases 3 and 4), the compression component in the shallow part dominates (e.g., see box iii in Fig. 12b-3) and subsidence occurred at the surface (e.g., Fig. 8 case 4).
• In case 4, the relative magnitude of ṁ(0−1) zz in the LAB was so large (Fig. 12b-4 iv) that the dilation and compression pattern in the LAB almost controls the surface uplift and subsidence pattern directly above (Fig. 13).
• In cases 3 and 4, the shear component of stress glut tensor rate ṁ(0−1) dz was negative in most of the LAB (e.g., see box v in Fig. 12b-4), which contributed to the landward displacement rate above the LAB and to the seaward displacement rate below the LAB (e.g., Fig. 11 case 4).

Discussion
We now discuss the possibility of detecting among-case differences in displacement rate associated with the modeled differences in LAB viscosity (Table 3, Figs. 9 and 10).The ∼ 35 mm/yr maximum difference in cases 2 and 3 was significantly larger than the observation error, assuming the displacement rate error based on the land GNSS measurements to be 4 mm/yr for the horizontal component.For the vertical component, the maximum difference in vertical displacement rate between cases 2 and 3 was 9 mm/yr.This difference is similar in magnitude to the observation error, which can be assumed to be 8 mm/yr (twice the horizontal component) and thus is difficult to detect with any confidence.These thresholds were assumed with reference to values from previous studies of transient crustal deformation (e.g., Ochi 2015;Iinuma et al. 2016;Ozawa et al. 2016).The difference in vertical displacement rate between cases 3 and 4 was 40 mm/yr at maximum, which was substantially larger than the observational error and thus offered a strong possibility of detection.
The accuracy of horizontal displacement rate based on seafloor GNSS-A measurements can be assumed to be 60 mm/yr (Yokota et al. 2021).The maximum difference in displacement rate associated with the LAB viscosity difference between cases 2 and 3 was up to 50 mm/ yr (SIOW: observation point 7 in Fig. 1) and thus difficult to detect with any confidence.On the other hand, the maximum difference between cases 2 and 4 was 350 mm/yr (SIOW), which was sufficiently large, compared to the observational accuracy, to be detectable.The  observational accuracy of GNSS-A measurements for the vertical component can be assumed to be 120 mm/ yr (again, twice the horizontal component).The maximum difference in displacement rate due to the difference in viscosity of the LAB between cases 3 and 4 was 180 mm/yr (TOK3: observation point 8 in Fig. 1), and thus detectable.On the other hand, the maximum difference between cases 2 and 3 was 28 mm/yr (TOK3) and thus difficult to detect.The accuracy of the vertical displacement rate as estimated from OBP observations at DONET stations likely reached 10 mm/yr (Machida et al. 2020), making even the difference between cases 2 and 3 likely detectable.
Even when the calculated surface displacement rate differences due to rheological differences exceed the observational error, a low-viscosity LAB was not always detectable because the effects of the afterslip and viscous relaxation needed to be estimated simultaneously.Afterslip mainly occurs on the plate interface deep below the coseismic rupture, namely 20-45 km in depth, and thus most of the land GNSS stations were located directly above this area.As for observations of the horizontal displacement rate at onshore GNSS station, most of these displacements were oriented toward the trough axis side, in keeping with the direction of the displacement due to the afterslip.For this reason, land GNSS data alone are not likely to be capable of distinguishing between displacement rate differences associated with afterslip differences from those associated with LAB viscosity differences.However, some land GNSS stations (Fig. 9(4)) showed large landward displacement rates; these may be effective in constraining the viscosity of the LAB.Conversely, the determination of the viscosity may change the estimation of the afterslip.On the other hand, the seaward side showed significant among-case differences in deformation pattern due to LAB viscosity differences, and, crucially, the effects of low LAB viscosity and of the displacement rate due to afterslip at the GNSS-A stations were in opposite directions.Therefore, at very low LAB viscosities (case 4), crustal deformation exceeded the observational accuracy of GNSS-A, offering a strong possibility that even the current seafloor observation networks can obtain observational data uniquely identifying signals of low LAB viscosity.For LAB viscosities that were not extremely low (case 3), the calculated displacement rates at the GNSS-A stations were comparable in magnitude to the observational accuracy of the horizontal component, making a low-viscosity LAB undetectable.However, it is possible to improve the observational accuracy of the displacement rate at the GNSS-A stations by increasing the observation frequency.Therefore, conducting more frequent measurements immediately after a megathrust earthquake would be effective to detect a low-viscosity LAB, since the difference in displacement due to the difference in viscosity of the LAB would be large.In addition, the density of the current network of observation points is low given the complexity of deformation pattern on the seaward side up to 1 year after the earthquake.Therefore, installation of additional GNSS-A stations may be effective in constraining the viscosity of the LAB based on the spatial pattern of postseismic deformation.Furthermore, it is highly likely that a lowviscosity LAB in case 3 would be detected given an observational accuracy of 10 mm/yr for the vertical displacement rates estimated from the OBP data observed at DONET stations.
However, in practice, it is necessary to estimate distributions of the afterslip and coseismic slip, and moreover it is necessary to determine how strong a constraint can be applied to the viscosity of the LAB given these factors.Therefore, it is necessary to discuss the sensitivity of observations at observation points in   Johnson and Tebo 2018).The contribution of afterslip to the postseismic vertical deformation between 0 and 100 km from the trough axis is suggested to be small (Johnson and Tebo 2018).Therefore, our discussion on the detectability of low viscosity in the LAB still holds, even when afterslip is taken into account.On the other hand, the contribution between 100 and 200 km from the trough axis is suggested to be large, and the calculated displacement rate may be significantly different from the actual possible displacement rate.Since viscous flow is initially conditioned by the coseismic slip distribution, it is impossible to discuss the details of displacement rates for different settings of slip distributions from our results.Therefore, the response to past Nankai earthquakes and other possible earthquakes is a future issue.On the other hand, our results shows the general patterns of the viscous flow considering the low-viscosity LAB: a large downdraft occurs near the bottom of the slip (Fig. 11), where the hinge is located and the landward displacement rate is enhanced by the low-viscosity LAB on the seaward side.This result is similar to the postseismic displacement of the Tohoku-Oki earthquake, and such a trend is expected for different slip distribution.
The present study focused on the deformation pattern from the trough axis to the coast, disregarding the influence of the elastic thickness of the upper plate.Although greater detail regarding crustal structure can be obtained by incorporating the results of seismic wave tomography and structural exploration, viscoelastic deformation is dominated by viscous parameters and thickness of elastic plate.The influence of detailed elastic structure can thus be considered low.The thickness of the subducting plate is 30 km, whereas the thickness of the LAB is as thin as 10 km.LAB thickness trades off with its viscosity ( Freed et al. 2017).Setting of LAB should be refined by constraints from other observations.The LAB and subducting plates are only modeled up to the depth of 75 km as the geometry of their deeper part is unclear.Even if we model the slab and LAB extending into deeper sections, the change in stress state due to coseismic slip is expected to be small, and the effect of viscoelastic relaxation on surface displacement in the corresponding sections is also expected to be small.
In this study, we performed crustal deformation analysis considering Maxwell rheology, and showed that the use of crustal deformation measurements obtained during the first year after an earthquake are effective for detecting a low-viscosity LAB.Sherrill and Johnson (2021) has shown that the transient viscosity coefficient does not improve the results for crustal deformation data after 1 year from the earthquake.On the other hand, previous studies have shown that fast viscoelastic deformation in the very early postseismic period (a few months) should be modeled by Burgers rheology, which allows for a transient component (e.g., Qiu et al. 2018).In the case of the Tohoku-oki earthquake, LAB, mantle wedge, and oceanic asthenosphere are modeled by bi-viscous viscoelastic rheology, respectively, to explain the observed data, while noting that simple bi-viscous rheology does not explain all the observed points (e.g., Sun et al. (2014)).Other authors have suggested that nonlinear viscoelasticity modeling can explain fast landward crustal deformation in the seaward area due to a localized effective viscosity decrease caused by coseismic stress concentration (Agata et al. 2019).Therefore, future analysis should incorporate more complex rheologies.Further research on the mechanism generating low viscosity at the base of the subducting slab also remains desirable.

Conclusions
In this study, we used a GPU-based viscoelastic crustal deformation analysis method that significantly reduced computational cost to calculate viscoelastic deformation around the Nankai Trough.Using a high-resolution 3D model that reflects the realistic crustal structure, we examined the effects of hypothesized differences in viscous structure on deformation patterns at the ground and seafloor surfaces.The surface response to coseismic slip for 4 years after an earthquake was calculated for cases in which various viscosities are assigned to the mantle wedge, oceanic mantle and LAB.The results show that for very low LAB viscosities ( 2.5 × 10 17 Pa s), very large landward displacement occurs on the offshore side, while in the vertical direction, complex behavior is predicted with varying distributions of uplift and subsidence in different cases.Especially in the first year after the earthquake, the presence of a low-viscosity LAB produces a magnitude of deformation that is easily distinguishable from the case without a low-viscosity LAB, given the observation error levels of the seafloor stations.Hence, if an LAB of very low-viscosity exists, the probability of detecting it is high during the first year.On the other hand, for moderately low LAB viscosities, e.g., 2.5 × 10 18 Pa s, the deformation pattern is similar to the case without a low-viscosity LAB.Since the difference in such cases is less than the observational accuracy levels available in current observation networks, it is expected to be difficult to detect an LAB of moderately low-viscosity, if indeed it exists.However, such low-viscosity LAB may be detected by improving observational accuracy in the displacement rate at GNSS-A stations by increasing the frequency with which observations are made and by achieving expected improvements in observational accuracy in the vertical displacement rate as estimated from OBP observations at DONET stations.It should be noted that the surface response was calculated assuming the maximum rupture in the expected source area, but a partial rupture or a two-stage rupture are also possible.In these cases, the response may therefore differ from the expectation set herein depending on the actual rheological structure and its more complex distribution.Future work should investigate the detectability of a low-viscosity LAB by examining response behavior in the context of these alternative rheological and earthquake scenarios.

Fig. 1
Fig. 1 Reference coseismic slip.Black contours represent the trough axis and the depth to the subducting Philippine Sea Plate at intervals of 10 km.Black circles, squares, and triangles represent onshore GNSS, GNSS-A (JCG) and GNSS-A (NAGOYA) stations, respectively.The red circle indicates the reference point.Line AB is the section profiled in Figs.7 and 8. Blue numbered points are the positions of the observation points shown in Figs. 9, 10 (points 1-6) andTable 3 (points 1-8)

Fig. 2
Fig. 2 Rheological structure of the vertical section along the line AB in Fig. 1.The elastic structure is based on the proposed Japan integrated velocity structure model (JIVSM) version 1(Koketsu et al. 2009(Koketsu et al. , 2012)), and the rheological structure and elastic structure were combined to create a crustal structural model.Numbered points (1 to 6) are the approximate locations of the observation points shown in Fig.1

Fig. 7
Fig. 7 Profiles along line AB (see Fig. 1) of a horizontal coseismic displacement and b horizontal displacement rate, with displacement toward the offshore side plotted as positive.The solid line, dashed line dot-dashed line, and dotted line indicate displacement rates in years 1, 2, 3, and 4, respectively.The purple line in b2-b5 shows the displacement rate 1 year after the earthquake in case 1.The d = −200 km point represents the coastline

Fig. 8
Fig. 8 Profiles along line AB (see Fig. 1) of a vertical coseismic displacement and b vertical displacement rate for each case.The solid line, dashed line, dot-dashed line, and dotted line indicate displacement rates in years 1, 2, 3, and 4, respectively.The purple line in b2-b5 shows the displacement rate for 1 year after the earthquake in case 1.The d = −200 km point represents the coastline

Fig. 9
Fig. 9 Time history of horizontal displacement rate of the ground surface at observation points (see Fig. 1 for the locations of observation points)

Fig. 10
Fig. 10 Time history of vertical displacement rate of the ground surface at observation points (see Fig. 1 for the location of observation points)

Fig. 12 a
Fig. 12 a Instantaneous deviatoric stress and b stress glut tensor rate in the vertical section along line AB (shown in Fig. 1).Color scales are different among cases

Fig. 13 a
Fig. 13 a Profiles along line AB (see Fig. 1) of vertical displacement rate (same as Fig. 8b-4) and b vertical dilation component of stress glut tensor at the middle of the LAB in the vertical section along line AB during the first year in case 4

Table 1
Material properties ( V p , V s , ρ and η indicate primary wave velocity, secondary wave velocity, density, and viscosity, respectively)

Table 3
Among-case differences in displacement rate (mm/ year) associated with the modeled differences in LAB viscosity Location of observation points are shown in Fig.1