{'fraction_non_alphanumeric': 0.02990622418735727, 'fraction_numerical': 0.006875273310334775, 'mean_word_length': 4.637222678718159, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dynamical system analysis and forecasting of deformation produced by an earthquake fault

We present a method of constructing low-dimensional nonlinear models describing the main dynamical features of a discrete 2D cellular fault zone, with many degrees of freedom, embedded in a 3D elastic solid. A given fault system is characterized by a set of parameters that describe the dynamics, rheology, property disorder, and fault geometry. Depending on the location in the system parameter space we show that the coarse dynamics of the fault can be confined to an attractor whose dimension is significantly smaller than the space in which the dynamics takes place. Our strategy of system reduction is to search for a few coherent structures that dominate the dynamics and to capture the interaction between these coherent structures. The identification of the basic interacting structures is obtained by applying the Proper Orthogonal Decomposition (POD) to the surface deformations fields that accompany strike-slip faulting accumulated over equal time intervals. We use a feed-forward artificial neural network (ANN) architecture for the identification of the system dynamics projected onto the subspace (model space) spanned by the most energetic coherent structures. The ANN is trained using a standard back-propagation algorithm to predict (map) the values of the observed model state at a future time given the observed model state at the present time. This ANN provides an approximate, large scale, dynamical model for the fault.

deterministic or stochastic relationships involving only a few irreducible degrees of freedom, which may be used for short-term prediction. We show that for earthquake forecasting this can be done using the spatiotemporal strain patterns embedded in the observable surface displacements. Such an approach is based on the observation that the large scale dynamics of the system often evolves on a manifold (or invariant measure for the case of strange attractors) with a dimension that is significantly smaller than that of the system's phase space. In such cases, a few macroscopic observables can approximate very well the present state of the system and predictive models based on the dynamics of a reduced number of macroscopic observables can then be constructed.
In order to identify meaningful low-dimensional structures, we apply the Proper Orthogonal Decomposition (POD) -also known as the Principal Component Analysis (PCA) or Karhunen-Loève expansion -to an ensemble of surface deformation data generated by the system's dynamics. The POD provides the most efficient way of capturing the dominant components of a dynamical process with only finitely many, and often surprisingly few, "modes" [17]. Using synthetic calculations for a strike-slip fault system, we show that a reduced number of deformation modes can explain on the average the large scale dynamics of elastic surface deformations. The dynamics of these modes live on a low-dimensional "reduced" attractor in the neighborhood of which the system spends most of its time. The state of the system in this reduced spacemodel state space -is represented by a set of modal coefficients that measures the projection of the ground surface deformation onto each dominant mode. Our goal is to extract the nonlinear modal dynamics in this model space by constructing a map of the observed model state at a future time given the observed model state at present time. We will present preliminary results in which an artificial neural network has been used with promising success to learn the dynamics of the reduced model from these modal time series. The reconstructed map has been iterated to obtain predictions for the long term model dynamics starting from its current model state. A rough test of the reliability of the model forecast to approximate the future of the fault system is also discussed.
Our method may be compared to the linear pattern dynamics introduced by Rundle and his coworkers [25].
Their technique is based upon a Karhunen-Loève expansion of the spatio-temporal seismicity data and is used to estimate a linear stochastic model for the evolution of a probability density function for seismic activity.
In contrast to that method, which provides a local linear approximation in a probability space, we propose a global nonlinear approximation that describes the effective large-scale dynamics in a low-dimensional phase space.

Description of the Earthquake System
We distinguish between the assumed physical system and the model which is an approximate representation of the system. The assumed system is itself a simplified representation of fault systems in the earth, but will be considered a valid description of reality if its behavior is typical for the dynamics of an isolated fault. Our goal is to build data driven models, with useful predictive skills, which approximate the large scale dynamics of the system. If this is possible, we can conclude that the methods we propose can then be extrapolated to real isolated faults and, perhaps, even to complex fault systems.
Our assumed system corresponds to a discrete strike-slip fault of length L = 70 km and width W = 17.5 km embedded in a 3D elastic continuum [3,5]. The fault consists of a uniform grid of dynamical cells where slip is governed by static/kinetic friction processes, surrounded by regions with imposed constant slip rate of V pl = 35 mm/year, representing the tectonic loading. We divide the 70 km × 17.5 km computational grid into 128 × 32 square cells, with length ∆x and depth ∆z having equal dimensions of approximately 550 m, that corresponds to the dimension of small, effectively disconnected, slip patches in mature fault zones [5]. The brittle deformation at any fault position and time is governed by quasi-static, 3D elastic dislocation theory and spatially varying "macroscopic" constitutive parameters that describe the static and kinetic friction processes. The stress at any fault position (cell) increases with time, t (measured in years), due to the gradual tectonic loading and the time-dependent brittle deformation at other fault locations. It can be written as a boundary integral over the slip deficit φ kl = V pl t − u kl , where F ij,kl is the elastic stress transfer function based on the solution of Chinnery [10], u kl is the right lateral slip at fault cell (kl) ( measured in mm), and the summation is over the brittle area, B =70 km × 17.5 km, of the fault plane. The cell indexes i, k and j, l along the horizontal and vertical directions, respectively, define the spatial cell location. We assume that the static brittle strength is uniform along the fault and is set to a value of τ s = 100 bars. If the stress τ ij reaches the static strength τ s , a brittle failure occurs at this location and the cell slips an amount ∆u ij , till the local stress drops to a prescribed arrest stress level, τ a,ij .
The stress transferred from the failed cell can lead to subsequent brittle failures (i.e., rupture propagation) if the stress anywhere increased to the failure threshold. These failures may, in turn, induce or reinduce more brittle slip events. After an initial slip event the strength drops to a dynamic level, τ d,ij < τ s,ij , and reinitiation of brittle slip on an already failed cell occurs when τ ij ≥ τ d,ij there. The brittle slip associated with a stress drop can be described by where F ij,ij is the self stiffness of a cell and the failure stress τ f,ij equals the static failure stress, τ s , for a cell which has not slipped or is adjusted to its dynamic value, τ d,ij , every time a cell has failed. Figure 1  This marks the end of an earthquake event whose strength is measured by its potency P defined [7] as the integral of slip over the rupture area, where ∆A = ∆x × ∆z = 0.29 Km 2 and ∆u ij is the total slip at cell (ij). During a seismic event the tectonic slip, V pl t, remains constant. At the end of each event the failure strength everywhere on the fault recovers back to τ s . In order to initiate the next seismic event we advance in Eq. (1) the tectonic slip (or equivalently, advance the time t) till the fault cell closest to failure reaches its static failure strength. This triggers a new seismic event and starts again the failure iteration loops.
An important component of the fault rheology is the distribution of arrest and dynamic stresses along the fault. In the simulations used in this paper we chose random stresses, where < τ a > is the fault-averaged arrest stress, ξ is a random number uniformly distributed in the range [−0.5, 0.5], and A is the noise amplitude. The static strength, dynamic strength, and arrest stress are related to each other, at each point on the fault, as where D is a dynamic overshoot coefficient.

Qualitative dynamics of slips, stresses, and seismicity
We start by analyzing the evolution of seismicity and the dynamics of the first and second order moments of the stress and slip fields. Here we only investigate the dynamics of two fault realizations that differ largely in the value of the dynamic overshoot coefficient. System A with D = 1.5 has a large dynamic weakening, while system B with D = ∞ has no dynamic weakening. Both systems have random arrest stress distributions, as in (4), with < τ a >= 80 bars and A = 20 bars amplitude. A detailed investigation of the dynamics, covering the system space between these two limits, will be presented elsewhere [2]. In general, each system observable contains a different amount of information for understanding the underlying dynamics. As discussed below, for prediction purposes we are only interested in the dynamics at the large length and time scales. As the small scale dynamics overwhelmingly dominates the statistics, it is not obvious a priori which observable provides information about the large scale dynamics. In Fig. 2 we present the time evolution of the potency released over a 30 years time interval. The fault motion in system B (D = ∞) seems to be dominated by erratic, random appearing evolution. There are no detectable time patterns. When significant dynamic weakening is present, as in system A (D = 1.5), a clear signature of a quasi-regular evolution emerges. This can be identified in the presence of very large events (these are large scale events in which a large fraction of the fault slips), followed by intervals of seismic quiescence. There is a quasi-regular build up of correlations [6] which ultimately generates very large events, much larger in size than in system B. Although the evolution of seismicity is dominated by small events, there are significant evolution patterns reflecting the dynamics of the large scale events. These can be better seen in Fig. 3 which shows the evolution of the average slip deficit, and slip deficit variance, Due to the inhomogeneity of the slip deficit, the slip fluctuations are measured from a local time average, < φ ij >, defined as where the length of the averaging time used is T = 100 years. For the D = 1. the dynamics for cases with and without dynamic weakening are given in [6].
The emergence of a clear separation of length scales is also revealed in the frequency-size statistics shown in Fig.4. The D = ∞ system has a power-law distribution with an exponential cutoff, characteristic of systems close to a critical point. This confirms the analysis of Fisher et al. [11], who showed that the dynamic weakening τ s − τ d,ij is a tuning parameter and found that the system has an underlying critical point of second order phase transition at zero dynamic weakening (D = ∞) and full conservation of stress transfer during failure events. The frequency-size statistics for the D = 1.5 system shows the existence of two separate earthquake populations. At small scales we observe a power-law distribution of event sizes. The other events cover a broad range of scales and their frequency is significantly higher than an extrapolation based on the power-law statistics of the small events. This behavior is reminiscent of the characteristic earthquake model [30,7], which has a distinct probability density peak at the large size end of the frequencysize distribution. Moreover, a careful analysis [2] has shown that most stress dissipation is produced at the large size end of the potency spectrum and is due to the large scale motions. Averaged over small time intervals, the small events that dominate the frequency-size histogram do not produce any stress dissipation on average and only transfer stresses internally. This suggests a possible description of the dynamics in terms of large scale dissipative motions, that balance the stress increment due to an effective loading rate, which absorbs most of the small scale events into a renormalized plate velocity. This coarse-grained dynamics is inherently of lower dimension than the original dynamics. The question is whether the dimension collapse is significant to render data driven model reconstruction algorithms practical. scales. In such cases, we can in principle isolate the large scale dynamics which, evolving slowly, will drive the rapidly relaxing small scale dynamics. But, are real faults characterized by a large dynamic weakening?

Effective dimension of the large scale motions
Full elastodynamic simulations of rupture propagation in a 3D elastic solid [20] indicate that D ≃ 1.25, corresponding to large dynamic weakening. Moreover, we have shown elsewhere [2] that this separation of scales holds for a broad range of D values and is, therefore, generic. Consequently, restricting our study to system A, we shall address for the rest of this paper the following issues ( In order to model the large scale motions we have to find first if their dynamics is deterministic in character.
In principle, the coarse grained dynamics of any deterministic system is stochastic. When a separation of scales exists, we can assume that only the statistical properties of the small scale motions influence the larger scales, through coefficients of turbulent viscosity, and that at any moment these statistical properties are determined by the larger scale motions [19]. A model describing only the large scale dynamics is assumed to be effectively deterministic. However, we have to be careful, because, as pointed out by Lorenz [19], there remains uncertainties in these statistics, and their progression to the very large scales will ultimately produce large scale errors and limit the forecasting skill of our model. We will therefore assume, as a working hypothesis, that a model describing only the larger scales motions is effectively deterministic. In order to test our assumption and find an upper bound on the dimension of the model, we start with a phase space analysis of the large scale motions as is reflected in the scalar dynamics of the slip deficit variance.
The key to unraveling the dynamics from a scalar time series is to reconstruct a vector space, formally equivalent to the original attractor of the dynamics. For a deterministic system, Takens [18] embedding theorem guarantees that such a reconstruction is possible from a time-delay reconstruction of only one scalar observable. In our case, the embedding is performed using Cartesian coordinates made out of the observed slip deficit variance and its time delayed copies, where σ(n) = σ(n∆t) is the slip deficit variance measured at equal sampling times ∆t = 0.1 years, d is a time lag (an integer multiple of the common lag ∆t), and m is an embedding dimension. The optimal time lag, d = 3, has been chosen at the first minimum of the time-delayed mutual information as suggested in [12] and motivated in detail in [18]. The problem is to determine the optimal embedding dimension m, for which the reconstruction by delay vectors provides an acceptable unfolding of the attractor so that the orbits composing the attractor are no longer crossing each other in the reconstructed phase space [1]. The practical solution is to look for false neighbors in the embedding phase space at a given value of m [18]. To understand this concept, consider the situation that an m dimensional delay reconstruction is an embedding, but an (m − 1) dimensional delay reconstruction is not. If the embedding dimension is too small to unfold the attractor, a small R m−1 neighborhood will contain points that belong to different parts of the original attractor. Therefore, at a later time, the images of these points under the system's dynamics will split onto different groups, depending on which part of the attractor the points are originally coming from. This lack of a unique location of all the images in (m − 1) dimensions is reflected in finding false neighbors, meaning that determinism is violated. When increasing m, starting with small values, one can detect the minimal embedding dimension, m e , by finding no more false neighbors. Then, we can invoke the result of Sauer et al. [27] who showed that the attractor formed by X n in the embedding space is equivalent to the attractor in the unknown space in which the system is living if m is larger than twice the box counting dimension, d 0 , of the attractor. Often, an m = d 0 embedding dimension is enough for unfolding the attractor. Therefore, the minimal embedding dimension sets the following bounds, (m e − 1)/2 ≤ d 0 ≤ m e , for the dimension of the attractor in the true phase space. Figure 5 shows that the percentage of false neighbors for the D = 1.5 system falls under 1% when the phase space embedding of the slip variance data reaches dimension seven.
This behavior suggests a topological dimension in the true phase space between three and six. We note that the percentage of false neighbors is not completely reduced to zero. The explanation for this behavior is discussed in Section 6. To appreciate this result, we show for comparison in Fig. 5 the behavior of the system in the limit of no dynamic weakening. Even for an embedding dimension as high as ten, the percentage of false neighbors does not drop below 4%, and is an order of magnitude larger than the percentage for the D = 1.5 system at the same embedding dimension. This confirms our previous analysis and proves that in the limit of low dynamic weakening the system has high dimensional dynamics.
The foregoing analysis shows convincingly the deterministic structure of the large scale motions: for stochastic or very high dimensional dynamics the number of false neighbors will never effectively drop to zero. The low-dimension result is consistent with the embedding theorem [27] which asserts that it is not the dimension of the underlying true phase space that is important for the minimal dimension of the embedding space, but only the fractal dimension of the support of the invariant measure generated by the dynamics in the true phase space [18]. These results justify our earlier remarks that due to some collective behavior only a few dominant degrees of freedom remain and our conjecture that their dynamics is effectively deterministic.
Before addressing the problem of identifying these collective degrees of freedom, we would like to know if their dynamics is chaotic and, if so, what are the implications for their predictability. To answer this question we compute the maximal Lyapunov exponent. This exponent measures the exponential divergences of nearby trajectories and is an average of these local divergences over the whole data. A positive maximal Lyapunov exponent is a signature of chaos. Its computation is based on the algorithm described by Kantz and Schreiber [18] and uses the implementation provided in the publicly available software package TISEAN [16](all nonlinear time series algorithms described in this paper use the same excellent implementation). For a point X n of the time series in the embedding space, the algorithm determines all neighbors X n ′ within a neighborhood U n of radius ǫ. Then, for each neighbor, we consider the distance δ(0) = |X n − X n ′ | and we read its evolution δ(∆n) = |X n+∆n − X n ′ +∆n | from the time series. If for all neighbors δ(∆n) ≃ δ(0)e λ∆n , then λ is the local divergence rate. The local divergence rate varies along the attractor and the true maximal Lyapunov exponent is an average over many reference points n. Therefore, in practice we compute the average over the distances of all neighbors to the reference part of the trajectory as a function of the relative time ∆n. The logarithm of the average distance at time ∆n measures the effective expansion rate over the time span ∆n: where the outer bracket denotes the averaging over the reference point X n , |U n | is the number of neighbors in U n , and the argument of the logarithm is the average local expansion rate. If S(∆n, ǫ, m) vs ∆n exhibits a robust linear increase with identical slope for a range of ǫ values and for all m larger than some minimum embedding value, then its slope is an estimate of the maximal Lyapunov exponent λ max per time step [18]. sensitive dependence on the initial conditions. Does this dependence sets a fundamental limitation to longterm forecasting? Two initially close trajectories will diverge exponentially in the phase space with a rate given by the largest Lyapunov exponent λ max . Therefore, for any finite uncertainty δ in the initial conditions, we can forecast the future state of the system only up to a maximum time, where ∆ is the accepted tolerance level. For the earthquake system, an approximate estimate for the predictability time is T p = 1/λ max ≃ 1.25 years. This result is disappointing and, due to expected model reconstruction errors, which might overshadow the uncertainty in the initial conditions, raises questions on the value of any forecasting endeavor. But this estimate is naive and holds only for infinitesimal perturbations and in nonintermittent systems. Generally, the predictability time is scale dependent [8] and can be much longer than the rough estimation T p ≃ 1/λ max . As we will shortly see, at the large length scale of interest for forecasting, a better estimation of the predictability time is T p ≃ 10 years.

Proper Orthogonal Decomposition
The dynamics of the first and second order moments of the slip and stress fault fields provides an adequate description of the large scale motions. Unfortunately, their measurement poses very difficult problems: the stress field is not directly observable while the slip field requires the solution of an ill-posed inverse problem.
Therefore, we propose an alternative analysis of the large scale motions based on the spatio-temporal strain patterns embedded in the surface deformation fields that can be accurately measured by InSAR and GPS observations [9]. Assuming uniform slip discontinuities over each rectangular cell, the components of the surface displacement vector U α are found to be [22], where Γ α (x; x ij ) is the surface displacement in the direction α = x, y, z at the observation point x = (x, y, 0) due to uniform unit slip on the fault cell (ij) located at x ij = (x ij , 0, z ij ). The space-time signal is obtained by simultaneous measurements of the surface displacements on a 64 × 32 uniform rectangular grid which covers a 100 km × 50 km surface area centered around the fault. For each surface deformation we compute an ensemble of N snapshots, U α (x, n) = U α (x, n∆t), n = 1, . . . , N, and α = x, y, z , every ∆t = 0.1 years. The analysis of each deformation direction proceeds identically, therefore we will henceforth drop the deformation index α. Since we are only interested in decomposing the dynamics of fluctuations, it is convenient to separate the flow U (x, n) into a time-independent mean and a fluctuating part, i.e.
The identification of the active degrees of freedom in the surface deformation fields uses the Principal Orthogonal Decomposition (POD) [17]. The POD seeks an optimal representation (in the least square sense) for the members of the ensemble {u(x, n)} N n=1 . It searches for generalized directions φ(x) in the configuration phase space, such that most of the ensemble fluctuations will be directed along φ(x). Mathematically, this is achieved by maximizing the average projection of the where (·, ·) is the L 2 inner-product in the configuration space and the L 2 norm constraint, φ = 1, is required for the maximum to be defined. The solution to this variational problem is given by the eigenfunctions {φ i } M 1 of the following integral equation [17], to the smallest scales, described by the higher order modes. This observation is hardly surprising, since for a translationally invariant system the POD modes are simply the Fourier modes [17].
The basis {φ i } M i=1 is optimal in the sense that any reduced representation of the form describes typical members of the ensemble better than any linear representation of the same dimension m in any other basis: the leading m POD modes contain the greatest possible "energy" on average [17]. Therefore, due to the coherence and optimality of the POD basis, Eq. (16) provides the most efficient Euclidean (linear) embedding for the large scale motions of the system. Moreover, the modal coefficients are uncorrelated on average, i.e.
which reflects the fact that orthogonality in the embedding space, (φ i , φ j ) = δ ij , is related to the statistical properties of the time series.
The choice of the embedding dimension m is based on the computation of the cumulative normalized eigenvalue spectrum, Λ m , defined as: The cumulative spectrum can help us define an effective POD embedding dimension by finding the minimum number of modes needed to capture some specific fraction f < 1 of the total variance of the data: Figure 8 shows the cumulative normalized spectrum for each of the surface deformation modes. For system A (D = 1.5), in order to explain f = 95% of the variance in the u x , u y , u z data sets we need only 2, 3, and 5 modes respectively. This is consistent with our earlier estimates of the embedding dimension. In contrast, for system B (D = ∞), we need 9, 27, and 42 modes respectively.
It is important to realize that the POD modes carry no dynamic information. Reshuffling the ensemble of configurations -and, therefore, destroying the dynamic content hidden in their time ordering -produces the same POD modes. Therefore, in order to identify the dynamics of each POD mode we calculate the projections of the surface deformation fields onto the dominant deformation modes, i.e., This way, we generate a small number of modal time series A i (n), i = 1, . . . , m, that encode the evolution, interaction and dynamics of the spatial modes [17]. They encapsulate the projection of the system's dynamics onto the m-dimensional model space defined by Eq. (16).
In Fig. 9 we represent the modal time series (red lines) corresponding to the first four u x spatial modes for fault system A. The blue lines describe the time evolution of the binned potency released (the total potency released over the time interval ∆t). For the first time series, A 1 (n), the location in time of the amplitude jumps coincides with the time of large events (within the temporal resolution defined by ∆t). Moreover, we can show that the size of the amplitude jumps is proportional with the size of the event, i.e. ∆A 1 (n) = where P max (n) is the largest event in the n'th time interval. Therefore, if we can model the evolution of this mode, the time and size of the amplitude jumps will give us useful information about the time and size of the large earthquake events. Moreover, we notice that the higher order modes evolve on faster time scales, as is already noticeable from the faster decay of the temporal autocorrelations shown in Fig. 11.
The picture that emerges from this decomposition evolves gradually from the slow, coherent large scale motions, whose dynamics is captured by the most dominant POD modes, to the fast, incoherent small scale motions described by the higher order modes. In Section 7 we describe how a neural net can be used to process these time series in order to extract a low-dimensional nonlinear dynamic model with short-term predictive capabilities.

Geometry of the attractor in the POD basis
The POD modes provide an optimal basis for the phase space embedding expressed in Eq. (16). The phase space vectors in this reconstruction, X n = (A 1 (n), . . . , A m (n)), approximate the location of the system on the attractor at discrete time moments. Compared to the time-delay embedding, the POD basis unfolds the attractor with increased resolution as the embedding dimension increases. This property allows us to probe the structure of the attractor at different length scales by computing its correlation dimension. This measure was introduced by Grassberger and Procacia [14] to quantify the self-similarity of geometrical objects. It is based on the definition of the correlation sum for a collection of points X n in some vector space R m . This is the fraction of all possible pairs of points which are closer than a given ǫ in a particular norm. The basic formula is where we use the maximum norm, X = max |x i |, Θ is the Heaviside step function, N is the total number of data points, and R is the total number of reference points. The sum counts the pairs (X i , X j ) whose distance is smaller than ǫ. Note that data points within a time window 2W of any reference point j are not included in the correlation sum. An appropriate choice of W = 100 reduces the spoiling effect of autocorrelations from the correlation sums as suggested by Theiler [29]. In the limit of an infinite amount of data and small ǫ, the attractors of deterministic systems show power law scaling, i. e., C(ǫ) ≃ ǫ d2 , and we can define the correlation dimension d 2 by For each embedding dimension m, d 2 (m, ǫ) measures the local slopes of the correlation sum at different length scales ǫ, In For system A, in region III located at larger length scales between regions II and IV, the dynamics evolves on a low-dimensional self-similar attractor whose presence is detected by the plateau in the correlation dimension at d 2 (m, ǫ) ≃ 1.3. Note that for system A the scaling regime, ǫ → 0, is not reachable. The breakdown of scaling at the smaller length scales is dynamic in nature and is not due to the lack of good statistics. The closer we look at the system, the more degrees of freedom become visible and the dimensionality of the dynamics is higher than our largest embedding dimension. This also explains why the percentage of false nearest neighbors never drops completely to zero with increased embedding dimension. On the same figure, within the range of embedding dimensions that were numerically accessible, system B shows no signature of a large scale dimension collapse, in dramatic contrast with the behavior of system A.
The length dependent dimensions lead us to conjecture that we observe different subsystems on different length scales. From the self-similar and low-dimensional structure at the large length scales the structure of system A attractor crosses-over to a high dimensional, stochastic like structure when probed at very small length scales. Our strategy for building a model with good forecasting skills should be obvious now. The underlying idea is to approximate the large scale motions of system A through a low-dimensional embedding

Model Reconstruction and Short-Term Earthquake Forecasting
Using the POD decomposition we have identified a low-dimensional linear space in which system A evolves most of the time and we have reduced its dynamics to a small set of time series, A i (n), i = 1, . . . , m, describing the evolution of the system in this reduced linear space. We now face the problem of determining the underlying dynamical process from the information available in these time series. They are assumed to be governed by a nonlinear set of ODEs and our modeling approach relies on the ability to identify an approximate m dimensional model, that describes an explicit Euler approximation to the evolution and interaction of the spatial modes. To identify this nonlinear mapping we employ an artificial neural network (ANN). Generally, the neural network approach is used as a "black-box" tool in order to develop a dynamic model based only on observations of the system's input-output behavior [23,24]. In the learning process the network adjusts its internal parameters to minimize the squared error between the network output and the desired outputs. A typical learning method is the error back-propagation algorithm which is a first order gradient descent method [15].
All reconstruction results described here refer to the dynamics of the u x surface deformation modes for model A with D = 1.5. We found it difficult to identify a simple model having the structure defined by Eq. (24). In order to improve the model forecasting skill it is useful to enlarge the structure of the model to include information about the past history of the modes. While the best model structure is still a subject of investigation, an analysis of the time patterns present in the modal evolution provides a partial understanding of this result. As Fig. 9 shows, one essential feature of the modal dynamics is the presence of two time scales: within each earthquake cycle there are intervals of slow and fast motions with detectable quasi-regular behavior. This is reflected in the two-point autocorrelation functions (Fig. 11) defined as where the < · > denotes the time average. The resulting plots exhibit oscillatory behavior with a slow amplitude decay over a longer time scale indicating that the system has two correlation time scales. We have observed that providing the ANN with information about these long-term correlations of the modal dynamics produces nonlinear models with increased forecasting skill. To include this information, we have modified the structure of the model in Eq. (24) and replaced the modal coefficients A i with time-delayed vectors X i ,  (d − 1). The ANN output then provides a prediction of the mode amplitude A i at time (n + 1),

short-time memory
based on input information describing the past mode histories.
The ANN is trained using a standard back-propagation algorithm [15].
where i, j = 1, . . . , m denote the input modes, r, s = 1, . . . , K + (d − 1) are the components of the timedelayed mode vectors, and < · > is the average over the ensemble of input vectors. Due to the statistical independence of the modal coefficients, Eq. (17), the correlation matrix is to a good approximation block diagonal, i.e. K ir,js ≃ 0 for i = j (any modal cross-correlations are due to finite size effects). The eigenvalues and normalized eigenvectors of the correlation matrix, will now provide an optimal representation of the dynamics of the spatial modes: where We found that 99% of the variance of the input ensemble (m = 2, K = 6 and d = 6 model) can be represented by the first six temporal modes Ψ k , k = 1, . . . , 6. We have therefore chosen to approximate the input vectors presented to the ANN by their six dimensional POD projection B(n) = (B 1 (n), . . . , B 6 (n)) onto the space spanned by the first six temporal eigenvectors. Denoting by P X→B this projection operator, the structure of the model has now become: The best forecasting performance was obtained for an ANN with two hidden layers of 10 neurons each. The input to the network is the 6-dimensional projection B(n) and its 2-dimensional output, (A 1 (n+1), A 2 (n+1)), is the model forecast for the first two mode amplitudes at the next time step. Because the large events have long recurrence times, the sequential selection of the training set will be dominated by input-output pairs that only express the quasi-linear behavior between two consecutive large events: the ANN will learn a linear model and will fail to forecast the intrinsically nonlinear large earthquake events. Therefore, we have designed the training set to include enough information describing the dynamics around the scarce large events. "true" modes, describing the large scale motions, will never be exactly known. How can we then evaluate the model forecasting skill? A detailed answer is beyond the scope of this paper and our goal here is only to show that despite all these difficulties, and many others not mentioned here, our imperfect models can provide stable and robust forecasts. Our remarks here follow Smith [28], to which we also refer for a clear discussion of uncertainty in initial conditions, model errors, ensemble verification, and predictability in general.
For nonlinear systems, uncertainty in the initial conditions severely limits the utility of single deterministic forecasts. Internal consistency requires that all nonlinear forecast should be ensemble forecasts [28]. In this approach to forecasting, a collection of initial conditions, each consistent with the observational uncertainty, are integrated forward in time. When the system evolves on an attractor, the selection of the initial conditions is far from trivial. For a perfect model ( a model that has the right dynamics and whose phase space is identical with the system's phase space) the members of the ensemble should be restricted to live on the system's invariant measure (attractor). Obviously, this is not a priori known. Unrestricted initial conditions, consistent only with the observational uncertainty, will extend into the full phase space and generate overdispersive ensembles. When this is the case, the predicted probabilities will not match the relative frequencies as demonstrated by Gilmour [13]. Generally, all models, including ours, are imperfect and, therefore, a perfect ensemble and an accountable forecast method do not exists. Nevertheless, when the initial conditions are uncertain, ensemble forecasts are still required. A reasonable constraint would require the members of the ensemble to live on the projection of system's invariant measure (attractor) into the model phase space, although this choice is not necessarily optimal or even unique [28].
Assuming no a priori knowledge of the true attractor projected into the model phase space, we present in Fig. 12  time that coexist with regions of relatively short predictability time [8]. Significantly, the ANN has shown consistently the ability to improve its forecasting as the system approaches a large event. The predictability time is significantly larger than the microscopic expansion rate which is controlled by the maximal Lyapunov exponent. As hinted earlier, the predictability is scale dependent and the large scale motions are more predictable than the small scale motions.
The predictability also depends on the location of the system on its true underlying attractor. But unlike the growth of initial uncertainty in model phase space, the local nonlinear structure of the system's attractor controls how rapidly the small scale motions, not included in our model, progress to reach the large length scales. This expansion rate defines how rapidly the truth, red line in Fig. 12, diverges from the best guess model trajectory. This is the trajectory starting from the true system state projected onto the model phase

Conclusions
We describe a conceptual framework for modeling and forecasting the evolution of a large strike-slip earthquake fault. The approach relies on the detection of spatio-temporal strain patterns embedded in the observable surface displacements: no detailed knowledge of the fault geometry, dynamics, or rheology is required. There are many difficult problems that are currently under investigation. For example, we presently study how robust the large scale behavior is to changes in the model parameters. Preliminary results confirm the controlling role of the dynamic weakening and show that the deterministic character of the large scale motions is robust to large changes in this parameter. We are also studying the effects of correlated heterogeneities, smooth brittle to ductile transitions, and continuous transition from static to kinetic friction. We do expect that all these effects are averaging out the small scale dynamics and strengthen the deterministic structure of the attractor. Probably, the most difficult problems to address are the construction of good models from short data observations and the problem of spatio-temporal chaos. In the last case, dimensions and Lyapunov exponents become intensive quantities and we want to understand how the method scales with the size of the fault system. The current paper is concerned primarily with developing a methodology and the obtained results are preliminary. Continuing studies along the directions of this work may have a significant impact on the earthquake predictability problem.       of the variance in the u x , u y , u z data sets for system A we need only 2, 3, and 5 modes respectively.
In contrast, for system B we need 9, 27, and 42 modes. Generally, the surface deformations along the strike direction, u x , have the most efficient POD description.    We also notice a systematic bias in estimating the time to the next large event.