Rupture Model of the 5 April 2024 Tewksbury, New Jersey, Earthquake Based on Regional Lg‐Wave Data Open Access

On 5 Friday April 2024, at 10:23 a.m. (EDT), a moderate earthquake of moment magnitude $Mw$ 4.8 (U.S. Geological Survey [USGS]) occurred at Tewksbury Township in northwestern New Jersey, United States (Fig. 1a). The ground motion reached the maximum intensity VI on the modified Mercalli intensity (MMI) scale, and over 5300 residents up to a distance of 25 km from the epicenter reported intensity V. Millions of people around the metropolitan New York City, about 65 km northeast from the epicenter felt intensity IV shaking. The earthquake is the largest instrumentally recorded event since 1900 in New Jersey (Stover and Coffman, 1993). The mainshock was preceded by a foreshock of $mbLg$ 2.2 on 14 March 2024 and followed by 150 aftershocks with magnitudes > 1.0 until 10 May 2024. About 7.5 hr after the mainshock, the largest aftershock occurred (⁠$Mw$ 3.7; USGS). The Community Internet Intensity Map received over 183,000 responses from the public though the “Did You Feel It?” (DYFI) website maintained by the USGS (Wald et al., 2011). The intensity is II around Washington D.C., ∼300 km southwest of the epicenter, whereas an intensity of III is reported around White Mountains, New Hampshire, ∼450 km northeast of the epicenter. Hence, there is a sense of higher intensity observations toward the northeast than the opposite southwest direction. The observed peak ground velocity (PGV) distribution of the mainshock shown in Figure 1b is consistent with the intensity reports; that is, higher PGV values are distributed along the east‐northeast direction.

Tewksbury Township is a leisurely rural hamlet located in the foothills of the New Jersey Highlands (Fig. 1a). Although the well‐known Ramapo fault system runs close to the area with a northeast trend, the focal mechanism suggests no correlation with the fault system (Fig. 1a; Aggarwal and Sykes, 1978). Although there are no large population centers within 10 km of the epicenter, an intensity VII is expected for the $Mw$ 4.8 event within a distance of ∼10 km from the magnitude–intensity relationship established in the central and eastern United States (Boyd and Cramer, 2014). However, even maximum intensity VI was hardly reported with significant property damages, but minor damages, such as cracked walls and objects being thrown from shelves. So far, no surface trace of faulting has been reported for this shallow earthquake of ∼5 km depth.

Therefore, clear focusing of ground motion toward the northeast as manifested by waveform data and PGV (Fig. S1, available in the supplemental material to this article, Fig. 1b) and relatively low near‐field strong ground motion excitation is puzzling; perhaps an unusual rupture process may have taken place. The site effects are also an important factor contributing to ground shaking. However, we found no correlation between known site effect distribution and the felt intensity around our study area (McNamara et al., 2015). In addition to these questions, reported focal mechanisms and our own analysis given in the next section indicate that the observed directional ground motion makes an angle of ∼45°–50° to the strikes of both nodal planes (NPs) of this oblique faulting event.

This report aims to investigate the detailed rupture process of the mainshock to answer the aforementioned questions. Because of the lack of local stations—only four stations within 100 km, we will employ the empirical Green’s function (EGF) approach and utilize regional Lg waves in retrieving relative source time functions (RSTFs).

The regional waveform‐based focal mechanism solution of the mainshock reported by the USGS showed oblique thrust faulting with a significant strike‐slip component at 7 km depth and $Mw$ 4.76 (Table 1). Global Centroid Moment Tensor (Global CMT) solution (Ekström et al., 2012) indicated $Mw$ 4.7 strike‐slip event at a fixed depth of 12 km (Table 1). To reconcile the two different focal mechanisms and focal depths reported, we carried out waveform modeling using the frequency–wavenumber integration method for the focal mechanism inversion (Saikia, 1994). We employed a crustal velocity model that has been used to locate earthquakes in the southern New York and northern New Jersey regions (Yang and Aggarwal, 1981). Observed waveforms from a dozen stations at 78–394 km distances are modeled. Period bands of 10–33 and 6.7–20 s were used to process waveforms depending on the epicentral distances. We determined focal mechanism parameters (strike, dip, rake, and seismic moment) using a grid search under a constrained deviatoric moment tensor (e.g., Kim et al., 2010). The best‐fitting double‐couple parameters are listed in Table 1, with the moment magnitude of $Mw$ 4.84. Focal depth is constrained by waveform inversion and by searching depths that maximize variance reduction (VR). The focal depth is fairly well constrained between 3 and 8 km, with peak VR at 5 km.

We followed the same procedure as the mainshock to constrain the focal mechanism and depth for the largest aftershock. We obtained $Mw$ 3.82 for the aftershock. The focal mechanism solution is nearly identical to that of the mainshock with a slightly higher thrust component (Fig. 1a; Table 1), and the waveform data fit best with a slightly deeper depth of 7 ± 2 km. Our moment magnitude (⁠$Mw$ 4.84) for the mainshock is about 0.08–0.14 magnitude units greater than those reported by USGS and Global CMT. The shallow focal depth of 5 km constrained in our solution can account for this difference.

We employed an EGF method to retrieve the RSTF of the mainshock (e.g., López‐Comino et al., 2012; Seo et al., 2022) at stations of distances from 132 to 555 km (Fig. 1b). The largest aftershock, which has a common focal mechanism and is located ∼1 km from the mainshock, is used as the EGF event (Fig. 1a, Fig. S1). We used data from 52 broadband stations covering an azimuth range from −149° to 66° (clockwise direction). The waveforms are dominated by the Lg wave (Campillo et al., 1984). We used Lg waves on transverse‐component records in the time window of 5 s before and 15 s after the arrival at an apparent group velocity of 3.6 km/s, filtered at 0.3–3 Hz.

The Lg RSTFs were obtained by aligning the waveforms of the EGF and the mainshock in the same time window starting from the origin times of each event rather than using phase arrival times or cross‐correlation times (origin‐time preserved RSTF; Han et al., 2024). Offset times in the RSTFs, which reflect travel‐time differences between the mainshock and EGF, were used to determine the relative locations of the mainshock rupture against the EGF location. We used projected‐Landweber deconvolution with constraints of positivity and causality, allowing us to capture the approximate start and end times of the source pulse (Bertero et al., 1997).

Although the Lg RSTFs do not have good depth constraints, they preserve source properties in the azimuthal direction (Gallegos and Xie, 2020; Han et al., 2024). The Lg RSTFs with the $VR([1−∑t(target(t)−EGF(t)*RSTF(t))2∑ttarget(t)2]×100%)$ greater than 50% were selected, normalized so that each integrated area is equal to the seismic moment, and stacked in 15° azimuth bins to enhance the robustness of their shapes (Fig. 2a).

The estimates of the rupture start and end points are (−790 m, 70 m) and (90 m, 590 m), and their corresponding times are 0.23 and 1.26 s, respectively (Fig. 2b). Therefore, the mainshock rupture propagated with a horizontally projected length of ∼1000 m, mostly extended toward the east‐northeast (∼59°) in ∼1.0 s (Fig. 2b). The rupture propagating toward opposite direction cannot be precisely estimated due to resolution limits imposed by the limited azimuthal coverage (Abercrombie et al., 2017).

Because we observed variations in the duration of RSTFs, we further utilized the spectral ratios of the mainshock and the largest aftershock as an EGF pair to identify rupture directivity in the frequency domain (Prieto et al., 2009). We stacked the spectral ratios for stations used in the RSTF analysis along three directions to represent directivity—forward, orthogonal, and backward of the mainshock rupture direction presumed in RSTF analysis (Fig. 2c). The higher and lower amplitudes along the forward and backward direction in high frequency, respectively, show strong directivity with unilateral rupture characteristics (e.g., Calderoni et al., 2017; Trugman, 2022).

Preliminary inversion was performed for two strike directions of 13° (NP1) and 109° (NP2) by nonnegative least‐squares minimization (Lawson and Hanson, 1974) while searching for the optimal $VR′$ and $tr$ that maximize the VR (Fig. S2). Results show that this simple 1D rupture model, constrained to propagate along the strike direction, cannot fully capture the observed RSTF variability across different azimuths (Figs. S2c,d). Therefore, we performed the inversion by varying $ϕr$ from 0° to 180° in 5° increments, considering the directivity associated with rupture propagation in the dip direction. We investigated the best‐fit $VR′$ and $tr$ for each rupture direction. The VR is significantly improved to 95.7% at $ϕr=55°$⁠, with its uncertainty range of 40°–65° estimated at 1.1 × optimum misfit (Fig. 3a,c). This $ϕr$ is quite close to the direction estimated from rupture start and end points (∼59°; Figs. 2b, 3a). When we used $VR'$ of 1.36 km/s (⁠$∼0.4×VS$⁠) and $tr$ of 0.4 s at $ϕr=55°$ (Fig. 3b), the slip distribution reveals that the rupture propagated forward for 1000 m in ∼1 s (Fig. 3d). A slip pulse at −900 m, as shown in Figure 3d, could be attributed to either backward rupture propagation or a small subevent.

To infer vertical components of the rupture direction, which is not observable from horizontal measurements of the Lg RSTFs, we assume the rupture direction is within the fault plane of two NPs. Then, we project the horizontal direction $ϕr=55°$ onto the NPs. We estimate true rupture velocity on the fault plane by accounting for the apparent velocity reduction by $cosδ$ in the dip direction, which is reflected in $VR'=0.4×VS$⁠. As the results indicate, the inverted rupture velocities for NP1 and NP2 are $∼0.5×VS$ and $∼3.1×VS$ toward the down‐dip and up‐dip directions, respectively. The steep dip of 84° results in an abnormally fast (super shear) rupture velocity for NP2, whereas the value for NP1 is acceptable.

We constructed 2D model space on a 4 × 4 km fault plane centered at (⁠$x0,z0$⁠) with a grid spacing of 100 m in both x and z directions. We implemented linear regularization for spatially adjacent values of $A(x,z)$ to avoid ill‐posed problems during the inversion (Everett, 2013). Based on the L‐curve test result, we selected a regularization parameter $λ=400$⁠. The optimal rupture velocities were found to be $0.55×VS$ for NP1 and $1.75×VS$ for NP2, at a risetime of 0.2 s.

Figure 4a shows slip distribution on NP1 with the rupture direction constrained from 1D inversion (blue arrow). High‐slip density is southeast to east, which is somewhat shifted from the expected direction along the blue arrow starting from the hypocenter. This might be due to the lack of station coverage in the southeastern quadrant (Fig. 1b) and the large model space of 2D inversion cells with 1600 elements (40 × 40). The rupture direction is robustly constrained by the 1D inversion (Fig. 3d), whereas the 2D inversion yields reasonable estimates of the main rupture area, clearly indicating a circular rupture of radius ∼1.1 km (Fig. 4a). A static stress drop of ∼6.6 MPa can be estimated from the source radius using $Δσ=(7/16)M0/r3$(Eshelby, 1957).

We note that high‐slip density on the fault plane is predominantly distributed in the down‐dip direction, and the area extending the farthest from the hypocenter corresponds to the forward direction of the 1D rupture. The projected EGF on the fault plane is situated at the periphery of the slip area (Fig. 4a, cyan circle), indicating that the largest aftershock occurred at the edge of the mainshock rupture. Synthetic Lg RSTFs are nearly identical to observations, indicating the robustness of our model (Fig. 4b). For NP2, the model shows overly extended slip distribution (>5 km) in the fault dip direction with anomalously fast (super shear) rupture speed (Fig. S3).

We inverted Lg RSTFs (Fig. 2a) using three progressively evolved models. (1) By utilizing the start and end times of the origin‐time preserved Lg RSTFs (Fig. 2a), we obtained a first‐order solution for the mainshock rupture that propagated for ∼1 km toward east‐northeast (∼59°) in ∼1 s. (2) Using the 1D rupture model (Fig. 3d), we obtained a dominant rupture direction of 55° from a linear slip distribution. The 1D inversion is robust due to its simple parameterization. (3) Through 2D modeling, we obtained a slip distribution on the fault plane (Fig. 4a), providing a realistic rupture dimension and rupture velocity. The slip distribution favors NP1 as the rupture plane and, hence, the rupture propagation in the down‐dip direction.

Assuming that the EGF and the mainshock occurred on the same fault plane, the relative horizontal location of the EGF to the mainshock suggests that the EGF is 0.9 km deeper than the mainshock. The depth difference is consistent with the result of the moment tensor inversion, showing ∼1 km deeper for the EGF than the mainshock. The EGF and mainshock relocated by the double‐difference method (Waldhauser and Ellsworth, 2000) using body‐wave phases indicated that the EGF is 1050 m east, 170 m south, and 830 m deeper than the mainshock (Figs. 2b, 4a).

If we consider an NP‐striking east‐southeast (109°) and dipping steeply to the southwest (84°; NP2) as the fault plane, the depth of the EGF becomes ∼2 km shallower than the mainshock hypocenter, which is inconsistent with the waveform modeling and relocation results (Fig. S3). The anomalously fast rupture velocity (⁠$∼1.75×VS$⁠) is obtained for rupture propagation toward the up‐dip direction that cannot be easily reconciled.

Figure 4c summarizes the rupture process on the fault plane based on our inversion result. It shows the horizontal rupture length of ∼1 km in the direction of 55°–59°. Based on the fault geometry of NP1, this corresponds to the rupture length of ∼1.23 km with a rupture velocity of $∼0.55×VS$ in the down‐dip direction.

The rupture progressing down‐dip (Fig. 4c) might have resulted in the lack of strong shaking exceeding MMI VI, which is usually expected for a moment magnitude of 4.8 events in the near‐field around the epicenter (Fig. S4). Only six felt reports show VI, which contrasts with the frequent VI region near the epicenter of the 2003 Central Virginia earthquake with a lower $Mw$ of 4.3 and a deeper focal depth of 10 km (Kim and Chapman, 2005). The high PGV values and the intensity values of IV broadly felt in the east‐northeast direction, where New York City is located, must have been due primarily to the mainshock rupture directivity at 40°–65° direction (Fig. 1b).

We summarize our proposed rupture model for the 2024 Tewksbury, New Jersey, earthquake as follows: the rupture initiated at 5 km depth on the fault plane trending north–south and dipping moderately to the east. It propagated down‐dip for ∼1 s with a speed of 1.87 km/s and in azimuth of 55°–59°. The mainshock rupture area indicates a circular rupture with a radius of ∼1.1 km that yields a static stress drop of ∼6.6 MPa. The proposed rupture model accounts for the lack of maximum intensity greater than V (MMI) at the epicenter and focused broad‐intensity IV areas and high‐PGV values toward the northeast.

Focal mechanism solutions from the U.S. Geological Survey (USGS; https://earthquake.usgs.gov/earthquakes/eventpage/us7000ma74/moment-tensor) and Global Centroid Moment Tensor (Global CMT; https://www.globalcmt.org/) were used. Seismic intensity estimated by “Did You Feel It?” (DYFI; https://earthquake.usgs.gov/earthquakes/eventpage/us7000ma74/dyfi/intensity) was used. Broadband waveforms and station metadata were obtained via the International Federation of Digital Seismograph Networks (FDSN) web services (www.fdsn.org/webservices). Seismic Analysis Code (SAC; https://ds.iris.edu/ds/nodes/dmc/software/downloads/sac/) and ObsPy (v.1.4.0; https://docs.obspy.org/) were used for data processing; and Basemap (v.1.4.1; https://matplotlib.org/basemap/stable/) and Generic Mapping Tools (GMT, v.6.3.0; https://www.generic-mapping-tools.org) for making geological maps. Past seismicity in the eastern United States was obtained from the Advanced National Seismic System (ANSS) Comprehensive Earthquake Catalog (https://earthquake.usgs.gov/data/comcat/). All websites were last accessed in May 2024. The supplemental material includes four figures.

The authors acknowledge that there are no conflicts of interest recorded.

This work was supported by the National Research Foundation of Korea (2022R1A2C1003006 and 2022R1A5A1085103). W.‐Y. Kim received support from the Consortium for Monitoring, Technology, and Verification under the Department of Energy National Nuclear Security Administration Award Number DE‐NA0003920. The authors thank Editor‐in‐Chief Keith D. Koper and Associate Editor Pascal Audet for their support and appreciate Rachel E. Abercrombie and one anonymous reviewer for their thoughtful comments and constructive reviews.