• from Someya, Yomada, & Okazaki (2025:3)

To simplify the use of functions and to easily calculate gradients and Hessians, we provide a high-level wrapper class, OkadaWrapper. This class abstracts the low-level subroutines (SPOINT, SRECTF, DC3D0, DC3D) and offers a unified interface for forward modeling and derivative computation.

This class supports three main methods: compute, gradient, and hessian. Each method accepts the same basic inputs: observation coordinates (coords) and fault parameters (params), with additional arguments specific to each method. In the following subsections, we describe each method in turn, along with examples.

• from Someya, Yomada, & Okazaki (2025:3-5)

The compute method performs the forward calculation; given the fault parameters, the displacements and/or strains are calculated.

Required arguments are as follows.

• coords: A Python dictionary containing observation station coordinates. Allowed keys are "x", "y", and optionally "z" (all PyTorch tensors of the same shape). The coordinate system is right-handed: x (east), y (north), z (upward). Note that z must be negative for subsurface observation points.

• params: A dictionary specifying the fault parameters. Required keys include:

  – "x_fault", "y_fault", "depth": x, y coordinates and depth of the source (depth is positive). In the case of a point source, these values represent the location of that point. In the case of a rectangular fault, the flag fault_origin specifies which point these values represent. If fault_origin is "topleft", then "x_fault", "y_fault" and "depth" represent the coordinates of the top left corner of the rectangle. If fault_origin is "center", then "x_fault", "y_fault" and "depth" represent the coordinates of the rectangle’s center.

  – "strike", "dip", "rake": fault orientation.

  – "slip": slip amount (for rectangular faults) or potency (for point sources).

Optional keys: "length" and "width" for rectangular faults.

Other optional arguments are as follows.

• compute_strain: whether to compute strain components in addition to displacements. Default True.

• is_degree: whether strike/dip/rake are in degrees. Default True.

• fault_origin: which point the fault location parameters ( "x_fault", "y_fault" and "depth" ) refer to. Either "topleft" (default) or "center" can be specified.

• nu: Poisson’s ratio of the medium. Default 0.25.

A typical usage is as follows:

from OkadaTorch import OkadaWrapper

coords = {
    "x": x,  # torch.Tensor representing x-coordinate of the station
    "y": y   # torch.Tensor representing y-coordinate of the station
}

params = {  # All values are torch.Tensor (each is a scalar)
    "x_fault": x_fault,
    "y_fault": y_fault,
    "depth": depth,
    "length": length,
    "width": width,
    "strike": strike,
    "dip": dip,
    "rake": rake,
    "slip": slip
}

okada = OkadaWrapper()

out = okada.compute(coords, params)

If compute_strain is True, out is a list of 12 tensors: 3 displacement components and 9 strain components:



If False, only 3 displacement components are returned:



Figures 1 and 2 present a demonstration of the okada.compute method applied to a rectangular fault model. The fault parameters are taken from the model 10 of Table S1 in [32]. Figure 1 illustrates the surface displacement field, while Figure 2 shows the corresponding strain distribution computed at the surface. The result confirms that the wrapper produces physically consistent deformation patterns.

• from Someya, Yomada, & Okazaki (2025:5)

The gradient method computes the derivative of the model output with respect to a single input variable. PyTorch’s function jacfwd is used internally.

• arg: The variable to differentiate with respect to. This must be one of the keys in either coords (e.g., "x", "y", "z") or params (e.g., "depth", "strike", etc.).

If arg is one of the fault parameters, the resulting first-order derivatives represent the sensitivity of the output with respect to that parameter. These sensitivities can be directly used for gradient-based first-order optimization methods, as well as for local sensitivity analysis and parameter studies.

Other arguments (coords, params, compute_strain, is_degree, fault_origin and nu) are the same as in the compute method.

A typical usage is as follows:

out = okada.gradient(coords, params, arg="x")
out = okada.gradient(coords, params, arg="depth")

If compute_strain is True, out is a list of 3 displacement components and 9 strain components differentiated by arg:



If False, out is a list of 3 displacement components differentiated by arg:



Figure 3 shows the first-order derivatives of the vertical displacement u_z with respect to nine fault parameters (same parameters as in Figures 1 and 2). As expected, the derivative with respect to slip — being linearly related to the output — exhibits the same spatial pattern as u_z itself. In contrast, derivatives with respect to the other nonlinear parameters display spatially distinct and interpretable structures. Such visualizations may also be useful for educational purposes, for example, to help students new to seismology understand how each fault parameter influences observed surface deformation.

• from Someya, Yomada, & Okazaki (2025:6-7)

The hessian method computes the second-order derivatives of the model output with respect to two variables. PyTorch’s function jacfwd is used internally.

• arg1, arg2: The variables to differentiate with respect to. These must both be from the same category (either both from coords or both from params). For instance, ("x", "y") or ("depth", "rake") are valid pairs, but ("x", "slip") is not.

If both arg1 and arg2 are fault parameters, the resulting second-order derivatives quantify the local curvature of the forward model output with respect to those parameters. These second derivatives are useful for second-order optimization, uncertainty quantification, and Laplace approximation methods in Bayesian inference, where knowledge of the Hessian around the optimum is essential.

Other arguments (coords, params, compute_strain, is_degree, fault_origin and nu) are the same as in the compute method.

A typical usage is as follows:

out = okada.hessian(coords, params, arg1="x", arg2="y")
out = okada.hessian(coords, params, arg1="depth", arg2="depth")

If compute_strain is True, out is a list of 3 displacement components and 9 strain components differentiated by arg1 and arg2:



If False, out is a list of 3 displacement components differentiated by arg1 and arg2:



Figure 4 shows diagonal entries of the Hessian of u_z, corresponding to cases where arg1 = arg2. As expected, the second derivative with respect to slip is zero, since u_z depends linearly on slip. In contrast, the second derivatives with respect to nonlinear parameters exhibit distinct spatial patterns, reflecting the complex ways in which each parameter influences the surface deformation.

• from Someya, Yomada, & Okazaki (2025:7-10)

In this subsection, we demonstrate how fault parameters can be estimated from observed displacement data using gradient-based optimization. Rather than using real geodetic observations, we generate synthetic data by adding random noise to surface displacements obtained from forward modeling.

We begin by computing displacements using the same fault parameters as in Figure 1. Gaussian random noise with an amplitude of 0.05 m is then added to each component (u_x, u_y, u_z) (Figure 5). Given these synthetic observations, we estimate the fault parameters by minimizing the misfit between observed and calculated displacements. The optimization is performed using the Adam optimizer over 2000 epochs. The core part of the code is excerpted as follows:

import torch
from OkadaTorch import OkadaWrapper

okada = OkadaWrapper()

params = {  # initialization
    "x_fault": torch.tensor(
        0.0,
        requires_grad=True
    ),
    "y_fault": torch.tensor(
        10.0,
        requires_grad=True
    ),
    "depth": torch.tensor(
        1.0,
        requires_grad=True
    ),
    "length": torch.tensor(
        150.0,
        requires_grad=True
    ),
    "width": torch.tensor(
        60.0,
        requires_grad=True
    ),
    "strike": torch.tensor(
        200.0,
        requires_grad=True
    ),
    "dip": torch.tensor(
        45.0,
        requires_grad=True
    ),
    "rake": torch.tensor(
        300.0,
        requires_grad=True
    ),
    "slip": torch.tensor(
        10.0,
        requires_grad=True
    )
}

optimizer = torch.optim.Adam(
    [
        p
        for p in params.values()
        if p.requires_grad
    ],
)

for iter in range(2000):

    optimizer.zero_grad()

    ux, uy, uz = okada.compute(
        coords,
        params,
        compute_strain=False,
        is_degree=True,
        fault_origin="topleft"
    )

    loss = 0.5 * (
        (ux - ux_obs) ** 2 +
        (uy - uy_obs) ** 2 +
        (uz - uz_obs) ** 2
    ).sum()

    loss.backward()

    optimizer.step()

Note that we do not explicitly call the gradient method here. Since the entire model is written in PyTorch, gradients can be automatically computed by defining a loss function and calling loss.backward(). This allows for seamless integration into PyTorch’s optimization pipeline.

Table 1 (middle row) and Figure 6 show the initial fault parameter values and the corresponding displacement components. Table 1 (bottom row) and Figure 7 show the final optimized parameters and the displacement components they produce. The evolution of the loss function is also shown in Figure 8.

All estimated fault parameters are reasonably close to their true values, indicating successful convergence of the optimization. However, some trade-offs between parameters are apparent. For example, both depth and slip are overestimated, suggesting that similar surface displacements can be explained by different combinations of fault geometry and slip amplitude. We note that for some parameters, particularly the fault location parameters x_fault and y_fault, the choice of initial values can have a strong impact on optimization. If the initial guess is too far from the true values, the solution may diverge. In practical applications, it may be advisable to fix such sensitive parameters and optimize only the remaining ones.

The purpose of this example is not to propose a robust inversion framework, but rather to illustrate that gradient-based parameter optimization is straightforward when using a fully differentiable implementation. Applying this technique to real data would require more advanced treatment.

• from Hinzen (2009)

• from Hinzen (2009:2855)

Since the early days of modern seismology, toppled artifacts such as tombstones and single columns have been used in the aftermath of earthquakes to deduce parameters of site-specific ground motions. The artifacts were generally treated as rigid bodies. Later, the theory of rigid block movements was also applied to precariously balanced rocks toppled by earthquakes. While the movements of a single rocking block can be described analytically, slide-rocking movements, bouncing, and multiple block systems require a numerical approach. We use multiple rigid block models with viscoelastic coupling forces in combination with full 3D ground motions (measured and synthetic) to analyze the dynamic response of building elements, relevant for archaeoseismological studies. First, the numeric modeling results are verified by comparison with analytically determined rocking motions of a single rectangular block. Stiffness and damping parameters of the coupling forces are adjusted to results from analog experiments with a rocking marble block. A model of a monolithic column and one consisting of seven drums is used to test the influence of the geometry and friction on the toppling behavior. The main question addressed in this study is whether toppled columns give a clear indication of the back azimuth toward the earthquake source. Input motion from 29 strong-motion records indicates little correlation between downfall directions and back azimuth. Clearly directed horizontal ground movements tend to topple the columns in the transverse direction. More complex ground motions result in quasi-random downfall directions. The friction coefficients have a minor influence on the downfall directions. Synthetic ground motions for two earthquakes with different source mechanism show toppling directions toward and away from the source as well as in the transverse bearing. However, it is not straightforward to deduce a reliable source location from the inversion of the toppling directions.

• from Hinzen (2009:2855-2856)

In the study of preinstrumental earthquakes, historical seismology and palaeoseismology are well-established branches of seismological sciences. Techniques to evaluate ground motions and parameters of causing earthquakes that have left their mark in written documents and in the near-surface geology were developed. Ever since man-made structures have been erected, earthquakes have also left their marks on these constructions. However, damages in archaeologically excavated buildings or continuously preserved monuments are often hard to unravel in terms of the causative effects. The new branch of seismological sciences, archaeoseismology, is defined as “the detailed study of preinstrumental earthquakes that, by affecting locations of human occupation and their environments, have left their mark in the archaeological record” (Buck and Stewart, 2000). Following this definition, the detailed study of earthquakes is the focus, and compilation, modeling, and interpretation of damage data is a means to an end.

The main questions to be answered by archaeoseismic investigations are (1) how probable is seismically induced ground motion as a cause of damage observed in man-made structures from the past, (2) when did the damaging ground motion occur, and (3) what can be deduced about the nature of the causing earthquake (Galadini et al., 2006). All three questions should be answered before results from archaeoseismic case studies can be included in seismic hazard analyses. While the first problem requires input from multiple disciplines including civil engineering, geophysics, and geotechnics, the second is a task for the geological sciences and archaeometry. The third question requires input from several seismological specialties ranging from seismotectonics to source studies, wave propagation, and site-effect modeling to soil–building interaction. The two main parameters to be extracted from these studies are the location and a measure of the strength of the causing earthquake, where both are often intrinsically tied to each other. While the macroscopic degree of damage is directly connected to the nature of ground motion at the site, many techniques from macroseismic data analysis can be adopted; determination of the direction of ground motion that caused damage or the wave propagation direction, which is even more complicated, remains a challenge in archaeoseismology.

Several case studies have been published in which surface rupturing of the fault plane directly affected man-made structures (i.e., Galadini and Galli, 1999; Galli and Galadini, 2001; Meghraoui et al., 2003; Galli et al., 2008; Marco, 2008). In these cases the location of the activated fault section is evident and source parameters such as surface rupture displacement can, under favorable conditions, be determined within small ranges of uncertainty. However, in cases where the activated segment of a fault is remote from the site with archaeological findings of building damages (archaeo-damages), source identification is a challenging task. No systematic studies or established methods presently exist to deduce the location of the rupturing fault from archaeo-damages. Several case studies purported to infer back azimuth to the activated parts of a fault from directional features in archaeo-damages (Korjenkov and Mazor, 2003). However, these interpretations usually suppose a very simple behavior for site-specific ground motions and do not consider the complexity of an extended seismic source and the uncertainties imposed by either the complexity of or randomness in the reaction of building components.

In the seismic design of modern buildings, finite element models are frequently used to study the behavior of the complete structure or crucial components of the structure (i.e., Meskouris et al., 2007). Most important in these kinds of studies are the correct modeling of the elastic or nonlinear material behavior to deduce the capacity of concrete beams, steel frames, shear walls, or wooden beams (Hinzen and Weiner, 2009). While finite element calculations predict the dynamic load under which a certain component will fail, the movement of structural parts of a building in the collapse phase cannot be modeled. The latter parameter is important to analyze directional features in archaeo-damages. Well-preserved ancient buildings were often constructed from blocks of natural stone, sometimes without any cementation. This was often the case in classic Greek buildings. For such structures it is feasible to use rigid block models and Newtonian mechanics to study the block movement and interaction as a first modeling attempt (Sinopoli, 1995; Papantonopoulos et al., 2002; Psycharis, 2007). A related problem is that of the dynamic stability of precariously balanced rocks, which have been used as low-resolution strong-motion seismoscopes (Brune, 1996). Anooshehpoor et al. (2004) used a numeric approach to model the 2D dynamic response of rocks to arbitrarily complex acceleration time histories (Purvance et al., 2008).

In this study we adapt rigid block techniques to the needs of archaeoseismology. Specifically, we regard linear viscoelastic coupling forces with finite friction in multiple block systems and true 3D ground-motion excitation. We use a numeric rigid block model of two cylindrical columns regarded as simple archaeo-seismoscopes. One monolithic column and one column consisting of seven separate column drums is used to study systematically the effects of (1) the columns’ slenderness, (2) changes in frictional coefficients, and (3) variability of measured and synthetic ground motions on the downfall directions.

• from Hinzen (2009:2856-2857)

In his famous work about the Neapolitan earthquake in 1857, Robert Mallet (1862) not only prepared the ground for evaluating earthquake strength with macroseismic methods, he also tried to infer the earthquake location from directional damage features. Without a scientific basis, still being actively sought (Ambraseys, 2006; Marco, 2008), the practice suggested by Mallet (1862) should not be applied. While Mallet used fresh traces of directional damage, in archaeoseismology such features have gone through altering processes, making it more difficult to deduce accurate directions toward the earthquake source.

As summarized by Galadini et al. (2006, and references therein), typical earthquake effects on constructions are (1) cross fissures in the vertical plane due to shear forces, (2) corner expulsion due to orthogonal motion of walls, (3) lateral and rotational horizontal and independent motion of blocks within a wall, (4) height reduction due to vertical crashing, (5) deformation of arch piers including collapse of key stones, (6) wall tilting and distortion, and (7) rotation or toppling of pillars or parts of it and drums of columns. Additional photos of typical damages are given by Marco (2008). The direction of any of these seismogenic damage patterns always results from the coaction of the orientation of the structure or structural component and the orientation of the ground movement, which is influenced by the source characteristics and site conditions.

Korjenkov and coworkers (i.e., Korjenkov and Mazor, 1999a,b, 2003; Al-Tarazi and Korjenkov, 2007; Korzhenkov et al., 2009) have given examples how directionalities in archaeoseismic damage can be quantified for a certain site. This quantification works best when carried out during an ongoing excavation. Because this is not always feasible, measurements are taken from preserved ruins (Galadini and Galli, 2001) or from the documentation of former excavations (i.e., Hinzen and Schütte, 2003). Korjenkov maps the direction of cracks and rotations of blocks with respect to the trend of walls. At sites with numerous damaged walls, a statistical approach can reveal preferred directions in the damage pattern. While such preferred directions of block shifts, rotations, toppling of wall fragments or complete walls, and toppling of columns provide strong arguments for significant ground-motion amplitudes in a certain direction, it is not straightforward to deduce the back azimuth to the causing earthquake from this direction. Near-fault strong ground motions are influenced by the source mechanism, rupture process details, distribution of asperities, fault plane geometry and extension, wave spreading conditions, and site conditions (i.e., Erdik and Durukal, 2004, and references therein). Also, secondary earthquake effects due to deformation of soft subsoils that form directional features such as foundation cracks and inclined walls (Hinzen and Schütte, 2003) should not be mistaken for the bearings toward the earthquake source.

Among the most obvious and promising directional archaeo-damages are toppled columns (i.e., DiVita, 1996; Nur and Ron, 1996; Marco, 2008) and columns with displaced column drums (Stiros, 1996; Bottari, 2005; Psycharis, 2007). As rotation-symmetrical construction elements, columns should react similarly to ground motions in any direction when freestanding and not connected to neighboring structural elements. This feature makes columns a good universal seismoscope. Even though such freestanding columns are rare in cultural heritage, the independence from a fixed trend of the structure with respect to unknown ground motions makes them a versatile tool to study basic toppling effects. Even the simplest object, a monolithic block on a plane subsurface, demonstrates a complex dynamic behavior including stress discontinuity, structural damping, contact friction, and impact (Sinopoli, 1995); all of which can influence toppling behavior and downfall directions. Many open questions remain about the particular influence of each of these factors.

Several authors have studied the dynamic behavior of cylindrical structures and classical columns. Koh and Mustafa (1990) studied the free rocking motion of rigid cylinders for various initial conditions and a stationary foundation during the motion of the cylinder. They numerically integrated the exact equations of motion of the model and mapped the boundary between toppling and not toppling. Koh and Hsiung (1991a,b) extended the 2D model to 3D rocking, rolling, and uplift of a rigid cylinder when subjected to ground motions, showing that 3D motion is significant under earthquake-like excitations. Mouzakis et al. (2002) used a 1:3 analog model of a multidrum column from the Parthenon of the Acropolis of Athens, even made from the same material as the original. Scaled earthquake ground motions in two- and three-dimensions were used to drive a shaking table with forces insufficient to topple the model. They found large deformations during the shaking, which were not necessarily reflected by the residual displacements at the end of the tests. A significant influence of imperfections of the model specimen and a very high sensitivity to even small changes in the input motion parameters were observed. In an accompanying article, Papantonopoulos et al. (2002) successfully used the distinct element method to numerically simulate the behavior of the same column. Konstantinidis and Makris (2005) also used numerical models of multidrum columns represented by a 2D discrete element model allowing rocking, sliding, and slide-rocking to show that relative sliding between drums happens even when the g-values of the ground accelerations are less than the coefficient of friction. They concluded that typical classical columns can survive the shaking from strong ground motions near the causative fault of earthquakes with moment magnitudes 6.0–7.4. Additionally, they found a more controlled seismic response of multidrum columns than monolithic columns of the same size. Psycharis (2007) made a backward analysis of a column of the temple of Olympios Zeus in Athens to investigate the seismic history of the area during the 2000 yr since the building was erected. A 3D numeric model of a single and double column structure allowed the author to constrain the maximum ground velocities at the site by comparing the model behavior with the current status of the real structure.

By applying archaeoseismologic methods, the orientation of damaged structures should be recoverable in most cases (i.e., Korjenkov and Mazor, 1999a), although usually little or nothing is known a priori about the nature of the ground movement. For elongated structural elements (i.e., a simple freestanding wall) the angle between the polarization of the ground motion and the trend of the element is decisive for estimating dynamic reaction of the structure. In contrast, rotation-symmetrical building elements (i.e., freestanding columns) do not exhibit this dependency. Therefore, to investigate correlations between ground-motion polarizations, back azimuth, and toppling directions, we will limit this study to the case of freestanding cylindrical columns.

• from Hinzen (2009:2857-2859)

Long before a strong-motion instrumentation was available, Milne (1881, 1885), Perry (1881), and others used the theory of dynamic block structures to deduce earthquake ground accelerations from toppled monumental columns and tombstones. A fundamental article on the theory of rigid block movements by Housner (1963) helped to explain observations made during the large Chilean earthquake of May 1960. Augusti and Sinopoli (1992) presented a comprehensive summary on the dynamic modeling of large block structures and Sinopoli (1995) reviewed studies of large block structure dynamics. Brune and Whitney (1992), Brune (1996), Anooshehpoor et al. (1999, 2004), and Zhang and Makris (2000) applied rigid block movement models to interpret precariously balanced rocks and a steam engine, the latter overturned during the great San Francisco earthquake of 1906.



We briefly introduce the problem of rocking rigid block dynamics in order to describe the analog and numeric experiments presented in the main part of this article; a comprehensive description of the theory can be found in Housner (1963) and Augusti and Sinopoli (1992). Figure 1 shows the cross section of a rigid body in a rocking motion. The block has height h and width b; the half-diagonal,

R = √[(h/2)² + (b/2)²],

is defined as the distance from the center of gravity, cg, to the actual center of rotation of the block, A and A′, respectively. The angle between the longer block side, which coincides with the vertical direction when the block is at rest, and the line R is α, and the inclination (rotation) of the block is θ. Housner (1963) solved the equations of motion for the rocking block for one degree of freedom, namely the rotation around the corner points A and A′ (Fig. 1d). This movement constraint specifies that (1) the friction at the corner points is large enough that no sliding between the block and the base occurs as in Figure 1c and (2) the block does not bounce during the movement through the static equilibrium (Fig. 1f). The equation of motion of the free rocking response of a 2D rigid block without sliding, slide rocking, or free flight behavior was given by Housner (1963):

I₀ θ̈ = −WR sin(α − θ) (1)

in which

I₀ = (4/3)mR²

is the mass moment of inertia of a homogeneous rectangular block about corner A and

W = mg

is the weight of the homogenous block, m and g are the block mass and acceleration of gravity, respectively.

For slender blocks (sin α ≈ α), equation (1) can be written as (Housner, 1963; Aslam et al., 1980):

θ̈ − p²θ = −p²α (2)

in which

p = √(WR/I₀) = √(3g/4R).

In order to start a forced rocking motion, a horizontal acceleration a is necessary that fulfills the condition (Housner, 1963; Augusti and Sinopoli, 1992):

a/g ≥ b/h (3)

Under more realistic conditions, energy dissipates during each impact of the block on the base. For inelastic impact (no bouncing), the reduction of kinetic energy, r, depends on the square of the ratio of angular velocities before, θ̇_i, and after, θ̇_i+1, an impact, respectively (Housner, 1963; Aslam et al., 1980; Augusti and Sinopoli, 1992):

r = [(1/2)I₀θ̇_i+1²] / [(1/2)I₀θ̇_i²] = (θ̇_i+1/θ̇_i)²  (4)

and for slender blocks the coefficient of restitution is:

η = √r = (2 − b²/h²) / [2(1 + b²/h²)] ≈ 1 − (2/3)α²  (5)

Housner (1963) showed that in this case the amplitude θ_n of the nth cycle can be written as:

θ_n = 1 − √[1 − rⁿ (1 − (1 − θ₀/α)²)]  (6)

and the half period of the rocking motion is:

T/2 = 2√(I₀/WR) tanh⁻¹ √[rⁿ (1 − (1 − θ₀/α)²)]  (7)

Amplitudes and the half period of the rocking motion decay rapidly with the number of cycles. Increasingly high-frequency movements follow a few slow rocking motions with large amplitudes. Under less ideal conditions (bouncing during impacts and sliding due to reduced friction), amplitude decay proceeds even more rapidly.

The sliding component of the movement is governed by the static coefficient of friction, μ_s, and the slenderness of the block. For the free motion of the block Augusti and Sinopoli (1992) showed that the inequality:

μ_s ≥ [3(b/h)] / [4 + b²/h²]  (8)

allows separating the conditions under which pure rocking and slide rocking exists. If the condition of equation (8) is not fulfilled, slide rocking, controlled by the kinetic friction coefficient μ_k, will start. For small angles for each b/h, a value μ_s exists above which only rocking will occur until the block returns to a state of static equilibrium (Augusti and Sinopoli, 1992).

• from Hinzen (2009:2871-2873)

Depending on the perspective (i.e., archeological versus seismological), different expectations from archeoseismological studies usually result. Archaeology seeks explanations for damage horizons in excavations, interprets the impact of the potential earthquake on cultural development, and recovers as much history as possible from a certain place or region. Seismology focuses on the causing earthquake itself and the mechanisms that produced the excavated damage pattern. In order to make archaeoseismic data useful for seismic hazard analysis, the questions of when, where, and how strong the causative earthquake was, are the most important ones to be answered. While the "when" is usually beyond the predictive capacity of seismological techniques, the "where" and "how strong" are challenging but appropriate tasks for seismologists.

The distance and strength of an archaeo-earthquake are intrinsically coupled and cannot be resolved in general. Even if this obstacle could be overcome, the question of a direction toward the epicenter remains open. Past attempts, including several case studies, concentrated on surmising this direction from directional damage patterns. Unfortunately, systematic approaches to develop archaeoseismic methods that make use of the modeling capabilities of engineering seismology are rare. Two arguments favor the use of rigid blocks as a first step in such models. (1) Clear directional archaeo-damages are often found in displaced large block structures, especially those that might have been constructed without cementing. Even if mortar was used to cement blocks, the tensile strength of the cementation is usually low and after an initial breaking, the mortar mainly influences the frictional parameters of interblock movements. (2) As only the final resting or downfall positions are preserved in the archaeological record, Newtonian mechanics are necessary to track the path of disintegrated wall or column elements and determine their impact positions. A clear shortcoming of rigid block models is the inability to forecast direction and degree of fracturing of intact blocks or structural elements. Additionally, it is usually much harder to link these damages to seismogenic causes than laterally or rotationally displaced blocks and toppled features.

We use a numeric model of two columns, monolithic and structured in seven drums, composed of rigid blocks and coupled by viscoelastic forces. Stiffness and damping parameters of the coupling forces are determined through the comparison with the measured rocking movement of a small analog marble block. Overturning of simple block structures by simple ground motions like single-sine pulses can still be handled analytically, provided that friction between the block and base is large enough to avoid sliding. Our numeric model successfully reproduces results of such analytical solutions. Effects of more complex ground motions (true earthquake movements), especially on multiple block systems, including limited friction, require the application of numerical models (i.e., Konstantinidis and Makris, 2005; Psycharis, 2007).

Test calculations with a variation of the slenderness of the columns in a wide range and subjected to measured strong ground motions of different character show clear differences for strongly and weakly linear polarized ground motions in the horizontal directions. If a clear directed pulse in a certain horizontal direction dominates ground motion, this pulse tends to topple the columns independent of the slenderness roughly in the same direction. In a case where the ground motion does not show a clear directional impulse, downfall directions are scattered and large changes occur in downfall directions due to small differences in the column slenderness, which supports the conclusions from Mouzakis et al. (2002). As larger horizontal strong ground motions are usually bound to the S phase, toppling directions tend to occur in these cases at angles roughly orthogonal to the back azimuth.

A variation of the frictional conditions with otherwise constant parameters shows that at least for more impulsive or clearly directional ground movements the friction has little influence on downfall directions. Directions vary more for the structured column than for the monolithic one. This makes sense, as the total number of degrees of freedom in case of the seven-drum column is larger. In the test case with clear pronounced ground-motion directions, the downfall directions for both column types show an angle of roughly 90° with respect to the back azimuth. Even for these cases, an estimate of the correct direction towards the earthquake source based only on one observation would not have been unequivocal. For more complex ground motions, lacking a clear directional pulse, small changes in the friction coefficient determine the downfall direction of the structured column. The downfall of the monolithic column is then spread over a larger angular range and even includes the correct back azimuth. Without knowledge of the nature of the causative ground motion, almost always the case in archaeoseismic case studies, there would be no way to determine whether the toppled column fell in the direction of the back azimuth or at an ±90° angle with respect to it. These observations are in agreement with results from earlier studies. Mouzakis et al. (2002) and Papantonopoulos et al. (2002) found in analog and numeric experiments with marble columns a large sensitivity of the response to small changes in the input motion parameters. Shortcomings of the two-column model used in this study such as simplification of the contact forces and neglecting elastic deformation influence the dynamic reaction. An approximate prediction of overturning behavior, however, is possible.

Korjenkov and Mazor (2003) correlate directional damage patterns to the source distance. They assume that directly above the source vertical movements dominate, resulting in severe damage patterns of random directions. In contrast, for sites that are at some distance from the epicenter, they assert that lateral movements become significant, producing tilting and collapse towards the epicenter. For the simple column model and synthetic ground motions from two earthquakes with a strike-slip and normal faulting mechanism, no clear distance dependent trend of the downfall directions of column drum displacements was found. For earthquakes strong enough to excite significant archaeo-damages, radiation and directivity effects are crucial for determining the character of site-specific ground motion and associated directions of damage patterns. Local subsurface conditions, although not specifically addressed in this study, have demonstrable influence on strong ground motions. For soft sediment upper layers, the soil-structure interaction becomes another crucial factor, complicating the damage process and potential deformations as well as throw directions.

The previously proposed measure of the predominant acceleration pulse duration in form of the PGV/PGA ratio in the study of precariously balanced rocks (Purvance et al., 2008) also works as an indicator of the overturning potential of simple columns where PGD should also be taken into account. This is a promising approach to determine site-specific ground-motion parameters from archaeo-damages. The Arias intensity can be used as an additional parameter quantifying the overturning potential.

The current study is a basic approach to systematic parameter studies. With the proposed techniques applied to actual field cases, structured and monolithic column differences must be addressed as part of the analysis of downfall position as shown for example by Psycharis (2007). A structured column toppled domino-style indicates the original impact locations of the drums. Model calculations have shown that drums can make significant postimpact movements on a stiff surface, frequently resulting in disintegration of the structure. Monolithic columns may move after impact, for example, bouncing and/or rolling, but this behavior is strongly dependent on the nature of the impact surface. If no clear impact traces are visible, the archaeologically documented position could be misleading as the column may have moved after impact (Mallet, 1862).

In the numerical approach of this study, we assumed the idealized situation of a flat column base resting on a flat subsurface without any clamping devices. Comparisons with measured movements of a rocking marble block showed that imperfection of the contact surface alter the results, at least for small motion amplitudes as also shown by Mouzakis et al. (2002). In actual field applications, the situation is usually more complex. Stiros (1996) presents schematic examples of different column settings and Sinopoli and Sepe (1993) examined the coupled motion of a three-block structure with the geometric features of the colonnade of the Selinute Temple in Sicily. For the presence of clamping forces at the column base, a cracking and shear failure of the column is more probable than if no clamping forces are present. In these cases, rigid block models might not be sufficient. Konstantinidis and Makris (2005) showed that wooden clamping poles have only a minor effect and that replacing them by stiff metallic shear links during retrofitting might even be counterproductive. Colonnades of a temple, a frequent structural element in Greek and Roman architecture, exhibit different reactions to motions parallel or oblique to the trend (Psycharis, 2007).

The ground motions in this study were intentionally limited to translational. The concomitance of rotational ground-motion components is a complicating factor and needs further research particularly with regard to initiation of structural damage and the displacement of blocks. The importance of such results extends beyond the interests of archaeoseismology. Widening the tests to incorporate realistic synthetic ground motions of extended sources might help to better constrain source parameters from archaeo-damage patterns.

In conclusion, calculated downfall directions of simply structured columns with clear boundary conditions show that deducing the direction toward an earthquake source from the positions of toppled columns is not straightforward. Ground-motion variability is the most problematic element. Only when a clear horizontal pulse of sufficient amplitude dominates the ground motion can coherent toppling directions be expected for columns of different geometric structure. Even in these cases, however, the relation between the throwing ground-motion direction, indicated by the downfall, and the back azimuth is not necessarily directly evident because fault plane orientation and structure, directivity, rupture velocity, as well as local site effects, all have their share in contributing to the ground-motion characteristics.

• from Hinzen (2009:2873)

Strong-motion seismograms used in this study were retrieved from the Web page of the Consortium of Organizations for Strong Motion Observation Systems, www.cosmos-eq.org (last accessed November 2007) and the European strong-motion database CD-ROM, European Council, Environment and Climate Research Program ENV4-CT97-0397, 2000.

• from Somerville et al. (1997:202-204)

The ground motion parameters that are modified for directivity effects include the average horizontal response spectral acceleration; the average duration of the two horizontal acceleration time histories; and the ratio of strike-normal to strike-parallel spectral acceleration. Strike-normal refers to the horizontal component of motion normal to the strike of the fault. Strike-parallel refers to the horizontal component of motion parallel to the strike of the fault. Following the method of Husid (1969), duration is defined as the time between 5% and 75% of the cumulative squared acceleration, following the convention of Abrahamson and Silva (1997b). The three ground motion parameters that we analyze, and the model parameters that control them, are listed in Table 1.



The modifications for duration and response spectral accelerations are made with respect to the empirical attenuation relations of Abrahamson and Silva (1997a,b). To account for these effects, ground motions are calculated by first using the Abrahamson and Silva (1997a,b) relations, or equivalent relations, and then applying the modification factors prescribed by the relations that we have developed. Although the models were developed based on the Abrahamson and Silva (1997a,b) models, the fact that we removed the bias between our data set and their models should make our modifications applicable to other attenuation relations. The modifications to account for the difference between strike-normal and strike-parallel response spectral acceleration are completely independent of any particular attenuation model and can be applied to the average horizontal response spectral acceleration calculated by any empirical model.

• from Somerville et al. (1997:204-206)

All of the model parameters considered in the development of the directivity functions, and their ranges of applicability, are listed in Table 2, and the functional forms used to model the ground motion parameters are shown in Table 1. Based on the description of rupture directivity effects given above, we expect amplitude variations due to rupture directivity to depend on two geometrical parameters. First, the smaller the angle between the direction of rupture propagation and the direction of waves travelling from the fault to the site, the larger the amplitude. Second, the larger the fraction of the fault rupture surface that lies between the hypocenter and the site, the larger the amplitude. We expect the duration of strong motion to be controlled by the same two parameters, with an inverse relationship between duration and amplitude. The azimuth and zenith angles and length and width ratios are illustrated for strike-slip and dip-slip faulting in Figure 5.



For strike-slip faults, the angle θ and length ratio X are measured from the epicenter to the site in the horizontal plane. For dip-slip faults, the angle φ and width ratio Y are measured from the hypocenter to the site in the vertical plane oriented normal to the fault. For sites located between the ends of the fault, the angles θ and φ, the length and width ratios X and Y, and closest distance r_rup are related as follows: r_rup/X = tan θ; r_rup/Y = tan φ.

Under far-field conditions, the quantification of the angles θ and φ is straightforward since they are assumed to be constant for all parts of the fault. However, in the near-fault environment, the fault dimensions are generally larger than the distance between the fault and the site. This causes the angle between the fault plane and the direction to the site, and consequently the radiation pattern amplitude, to continually change. We have used simple functional forms of the angles to represent the directivity effect and acknowledge the ambiguity in their relationship to ground motion parameters.

We chose to model the effects of rupture directivity on ground motion amplitudes and duration using the function X cos θ for strike-slip faults and Y cos φ for dip-slip faults. The cosine function was chosen to give the ground motion parameters a smooth decay with increasing angle. This functional form is compatible with the far-field approximation of rupture directivity effects (Aki and Richards, 1980). Preliminary analyses of the data showed that separate terms for the distance ratio and angle parameters were not required to adequately model the data.

Between the ends of a dip-slip fault, the variation of ground motion parameters with φ is indistinguishable from its variation with rupture distance r_rup. In contrast, the variation of ground motion parameters with θ is independent of the rupture distance r_rup. Since rupture distance is a primary ground motion parameter already included in attenuation relations, we expect to find more spatial variability of strike-slip motions with cos θ than of dip-slip motions with cos φ.

In a strike-slip earthquake, there are strong directivity effects off the ends of the fault. However, directivity effects are not important off the ends of a dip-slip fault because the directivity effect in this case is operating updip, not along strike. An additional parameter would be needed to describe the position of the site along strike, in order to modulate the strong directivity effect that would be predicted by the small value of the angle φ associated with stations off the end of a dip-slip fault. We have chosen not to incorporate this additional parameter, and so our model for the absolute amplitude and duration effects of rupture directivity is confined to the region between the ends of the fault. This region, which is shown in Figure 5, is the same as that used by Abrahamson and Somerville (1996) in their analysis of footwall and hanging-wall effects in dip-slip faulting.

As illustrated in Figure 2, the ratio of strike-normal to strike-parallel motions is primarily a function of the radiation pattern of the earthquake. For example, we expect (and observe) large ratios of strike-normal to strike-parallel motions in both the forward directivity direction, where amplitudes are large, and in the backward directivity direction, where amplitudes are small. We therefore expect the strike-normal to strike-parallel ratio to be mainly controlled by the angles, and less controlled by the length and width ratios.



To model the strike-normal to strike-parallel ratio, we chose to use a cos² ξ dependence of the strike-normal to strike-parallel ratio, where ξ is θ for strike-slip and φ for dip-slip faulting in the range of 0° to 45°. This form is consistent with the far-field radiation pattern for strike-normal and strike-parallel motions. In the near-fault environment, the rapid variation of angle with location on the fault can produce modulation of the strike-normal to strike-parallel ratio that is more complex than cos² ξ. Although this more rapid variation was observed in the data, we chose not to model its more complex form. We also fixed the ratio at unity for angles ξ larger than 45° in light of this complexity.