Figure 4. Pre-attack seismic amplitudes as a function of time, and the seismic amplitude inter-station correlation. The progression of seismic amplitudes at three IS seismic stations during the 2 hours leading to the attack are shown in the left column, and the covariation of amplitudes during the last 30 minutes before the attack are shown in the right column. (a-c) The three-component averaged spectral amplitudes as a function of local time on October 7 are shown by the black curves. The AMZNI, YATR, and KZIT amplitude curves are computed by averaging over frequencies between 2.1 to 2.8 Hz, 5.5 to 6.3 Hz, and 2.7 to 3.3 Hz, respectively. The background noise curves in these frequency bands and their median are shown by the grey and blue curves, respectively. The red dashed curve indicates the initiation of the 2nd deployment stage, characterized by high amplitude seismic energy and strong inter-station correlations.  The amplitudes were computed for successive 300-second long windows with 75% overlap. (d-f) Amplitude covariation during the 2nd deployment stage, starting from 06:10 local time. Each panel shows the October 7 and background amplitudes in one station plotted against the amplitudes at another station. The station pair used for each panel is shown in the upper left corner. Dashed black and grey curves indicate linear fits to the data. The value of R, the correlation between the linear fit and the observations, is given in black for the pre-attack window, and in grey for the mean of the curves fitted for the windows covering the preceding Saturdays during 2021, 2022, and 2023.
    I used the background noise levels to estimate the likelihood that random fluctuations in the target frequency bands may have coincided at all three stations.  To test for significance I generated 106 1500-second-long time-series whose number of samples equals that of the observed time-series at each station. The random amplitudes follow a normal distribution whose mean and standard deviation are extracted from the background noise curves. I count the number of windows exceeding 150% of the background median, which is approximately the mean amplitude level observed during the last 1500 s before the attack, and the distribution of the average amplitude correlation coefficient for the random realizations obtained for each station pair (Figure S1). The tests suggest that the probability the amplitude of three random 1500-second-long signals (drawn from distributions characteristic of the background noise), will rise once or twice above 150% of the median pre-attack background level is 5% and 0.03%, respectively (Figure S1b and S1c). In reality, however, amplitudes during the 1500 s preceding the attack reached >150% of the median background levels during multiple 30-s windows in all three stations, establishing the statistical significance of the October 7 observations at a high level of confidence. 
    A further indication of the strong October 7 inter-station correlation is shown in panels 4d-f, which present the amplitude covariation in the 20 minutes leading to the attack.  The degree to which the October 7 correlation stands out above the background correlations is remarkable. For example, the pair AMZNI-YATR correlates at 92% on the day of the attack, but averages to only 10% during the previous Saturday mornings. The October 7 data at the other two station pairs are slightly less well-correlated, yet their pre-attack correlation level far exceeds the one expected given the correlation observed in the preceding Saturdays.  The statistical analysis suggests <10-6 probability the average of the station pairs' background amplitudes would correlate at >85% (Figure S1a), as was observed on October 7.
    The synthetic tests were conducted assuming the amplitudes are distributed as during previous Saturdays and in the target frequency bands observed on October 7. Because diurnal, weekly, and annual variations in the spectral shapes are common, it is important to estimate whether the strong correlations observed on October 7 (Figure 4d-f) can be observed on other days but at different frequency bands.  This likelihood is estimated by selecting the 2 to 8 Hz spectral maxima after smoothing the spectra using a 1 Hz window (close to the widths of the October 7 spectral bumps; Figure 2) and correlating the amplitudes recorded between 06:10 and 06:30 local time Saturdays starting from March 2021. The average AMZNI-KZIT, AMZNI-YATR, and KZIT-YATR correlation coefficients are 0.09\(\pm\)0.46, 0.05\(\pm\)0.46, and 0.14\(\pm\)0.41, respectively. These values are significantly smaller than the October 7 correlations. Moreover, only in 8 of the preceding 143 Saturdays did the average correlations lie between 0.6 and 0.7, however these days do not exhibit the October 7, 2023, high amplitude level and almost-monotonic increase with time towards the attack (Figure 4a-c).  This rules out the possibility that the October 7 correlations had emerged by chance.
Gaza traffic source location from tripartite array analysis
Seismograms from the station triplet may be used to locate sources giving rise to coherent surface waves travelling outwards from the Gaza Strip, for example due to ground impacts produced by heavy vehicles. These coherent events were located from tripartite array analysis applied to the three IS stations.  The IS seismograms are first band-pass filtered between 2 to 6 Hz and downsampled to a rate of 40 samples-per-second. Then, the vertical channel time delays are obtained by cross-correlating a sliding 30-second window with 50% overlap starting from 8000 s before the attack. For each candidate window at one station, I examine windows from a second station over intervals ranging from -Tr to Tr, where Tr  is the inter-station surface wave travel time, assuming the wave travels at 1.5 km/s. I then shift the horizontal channels at one station by the time-delay corresponding to the maximum vertical cross-correlation. Next, I search for the pair of rotations maximizing the horizontal channel cross-correlations. I rotate the two pairs of horizontal seismograms at \(10^{\circ}\) intervals over a range of possible angles. Because rotating by more than \(180^{\circ}\)causes the polarity to flip, one channel pair is restricted to rotate between 0 to \(180^{\circ}\), whereas the other pair is rotated between 0 to \(360^{\circ}\). The technique is similar to one that was previously applied in order to enhance catalogs of weak and emergent seismic signals attributed to tectonic tremors in Cascadia (Armbruster et al., 2014) and Mexico (Peng and Rubin, 2017). Similar to traffic induced noise, tectonic tremors give rise to signals containing multiple repeated narrow-band wavelets (Shelly et al., 2007, Inbal et al., 2018), which sometimes facilitate correlating seismograms recorded by stations distant from each other.
    For perfectly coherent arrivals, the sums of three time delays and the three rotation angles equal zero. Thus, the observed sums can be used to assess the detection quality. The tripartite array self-consistency criteria are defined as (Eisermann et al., 2018):
\[\begin{aligned} Q_{\Delta t} &= 1-\frac{\Delta t_{12} + \Delta t_{23} - \Delta t_{13}}{|\Delta t_{12}| + |\Delta t_{23}| + |\Delta t_{13}|}\\ Q_{\Delta\theta} &= 1-\frac{\Delta\theta_{12} + \Delta\theta_{23} - \Delta\theta_{13}}{|\Delta\theta_{12}| + |\Delta\theta_{23}| + |\Delta\theta_{13}|},\\ \end{aligned}\]
where \(\Delta t_{ij}\) and  \(\Delta \theta_{ij}\) are the time-delay and relative rotation between station j to station i, respectively. For high-quality detections, \(Q_{\Delta t}\) and \(Q_{\Delta\theta}\) are approximately equal to 1. To ensure the robustness of the locations, I retain time windows with \(Q_{\Delta t}\) and \(Q_{\Delta\theta}\) larger than 0.7, which amount to about 10% of the total number of windows.
    The time-delay measurements passing the self-consistency criteria are used for determining the source locations via a grid-search approach.  The difference between the distance from the i'th station to source and the distance from the j'th station to the source can be written as a function of the product between the wave speed and the observed time-delay: 
\[R_j-R_i=c\Delta t_{ij},\]
where \(R_j=||{\bf x_j^r}-{\bf x^s}||_2\), and \({\bf x^s}\) and \({\bf x_j^r}\) are vectors holding the surface coordinates of the source and the j'th station, respectively.  I solve for a uniform surface wave speed and for the surface location along the Salah Al-Deen Road, a major traffic artery crossing the Gaza strip from Rafah in the south to Beit Lahia in the north (Figure 1 and 5).  To find the best-fit wave-speed \(c\)  and source location \({\bf x^s}\), I set i=1 and perform a grid search over \({\bf x^s}\), restricting the search for locations along the Salah Al-Deen Road, and for \(c\) in the range between 0.6 to 3 km/s. I retain locations for which the misfit between the two observed and two calculated time-delays is smaller than 1 second (see inset histogram in Figure 5). These correspond to errors that are of the order of 2-3% of the source-to-receiver travel-times. Figure S2 presents an example of a misfit function for one detection, and the distribution of wave-speeds obtained from this analysis. The average wave speed is \(1.1\pm0.3\) km/s. To validate the wave speeds obtained via optimization, I compare them to speed of surface waves excited by known impacts (Figure S4), which give speeds between 1 and 1.8 km, slightly higher than the speeds recovered from minimizing Equation 2. This may be due to harder sedimentary rocks found between station YATR to JER and RMNI relative to the softer rocks typical of the coastal planes found between the Gaza strip and stations AMAZIN, YATR, and KZIT.
    The location resolution is obtained by analyzing the Equal Differential Travel Time curves (EDT ; Eisermann et al., 2015). In a uniform velocity model, the EDTs for each station pair are defined by hyperbolas that intersect at the source location, and whose width is determined by the uncertainties on the time-delay and velocity model. I assess the location resolution by perturbing the EDT curves according to the time-delay errors (estimated from the misfit distribution ; Figure 5 panel a) and wave speed range (estimated from the optimization results ; Figure S2).  The resolution for sources lying along the Salah Al-Deen Road is found to be about 1 km in the along-road direction, and about 12 km in the oblique direction (Figure S3).
    The optimal rotation angles in windows passing the detection criteria were used to infer the horizontal motion polarization orientation, and to constrain the surface source locations.  For Rayleigh and Love surface waves, the horizontal polarization angles are aligned in the radial and transverse directions, respectively.  Hence the motions excited by the passing wave are dependent on the azimuth of the vector pointing from the station to the source, commonly referred to as the back-azimuth.  For coherent arrivals, the back-azimuths from the three IS stations are expected to intersect at the source location. Thus, the polarization-derived back-azimuths may be used to determine source locations independently from the locations derived based on the time-delay measurements. 
    A map showing the pre-attack signal locations and optimal rotation angles is presented in Figure 5. Most of the sources are found to originate from a portion of the road extending between the city of Khan Younes in the south out towards the Al Zawayda and Nuseirat refugee camps, which are located just south of the City of Gaza. Fewer locations are resolved in the northern and southern most extents of the strip, near Beit Lahia and Rafah. These locations are compared with the locations inferred from the range of back-azimuth intersections.  The intersection of the range of back-azimuths pointing from each of the stations defines a polygon centered on the Gaza strip, encompassing locations that were independently derived from correlation-based time-delay measurements. This consistency lends further support to the array-based location approach.  At stations AMZNI and YATR the polarization directions are \(0\pm30^{\circ}\) and \(80\pm40^{\circ}\), respectively, and at KZIT it is \(280\pm25^{\circ}\).  Based on this observation, it is inferred that the AMZNI and KZIT stations are mostly sensitive to Love waves, whereas the YATR station is mostly sensitive to Rayleigh wave energy.
    The array analysis is applied to successive time windows, providing an opportunity to assess the space-time distribution of sources lying along the Salah Al-Deen Road. A time-series of the along-road locations during the last two hours before the attack began is presented in Figure 5b. Based on its space-time distribution, the activity can be grouped into two main phases, termed here the early and late Hamas deployments. During both phases, activity appears to have originated from Khan Younes and to advance mostly north towards the city of Gaza, and later towards Beit Lahia.  The first deployment is associated with slightly slower motion, whose exact speed is difficult to determine given the scattered locations. A line connecting the first detections occurring after 04:45 in each portion of the Salah Al-Deen Road gives a minimal velocity of 44 km/s, which is indicated in panel 5b. The later detections in that phase may be associated with motions in the road-oblique direction.
    The early phase consisted of activity that lasted for about an hour, and that was concentrated near Khan Younes and Al Zawayda. Between 05:45 to about 06:00, the rate of activity was diminished. Then, the activity resumed with fast advancement from Khan Younes north towards Beit Lahia and south towards Rahaf.  The second deployment persisted for about 30 minutes, until 06:30 when the attack began. The two phases identified in the time-space plot in Figure 5b correspond to intervals in the IS seismograms containing high-amplitude seismic energy.  The early deployment, which occurred between 04:45 to 05:45, is best observed on AMZNI and, to a lesser degree, also on KZIT (Figure 4a,c). The late deployment, between about 06:00 to 06:30, is manifested by an almost-monotonically increasing amplitude observed on the three stations during the minutes preceding the attack (Figure 4).
    To assess the significance of the October 7 event location results, I applied the tripartite array location scheme to 100 1-hour-long time windows occurring on Friday nights and early Saturday mornings local time.  By applying the same detection criteria used for the October 7 data to the 100 control windows, I find a regular hourly detection rate of approximately 5 events-per-hour. This rate is about 5 times lower than the one observed during the last 2 hours before the attack.  Thus, the spatio-temporal distribution of coherent events presented in this section, and the seismic noise amplitudes evolution with lead time presented in the previous section, are the seismic signatures of precursory activity leading to the October 7 terrorist attack.