Probabilistic characterization of seismic ground deformation due to tectonic fault movements

Tectonic ground deformations in the near-fault region cause major damage to buildings and infrastructure. To characterize ground deformation demands on structures, a novel stochastic approach to evaluate the ground deformations of tectonic origin is developed by combining probabilistic models of earthquake source parameters, synthetic earthquake slip models, and Okada equations for calculating the deformation field due to a fault movement. The output of the method is the probability distribution of ground deformations at a single location or differential ground deformations between two locations. The derived probabilistic models can be employed as input to advanced structural models and analyses. The method is illustrated for the 16 April 2016 Kumamoto earthquake in Japan. By comparing simulated ground deformations with observed deformations at multiple sites, a set of refined source models is first derived and then used to investigate the detailed earthquake characteristics of the event and to develop probability distributions of tectonic ground deformations at target


Introduction
Movement of a fault triggers various effects on the ground, such as strong shaking and permanent deformation, and causes major destruction to buildings and infrastructure. In a near-fault region, both strong shaking and permanent ground deformation are strongly influenced by the characteristics of the fault movement, such as type of faulting (strike-slip, normal, and reverse) and slip distribution of earthquake rupture. The near-fault ground motions generate so-called killer-pulses [1] and are particularly damaging to structures [2]. Whereas careful seismic design considerations are needed when linear civil structures are subjected to large tectonic deformations [3][4][5]. On occasion, fault rupture reaches the ground surface and differential movement of the ground can be in the order of several meters to several tens of meters, depending on the earthquake characteristics, local site conditions, and various other factors, including earthquake-triggered settlement and sliding of soil [6,7]. Generally, the effects of the fault movement appear over extended areas along a fault with distributed ground deformation [4,8].
Predicting seismic ground deformations due to an earthquake has been conventionally carried out using empirical ground motion models [9]. Empirical ground motion models evaluate the amplitude of peak ground displacement as a function of earthquake magnitude, source-tosite distance, and site condition. Uncertainty of the estimation is characterized by a so-called logarithmic standard deviation (or sigma).
However, such conventional methods are faced with two difficulties. Firstly, in the near-fault region, estimating permanent ground deformations from observed accelerograms is not a simple task; careful baseline correction of the acceleration data is of critical importance [10,11]. Secondly, observed ground motion data are scarce in the nearfault region and the number of ground motion data that are usable for developing empirical ground motion models decreases drastically with longer spectral periods of ground motion intensity measures [12]. Moreover, empirical ground motion models cannot be applied directly to evaluate ground deformations at multiple nearby locations simultaneously, because spatial correlation models of ground deformations are not usually available [13]. In other words, conventional prediction models have been derived for a single site and spatial dependence of the seismic effects, which is important in evaluating differential ground deformations at two locations, is not fully characterized.
In light of preceding limitations of current tools for predicting spatially distributed ground deformation fields, it is necessary to develop an alternative approach. One way is to implement so-called Okada equations [14] to analytically compute ground deformations due to tectonic fault movements in an elastic half-space. Despite simplicity, Okada equations have been employed in various fields, including geodesy, volcanology, and tsunami modeling [15,16]. An important advantage of using Okada formulae is that the spatial field of ground deformations can be simulated and thus can be used for evaluating differential ground deformations between two locations. Input information of Okada equations is an earthquake rupture model that is defined in terms of source parameters, such as fault geometry and slip distribution. For a given earthquake scenario (e.g. earthquake magnitude and faulting type), source parameters can be estimated using empirical scaling relationships [17][18][19][20][21][22]. To capture the uncertainty of earthquake rupture, a stochastic synthesis method of earthquake slip distribution [23][24][25] can be adopted. It is highlighted that accounting for variability of earthquake slip is particularly important for evaluating near-fault ground motions and deformations. The requirements for capturing uncertainties of both geometrical parameters and earthquake slip distributions narrow the choice of applicable scaling relationships into those by [22]. This is because other equations are mainly focused upon geometrical parameters only, while the relationships developed by [22] are applicable to earthquake source parameters that specify the key features of slip heterogeneity over a fault plane.
This study develops a novel approach to characterize the ground deformations of tectonic origin probabilistically. Key components of the integrated method include: (i) probabilistic models of earthquake source parameters [22], (ii) stochastic synthesis of earthquake slip [25], and (iii) Okada equations [14,16] for calculating the deformation field due to a fault movement. The output of the method is the probability distribution of ground deformations at a single location or the probability distribution of differential ground deformations between two locations. The derived probabilistic models can be employed as input to advanced structural models and analyses (e.g. reliability analysis of a bridge subjected to differential ground deformations). In this study, the method is illustrated by focusing on the 16 April 2016 Kumamoto earthquake in Japan. This earthquake registered the moment magnitude (M w ) of 7.0 and was of a right-lateral strike-slip type. It struck rural areas of Kumamoto Prefecture in Kyushu Island of Japan and caused major destruction to buildings and infrastructure in the nearfault region [26]. Notably, the surface rupture was observed at many locations with horizontal dislocations exceeding 2 m at several locations [8,27]. Furthermore, for this event, several inverted source models have been developed [28,29], and geophysical observations, such as strong motion data [26], GPS deformation [30], and InSAR satellite images [8,31], are available. In the main part of the results of this study, numerous stochastic slip distributions having heterogeneous slips are generated based on the scaling relationships and stochastic synthesis method to identify a set of earthquake slip distributions with a better fit with near-fault GPS observations and strong motion data than the uniform-slip source model developed by the Geospatial Institute of Japan (GSI). The refined set of the source models can be used to develop the probability distribution of tectonic ground deformations at sites of interest. Finally, in a case-study application, sites near the Aso Bridge, which had collapsed due to the Kumamoto earthquake, will be focused upon, noting that the Aso Bridge is very near (within 0.1 km) from the fault rupture observed in the post-earthquake fieldwork [26,27] and from the fault strike of the causative Futagawa Fault. The illustration of the methodology will be followed by key conclusions from this study.
2. Estimating tectonic ground deformation through stochastic earthquake source modeling

Methodology
A procedure to evaluate the probabilistic characteristics of tectonic ground deformations due to an earthquake is developed in this section. The method is comprised of four steps: (i) definition of an earthquake scenario, (ii) stochastic synthesis of earthquake source models for the specified scenario, (iii) estimation of ground deformations at locations of interest, and (iv) post-processing of ground deformation results. A computational procedure is based on Monte Carlo simulation and is illustrated in Fig. 1. The main idea of the method is straightforward. For a specified target scenario (i.e. step (i)), numerous earthquake sources are generated stochastically based on seismological theories and models (i.e. step (ii)) and corresponding ground deformation fields are evaluated using Okada formulae (i.e. step (iii)). By adopting a subset of synthesized source models and deformation fields according to some acceptance criteria, probabilistic distributions of ground deformations can be characterized (i.e. step (iv)).
It is important to highlight that two types of investigations can be carried out using the stochastic procedure for estimating tectonic ground deformations outlined in Fig. 1. One way is to apply to a general forecasting situation, while the other is to apply to a retrospective case study. The former is applicable to situations where general seismological/geological information is available to define an earthquake scenario of interest but without specific observations. In this case, simulated stochastic sources and resulting ground deformations (i.e. steps (ii) and (iii)) are constrained on existing geological and seismological models (e.g. fault geometry based on regional geological and geomorphological investigations and scaling relationships of earthquake source parameters). For the latter case, when geophysical (field) observations are available for a specific earthquake, predictions based on the stochastic source models can be refined by comparing model predictions with the observations; source models that produce poor match with the observations may be discarded, while many more models that are in closer agreement with the observations can be generated. The selection criterion can be determined based on existing (crude/preliminary) benchmark earthquake source models and some other seismological and geological information. Essentially, the latter approach is equivalent to seismic source inversion [32]. From the perspective of inverse problems, the proposed method attempts to find a set of refined source models, rather than a single best model [33,34], that produce better predictions of geophysical quantities of interest. A search mechanism of the proposed method differs from conventional gradient-based methods

Stochastic source modeling and ground deformation modeling
The first step of the developed stochastic source and deformation modeling approach is to define a moment magnitude and source region model of the target earthquake scenario. The source region should be sufficiently large such that stochastic sources having various sizes and locations can be generated over the fault plane. To account for uncertainty of fault plane geometry, multiple fault planes, e.g. having different strike and dip angles, can be defined. The defined fault plane model essentially reflects the seismological knowledge of earthquake rupture in the target region. Subsequently, the whole source region is discretized into many sub-faults. For instance, a source region for a M w 6.5 event can be determined as a rectangular fault plane of 50 km long and 20 km wide, aligned with a target geological fault, and then 25×10 sub-fault grids can be defined by considering a sub-fault size of 2 km × 2 km. When uncertainty of the fault plane geometry is taken into account, one fault plane can be chosen randomly from numerous candidate fault planes.
Secondly, stochastic synthesis of earthquake rupture is conducted in two phases. In the first phase, eight earthquake source parameters, i.e. fault length (L), fault width (W), mean slip (D a ), maximum slip (D m ), Box-Cox parameter (λ), along-strike correlation length (C L ), along-dip correlation length (C W ), and Hurst number (H), are sampled from the probabilistic scaling relationships developed by [22]. These relationships have been obtained on the basis of 226 inverted source models in the SRCMOD database [34].
Physical interpretations of L, W, D a , and D m are straightforward, whereas those for λ, CL L , CL W , and H need to be explained. The Box-Cox parameter λ is one of the slip statistics parameters and together with D a and D m , determines the marginal distribution of earthquake slip over a fault plane. It is used to capture the positive skewness of earthquake slip values over a fault plane [25,36], which are typical features of the majority of existing earthquake source models [22]. More specifically, this parameter is used to transform normally distributed slip X into nonnormally distributed slip Y as follows: When λ = 0, the Box-Cox transformation corresponds to the logarithmic transformation. In the random field generation method that is implemented in this study, synthesized slip distributions have slip values that are normally distributed. Moreover, three source parameters, i.e. CL L , CL W , and H, define the spatial heterogeneity of earthquake slip over a fault plane, whose spatial features are specified by the von Kármán wavenumber spectra [23]. The von Kármán spectra can be expressed as: ) 0.5 , and k L and k W are the wavenumbers along strike and dip, respectively (note: wavenumber is proportional to the reciprocal of wavelength). Typically, wavenumber spectra have some plateau spectral levels in the low wavenumber range and decay gradually as the wavenumber increases. The correlation lengths CL L and CL W determine the absolute levels of the power spectrum in the low wavenumber range along the strike and dip directions, respectively, and capture the anisotropic spectral features of the slip distribution (when CL L and CL W are different). The Hurst number controls the slope of the power spectral decay in the high wavenumber range, and is theoretically constrained to range between 0 and 1. Further details of the von Kármán wavenumber spectra for characterizing the spatial distribution of earthquake slip can be found in [23].
For crustal earthquakes having M w < 7.5 and normal/strike-slip faulting mechanisms, the following equations can be used to sample values of L, W, D a , D m , CL L , and CL W : where ε terms in Eqs. (3)-(8) are the prediction errors and are modeled as standard normal variable (i.e. unit mean and zero variance). Because the above six source parameters are physically interrelated, the ε terms are correlated. The correlation structure of the error terms can be characterized by multivariate standard normal distribution with the linear correlation coefficient matrix equal to: On the other hand, the Box-Cox parameter λ is characterized by a normal variable with mean equal to 0.312 and standard deviation equal to 0.278. The Hurst number H is characterized by a bimodal random variable that takes a value of 0.99 with probability of 0.43 or a value sampled from the normal distribution with mean equal to 0.714 and standard deviation equal to 0.172 with probability of 0.57 [22]. Note that λ and H are not dependent on M w and are treated as independent variables.
In the second phase of the stochastic synthesis, once values of the eight source parameters are sampled, geometry of the fault plane can be fixed (based on simulated values of L and W as well as randomly chosen fault plane) and then it is positioned within the target source region (note: simulated values of L and W are usually less than the dimensions of the source region). Subsequently, a random slip field is generated using a Fourier integral method [37] from the simulated wavenumber spectra that are specified by the spatial slip distribution parameters CL L , CL W , and H. To achieve slip distribution with realistic positive skewness, the synthesized slip distribution is converted via Box-Cox transformation using the simulated value of λ. The transformed slip distribution is then adjusted to achieve the target mean and maximum slips D a and D m . The simulated slip values can be mapped to sub-faults on the fault plane. At this stage, consistency among the simulated values of L, W, and D a can be checked by comparing the target seismic moment (as specified by the scenario magnitude) and the simulated seismic moment (note: seismic moment = μD a LW, where μ is the rock rigidity). Due to the variability in L, W, and D a and adjustments for D m , random sampling of L, W, and D a may result in a seismic moment that is different from the target seismic moment or moment magnitude. In case a simulated source model does not achieve the consistency in terms of seismic moment, the generated source model is abandoned and a new source model is generated. In this study, the target moment magnitudes minus/plus 0.05 units are considered as a criterion for seismic moment consistency.
Thirdly, for a given acceptable source model, Okada equations are used to obtain surface deformations on the ground caused by an earthquake rupture of a finite-fault rectangular source. The Okada equations analytically compute tectonic ground deformations due to shear and tensile fault movements in a homogeneous semi-infinite medium (i.e. Earth), and have been widely used in various fields (including geophysics and volcanology). The slip vector is specified by the slip and rake angle. The results are typically obtained as EW, NS, and UD deformations. When the earthquake source model consists of multiple sub-faults with assigned (different) slip values, ground deformations due to each sub-fault rupture can be estimated by applying Okada equations individually and then are summed to obtain the total ground deformation. It is noteworthy that the calculated ground deformations are for tectonic movements of the fault and do not account for earthquake-triggered settlement and sliding of surface soil.
It is important to clarify that estimated ground deformations based on Okada equations are deterministic predictions of ground displacement due to a given earthquake rupture and those at different locations are physically consistent, unlike empirical ground motion models for ground displacement. In the developed stochastic source and deformation modeling method, uncertainty associated with predictions of tectonic fault rupture is characterized and captured in the stochastic earthquake source modeling phase.

2016 Kumamoto earthquake
This section provides a brief summary of the 16 April 2016 M w 7.0 Kumamoto earthquake. The information presented in this section will be used to illustrate the developed methodology for generating a set of refined source models, which are calibrated with measured ground deformation data, and for characterizing tectonic ground deformations due to earthquake rupture.

Earthquake characteristics
A moderate-size earthquake M w 6.2, originated from the northern segment of the Hinagu Fault, struck the Kumamoto region of Kyushu Island, Japan on 14 April 2016. This event was followed by a larger M w 7.0 earthquake on 16 April 2016 occurred along the Futagawa Fault (NE of the Hinagu Fault). Both earthquakes were of right-lateral strikeslip type occurring at shallow depths. The crustal deformation due to the mainshock was manifested as ground surface ruptures at many locations along the Futagawa Fault. At several places, ground deformations up to 2 m were reported [8,27]. Important aspects of the Kumamoto earthquakes are that various kinds of geophysical observations, such as strong ground motion [26], GPS deformation [30], and InSAR satellite images [8,31], are obtained and the results are made available publicly, facilitating detailed earthquake source studies [28,29]. Fig. 2 shows a digital elevation model (DEM) of the Kumamoto region based on GDEM2 (https://asterweb.jpl.nasa.gov/gdem.asp) together with fault geometry of two earthquake source models developed by the GSI. The GSI model 1 (solid-line rectangle) is a preliminary version, developed rapidly after the mainshock based on source inversion analysis of permanent GPS measurements and known geological information (e.g. fault geometry), whereas the GSI model 2 (three broken-line rectangles) is an updated version, developed by adding GPS measurements at temporary sites [27]. The source parameters of the GSI models 1 and 2 are summarized in Table 1. See Fig. 3 for the definitions of the fault geometry and parameters. In Fig. 2, the epicenter of the 16 April 2016 mainshock, based on the unified Japan Meteorological Agency catalogue is also included; the rupture was initiated near the SW tip of the Futagawa Fault and the rupture propagated towards NE along the fault strike. The elevation data clearly show the Aso Caldera (a horseshoe shape), and the fault rupture area of the mainshock (i.e. Futagawa Fault) stretches from the mouth of the Aso Caldera towards the Kumamoto plain areas in the approximately ENE-WSW orientation. The mainshock caused significant damage in the near-fault region along the fault strike, e.g. Mashiki Town, Nishihara Village, and Minami Aso Village [26].

Geophysical and geodetic observations
Numerous geophysical observations are available from seismological networks, e.g. K-NET and KiK-net, and geodetic networks, e.g. national GPS network GEONET; locations of 2 strong motion stations (Mashiki Town Office and the Nishihara Town Office [28]) and 6 GE-ONET stations in the near-fault region of the Kumamoto earthquake are shown in Fig. 2 with grey circle markers and blue circle markers, respectively. By integrating recorded accelerograms at the strong motion stations, reliable estimates of the permanent ground deformations can be obtained [28]. Moreover, researchers at the GSI conducted rapid GPS measurements (vertical deformations only) at 13 sites in the nearfault region immediately after the mainshock. The locations of the rapid GPS measurement sites are shown in Fig. 2 with red square markers. Fig. 4 shows the observed ground deformations at the above-mentioned locations. In the figure panels of Fig. 4, the observed deformations are indicated as vector arrows (see the right-bottom of each panel for the vector scale). It can be observed that for example, two strong motion stations were on the hanging-wall side of the fault rupture because they moved towards NE and subsided downwards, in agreement with the fault rupture, whereas the Choyo station (32.8707°N, 130.9962°E), which is a GEONET station inside the Aso Caldera (see Fig. 2) along the strike of the Futagawa Fault, was on the foot-wall side as it moved towards SW and was uplifted. This indicates that the fault strike of the Futagawa Fault runs through narrow areas surrounded by the two strong motion stations and the Choyo GEONET station. In fact, these three locations are along the inferred fault strike of the Futagawa Fault based on the previous geological investigations [38]. The strike of the Futagawa Fault almost coincides with the strike lines of the GSI models shown in Fig. 2. It is highlighted that these observed ground deformations are in agreement with deformation fields estimated from the InSAR images taken after the Kumamoto mainshock [8].
Using an available earthquake source model, ground deformation fields due to the earthquake rupture can be calculated based on Okada formulae. Such deformation fields (EW, NS, and UD) for the GSI models 1 and 2 are shown in Fig. 4a and b, respectively. See Figs. 2 and 3 and Table 1 for their source geometry and parameters. The GSI models are generally consistent with the observed deformation vectors estimated from the strong motion data and GPS data. The GSI model 1 assumes a uniform slip of 3.5 m over its fault plane (shown as a solid-line rectangular in Fig. 2). The GSI model 2 produces more complex deformation fields because it consists of three sub-faults with different geometry and slip values. It is important to note that the assumption of uniform slip over the fault rupture plane is usually too simplistic and other inverted source models for the Kumamoto earthquake produce heterogeneous earthquake slip distributions [28,29]. Moreover, uniform-slip source models tend to underestimate the fault dimensions because the earthquake moment is more or less constrained for a given event whilst a constant slip essentially reduces the fault rupture area.
More realistic earthquake source models are useful for evaluating the seismic effects on buildings and infrastructure in the near-fault region.

Case study for the 2016 Kumamoto earthquake
This section presents an application of the developed stochastic   procedure for estimating tectonic ground deformations to the 16 April 2016 Kumamoto earthquake. Section 4.1 is focused upon deriving a set of refined source models for the Kumamoto mainshock. In Section 4.2, results of the estimated ground deformations at several locations in the near-fault region based on the refined source models will be discussed.

Refined source models for the Kumamoto earthquake
In deriving a set of refined source models for the Kumamoto earthquake, the target moment magnitude is set to M w 7.0 minus/plus 0.05. The target source region is determined based on existing geological investigations by [38] as well as various post-earthquake studies [27][28][29][30] (see Section 3). The base fault plane model, which consists of two sub-fault planes, is shown in Fig. 5 (see thick black rectangles). The NE sub-fault plane essentially coincides with the inferred fault plane for the Futagawa Fault, while the SW sub-fault plane corresponds to the inferred plane for the Hinagu Fault. The fault plane parameters of the two sub-faults of the base model are listed in Fig. 5, which are in agreement with the GSI models as well as the inverted source model by [28]. The total length of the fault plane is set to 62 km (= 38 km and 24 km for the sub-fault planes 1 and 2, respectively), while its width is 30 km (same for the two sub-fault planes). Due to the differences in dip angles of the sub-fault planes 1 and 2, along-dip dimensions of the subfault planes 1 and 2, projected onto Earth's surface shown in Fig. 5, are different. It is noted that the dimensions of the base fault plane are larger than typical inverted source models for the Kumamoto earthquake (see Fig. 2); therefore, it can accommodate various source models within the target source region.
To consider a wide range of possible fault planes, in total, 567 fault plane models are defined. These variations are obtained by changing one of the geometrical parameters of the base fault plane. More specifically, for the sub-fault plane 1, 7 values of latitude of the upper-left corner, 3 values of strike angle, 3 values of dip angle, and 3 values of rake angle are considered (i.e. 189 cases), whereas for the sub-fault plane 2, 3 values of longitude of the upper-left corner are considered (see Fig. 5 for the varied parameters). Note that variations of the geometry for the sub-fault planes 1 and 2 are considered to be independent. Fig. 5 also presents varied geometry of the sub-fault planes 1 and 2; it covers a wide range of possible fault planes for the Futagawa-Hinagu Faults. Selection of the varied parameters and their values is based on inspections of various post-earthquake reconnaissance reports and studies. Many investigations clearly indicated that the mainshock rupture occurred along the Futagawa Fault. Among the geometrical parameters for the Futagawa Fault, it was indicated that the NE end of the fault rupture is more uncertain than other parameters, e.g. whether the rupture extended into the Aso Caldera [29], but the strike line of the Futagawa Fault is well constrained by the previous geological investigations [38] and also by various field observations [27] (see Fig. 2).
Subsequently, each of the fault plane models is discretized into subfaults of 2 km × 2 km (note: similar stochastic source and deformation   By implementing the stochastic source generation and slip synthesis, numerous candidate earthquake source models for the Kumamoto earthquake are generated. In the simulation, the GSI source model 2 (see Figs. 2 and 4 and Table 1) is adopted as a benchmark for deciding acceptance or rejection of a simulated candidate model. The error score of a candidate source model is calculated as follows: where u EW , u NS , and u UD are the predicted ground deformations based on the candidate source model and u Obs,EW , u Obs,NS , and u Obs,UD are the observed ground deformations. The prediction errors of ground deformations are evaluated at 8 sites for the EW and NS directions and at 21 sites for the UD direction. Sites 1 and 2 are the strong motion recording stations (Mashiki and Nishihara); sites 3-8 are the permanent GEONET stations, and sites 9-21 are for the temporary GPS measurements by the GSI researchers. In other words, a candidate model that produces a better fit with observed ground deformations at the strong motion sites and GPS measurement sites than the GSI model 2, is kept as a refined source model; the error score for the GSI model 2 is calculated as 8.41 (see Fig. 4b). In this study, the total number of refined models is set to 1000. To generate 1000 refined models, more than 86 million candidate models are synthesized and analyzed. A histogram of the total error score of the 1000 refined source models for the Kumamoto mainshock is shown in Fig. 7. The error score threshold for the top 100 models is 7.35. Figs. 8 and 9 show detailed results of two source models. Fig. 8 is for a candidate model that is considered to be unacceptable in terms of total error score (thus not included in the set of the 1000 refined models), whereas Fig. 9 is for the best source model in the set of the 1000 refined models. The error scores for the two models are 17.  slip is about 3.0 m). It predicts ground deformation fields that are consistent with the observed slips in terms of spatial distribution (Fig. 8b), however, at the 21 sites, mismatch of the predicted and observed ground deformations are significant, particularly in the nearfault region (Fig. 8c and d). This is caused by the asperity characteristics of the model, i.e. smaller slips and their concentrations at the center of the fault plane (rather than top edge along strike). In contrast, the best source model shows a more compact fault rupture area with greater maximum slip (Fig. 9a), asperity areas of which are near the top edge of the fault plane, concentrated near Nishihara Village (32.8348°N, 130.9030°E) where major subsidence was observed in the post-earthquake surveys [26,27]. The deformation fields of the best model (Fig. 9b) exhibit much greater deformations and sharper deformation boundary along the fault strike, in comparison with those for the unacceptable model (Fig. 8b). Better agreement of the predictions based on the best model, with respect to the unacceptable model as well as the GSI model 2, can be inspected in Fig. 9c and d. Similar evaluations have been carried out to derive the set of 1000 refined source models.
Advantages of generating multiple acceptable source models, rather than a single best model, are that source characteristics of the refined models can be studied in detail. Fig. 10 presents histograms of four varied geometrical parameters of the sub-fault plane 1, i.e. latitude of the upper-left corner, strike, dip, and rake, based on the 1000 refined models (see Fig. 5 11. Comparison of the simulated earthquake source parameters and the scaling relationships by [22]. The total number of accepted models is set to 1000 (blue circles). The top 100 models are shown with light blue circles. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article).

K. Goda
Soil Dynamics and Earthquake Engineering 100 (2017) 316-329 by [38]. Note that 0.005°shift corresponds to a movement of location by approximately 0.55 km. Essentially, the strike line determines whether a site is on the foot-wall side or hanging-wall side. For the Kumamoto earthquake case (i.e. right-lateral strike-slip event with normal fault components), sites on the hanging-wall side move towards NE horizontally along the fault and subside (see Fig. 4). The results shown in Fig. 10a suggest that slight shifts of the strike line (i.e. minus/ plus about 0.01°or minus/plus 1 km) are possible, whereas further shifts leading to more than 2 km are not strongly supported by the deformation data in the near-fault region. In other words, the strike line for the Kumamoto earthquake is well constrained within relatively narrow areas. This result is consistent with other studies that employed different geophysical measurements and techniques, such as InSAR images [8,31]. Fig. 10b shows that the strike angles of the 1000 refined models are either 230°or 235°(default value), indicating that more inclinations of the strike tend to generate consistent predictions of ground deformations in the near-fault region. Fig. 10c shows almost equal contributions from three dip angles, suggesting that observed deformations are not significantly influenced by the fault dip. On the other hand, Fig. 10d clearly shows superior fit of source models containing more normal components (which will increase the subsidence of the ground in the hanging-wall region). Fig. 11 compares the simulated values of the eight source parameters of the 1000 refined models with the scaling relationships by [22]. In the figure, results for the top 100 models are displayed with circle markers with light blue color. The results shown in Fig. 11 indicate that on average the simulated fault length is less than the expected length for a typical M w 7 event, while the simulated fault width is generally consistent with the empirical prediction for this earthquake magnitude (noting that significant variability of these two parameters). The shorter simulated fault length (typically 30-40 km) is related to the length of the Futagawa Fault and is consistent with the GSI models ( Fig. 2 and Table 1). Consequently, it can be concluded that the fault rupture area of the Kumamoto mainshock is relatively smaller, compared with a typical M w 7 event. Fig. 11c and d clearly show that the simulated mean and maximum slips of the 1000 refined models are greater than the empirical predictionsthese results are the consequences of the compact fault rupture areas of this event and are in agreement with other inverted source models for this event [28]. The asperities having relatively large concentered slips may be the direct cause of the significant ground deformations in the near-fault region [27]. The simulated values of the Box-Cox parameter (Fig. 11e) are similar to the assumed prediction model, indicating that the tail feature of the earthquake slip for the Kumamoto earthquake is similar to those of other historical events. The simulated values of the correlation lengths along strike and dip directions are less than those indicated by the empirical prediction models ( Fig. 11f and g). These results are expected because the correlation lengths and fault dimensions are positively correlated [22] (see Eq. (9)). On the other hand, the simulated values of the Hurst number are consistent with the assumed prediction model (Fig. 11h).
Furthermore, it can be observed from Figs. 10 and 11 that distributions of the simulated source parameters for the 1000 refined models and for the top 100 models are similar. It is important to emphasize that the effects of the earthquake source parameters (not only for fault geometry but also for slip distribution) are very complex; detailed results from the stochastic source modeling will provide useful insights regarding the complex interactions of these parameters.

Probabilistic characterization of tectonic ground deformations
The set of 1000 refined source models, derived in Section 4.1, can be considered as realistic representations of possible earthquake ruptures for the 2016 Kumamoto earthquake. On this premise, probability characteristics of ground deformations at target sites can be evaluated. In this study, two locations, i.e. Point 1 and Point 2, near the Aso Bridge, which collapsed during the mainshock, are focused upon (Fig. 12). The coordinates of Points 1 and 2 are (32.8838°N, 130.9871°E) and (32.8831°N, 130.9903°E), respectively; they are separated by 0.3 km. It is noted that the Aso Bridge (i.e. Points 1 and 2) Fig. 12. Elevation data near the Aso Bridge together with two photos from the reconnaissance work in Kumamoto [26]. are about 1.7 km NW of the Choyo GEONET station, where SW horizontal deformation and uplift are observed (i.e. the Choyo station is on the foot-wall side; see Fig. 4).
The Aso Bridge was an inverted ranger steel bridge, crossing over the Kurokawa River. Its total span length was 206 m consisting of 134 m long main structure, 54 m long three-span structure (supported by two piers) on the western side, and 18 m long single-span structure on the eastern side. The most likely cause of the collapse of the Aso Bridge was due to a massive landslide (see Photo A in Fig. 12) that affected the western section of the bridge; two bridge piers as well as the main steel structure had fallen into the river gorge. On the eastern side of the bridge, many surface fault ruptures were observed; Photo B shows the surface ground cracks in nearby areas of the bridge (about 0.1 km east). However, the exact failure mechanism/process of the bridge is unknown because the earthquake occurred at 1:25 am (local time) and there is no eyewitness. It is important to clarify the aim of the  investigations and results presented in this section. Rigorous examinations of the causes of the Aso Bridge collapse are beyond the scope of this study. Points 1 and 2 are used to illustrate how the developed stochastic method for estimating tectonic ground deformations probabilistically can be implemented. Fig. 13 shows histograms of the simulated ground deformations (EW, NS, and UD) at two target sites based on the 1000 refined source models. The results for Point 1 and Point 2 are similar due to the proximity of the two sites. It can be observed that significant ground deformations are experienced at both sides of the Aso Bridge across the Kurokawa River gorge. For the majority of the cases, horizontal deformations are in the NE directions and the vertical deformations are subsidence, indicating that Points 1 and 2 were on the hanging-wall side of the rupture. It is noted that there are some cases that result in SW horizontal deformations and uplifts.
It is important to note that the deformation results shown in Fig. 13 maintain physical dependency of the deformations at the two target locations. Therefore, differential deformation demands between the two sites can be calculated for each of the 1000 refined models and can be used to investigate the probabilistic ground deformations for linear engineering structures that pass through these two points. Such results are indeed valuable for more advanced engineering investigations. Fig. 14 shows histograms of the simulated differential ground deformations between Point 1 and Point 2 based on the 1000 refined source models. The horizontal differential deformations are calculated as the root mean square of the EW and NS deformations for each source model, whereas the vertical differential deformations are the absolute values of the differences of the UD deformations at the two locations.
The results indicate that the mean horizontal differential deformation is 0.08 m and for the majority of cases, horizontal differential deformations are less than 0.2 m, whereas the mean vertical differential deformation is 0.09 m and vertical differential deformations are typically less than 0.15 m. The results give quantitative estimates of the ground deformation demands to the Aso Bridge that are caused solely by tectonic fault movements. It is noteworthy that, although it is rare, there are extreme cases where very large differential deformations, exceeding 0.5 m, are predicted; these cases correspond to situations where the fault strike cut right underneath the Aso Bridge and western and eastern sides of the bridge move in the opposite directions.

Conclusions
A novel stochastic method for evaluating tectonic ground deformations was developed in this study. The computational framework integrates probabilistic models of earthquake source parameters, stochastic generation of earthquake slip, and deterministic estimation of surface ground deformations via Okada formulae. The method can be applied to both forecasting and hindcasting situations and are particularly suitable for applications where uncertainties of ground deformation demands at a single location and differential ground deformation demands at two locations need to be quantified. When applied to retrospective investigations, the method can be viewed as sampling-based source inversion methods to derive multiple source models that produce better predictions of the ground deformations at observed sites. An important advantage of the developed method is that spatial dependency of the ground deformation fields is fully accounted for.
The developed method was applied to the case study for the 16 April 2016 Kumamoto earthquake in Japan. During and after the earthquake, various geophysical observations of ground deformations and failures were recorded. Using a wealth of available ground deformations from strong motion networks and permanent/temporary GPS networks, existing earthquake source models developed by the GSI, which assign uniform slip values to a few rectangular finite-fault plane sources, were improved by generating numerous seismologically justifiable source models and comparing ground deformations at the observation sites. In total, 1000 refined source models were derived for the Kumamoto mainshock, and their source characteristics, i.e. fault plane geometry and spatial slip characteristics over the fault plane, were investigated in detail. Subsequently, probabilistic characterization of ground deformation demands was demonstrated for two locations near the Aso Bridge in Kumamoto Prefecture, which had collapsed due to the earthquake. The results can produce quantitative estimates of the ground deformations and thus are useful as input for advanced engineering analyses.
The developed tools can be applied to other case studies where ground deformation demands due to tectonic fault movements need to be estimated, and can also be extended. For instance, the method can be applied to quantitative risk assessments of linear engineering structures (e.g. long-span bridges and pipelines) that cut across active fault areas. The method can also be applied to subduction earthquakes by replacing the scaling relationships for shallow crustal earthquakes (Eqs. (3)-(9)) with those for mega-thrust subduction earthquakes [22]. For the latter, the method can be integrated with probabilistic earthquake source modeling (as typically done in probabilistic seismic hazard analysis) by considering a range of earthquake magnitudes/scenarios. Essentially, this will facilitate the development of an innovative fully-probabilistic procedure of characterizing tectonic ground deformations caused by earthquakes and can contribute towards the further development of performance-based earthquake engineering framework.