The GNSS stations used in this research were from the Geospatial Information Agency (BIG), or InaCORS, spread over Sulawesi, Kalimantan, West Nusa Tenggara, East Nusa Tenggara, Bali, and Java with 30 s sampling intervals including GLONASS and GPS satellites. A total of 16, 10, 10, 1, and 56 InaCORS stations are located in Sulawesi, Kalimantan, East Nusa Tenggara, West Nusa Tenggara, and Java with average distances of ∼96, ∼185, ∼87, ∼87 and ∼44 km, respectively. In total, 93 GNSS stations were used to examine the anomalies of the 2018 Palu Earthquake, as shown in Figure 1A. The condition of one-direction geometry between the epicenter and GPS receivers is shown in Figure 1B. This condition is important because it enables shallow LOS penetration with the CID wavefront.

A Slant TEC (STEC) value was obtained from the GNSS phase difference data conversion process between L1 (∼1.5 GHz) and L2 (∼1.2 GHz) using an ionospheric linear combination as follows: ∆ S T E C = 1 40.308 × f 1 2 f 2 2 ( f 1 2 − f 2 2 ) × ( ρ 1 − ρ 2 ) , where f ₁ and f ₂ are the carrier phase frequency, ρ 1 and ρ 2 represent the carrier phases from f ₁ and f ₂, respectively, and the phase difference ( ρ 1 − ρ 2 ) is expressed with the unit of meters (Heki, 2021). This research is focused on GPS satellites 21 and 25 due to the CID’s consistent appearance, as illustrated in Figure 2A. These GPS satellites are used, and both are visible from the studied region immediately before the mainshock (10:03 UT). However, in the reconstruction of tomographic 3D modeling, satellite GPS numbers 21, 25, and all the satellites in Figure 2B will be used.

The LoS will have an intersection point with the ionosphere when the GNSS satellites transmit microwave signals in both L1 and L2 frequencies. This intersection is called an ionospheric pierce point (IPP) at an altitude of ∼300 km, and its projection onto the earth’s surface is named a sub-ionospheric point (SIP). The SIP position is calculated using the formula developed by Klobuchar (1987). This research was conducted using satellites with no threshold or null elevation angles, as indicated by the SIP position, far from the zenith station.

After determining the STEC anomaly value, it is distributed to the cube block passed by the beam (LoS) by considering its individual length. The checkerboard test model area is made the same as the LoS distribution area, as shown in Figure 1. The trajectory plotting result of GPS and GLONASS satellites in the 2018 Palu Earthquake was found at 9,5—10,5 UT, as shown in Figure 2B. The GPS satellites which orbit during an earthquake are numbers 5, 12, 15, 20, 21, 24, 25, 29, and 31, while GLONASS are numbers 5, 9, 16, 19, 20, and 21.

The one-direction geometric condition between the epicenter and the GPS receiver was not found in this case. CIDs appeared after the mainshock and were detected by GPS satellites 21 and 25 with maximum amplitudes of 0.4 TECU (Cahyadi et al., 2022). This led to changes in STEC modeled using the reference curves obtained by fitting cubic polynomials of time to the vertical TEC (VTEC) (Ozeki and Heki 2010).

The input data for the calculation of the tomographic 3D block electron density anomaly are the derivative of the residual TEC (STEC) slant. The 3D blocks of ionospheric electron density were set up with the dimensions of 1 ° × 1 ° × 75 km over the Sulawesi, Kalimantan, Java, and surrounding areas, as shown in Figure 1A. The value of the STEC calculation results was distributed to the cube blocks passed through LoS by considering the length of each cube block with homogenous electron density. LoS penetrates multiple blocks, and the STEC residual can be expressed as the sum of the products of the penetration lengths and electron density anomalies of individual blocks, as shown in Figure 3.

In this case, the ionosphere is divided into voxel, therefore, the basic function is defined as follows: δ   ( λ , ∅ ,   h ) = { 1   i f   t h e   c e l l   i l l u m i n a t e d   b y   t h e   r a y 0   o t h e r w i s e   (1)

The penetration length is the distance between the two intersections of LoS with the block surface using simple geometric calculations in Eq. 2 from Fernandez (2004). Since the research area spans only a few degrees in latitude, the Earth is considered a sphere (its fattening is neglected) with an average radius. Figure 1B shows the geometry of LoS penetrating the blocks at an altitude of 300 km. ∆ S T E C i = ∑ j = 1 n A i j x j + e i (2)

The set of equation 2 for all LoS is written in a matrix form as y = A X + E (3) where y is a vector composed of STEC anomaly ( ∆ S T E C i ), A is a jacobian matrix composed of A j i , X is a vector composed of unknown parameters in the form of STEC anomaly in each block ( x ), and E is an error. Value of X can be estimated using the following equation: X = ( A T A ) − 1 A T y (4)

Although LoS is densely distributed, they do not penetrate all the blocks, specifically above the oceanic areas. This implies a certain constraint needs to be introduced to regularize the least-squares inversion. In this situation, a continuity constraint was applied, assuming that neighboring blocks have the same electron density anomalies with a specific allowance for the difference; and also use constraints around zero (Chapman function) to calculate neighbor voxel values based on altitude. By using this method, the voxel calculation value will be more realistic, following the observations by Muafiry and Heki (2020). One block normally consists of six neighboring blocks, namely up, down, north, south, east, and south, and all these pairs were added to the normal matrix as virtual observations (Nakagawa and Oyanagi, 1982). The non-juxtaposed lock pairs were not constrained because the tolerance corresponds to the virtual data’s “observation” error and the standard deviation of the actual differences between the adjacent blocks. It is assumed that the tolerance was 0.10 in 10¹¹ electrons/m³, equivalent to 1 TECU, or 10¹⁶ electrons/m², for a penetration length of 100 km. The influence of this value on the tomography results will be discussed in the next section. Coster et al. (2013) stated that the STEC observation error was also assumed to be 0.2 TECU, a few times as large as the typical error for differential GNSS VTEC measurements, which is consistent with the post-fit STEC residuals. The resolution of 3D tomography and its accuracy will be further discussed in this research.

A standard way to investigate the reliability of the 3D tomography solution is through the checkerboard resolution test. While conducting the checkerboard test, the real satellite/station geometry was kept while the STEC data were synthesized. The electron density anomalies were assumed at ±0.50×10¹¹ electrons/m³, as shown in Figure 4-upper panel.

Figure 4 (bottom panel) shows the distribution of the anomalies recovered using the synthetic data, suggesting that the spatial structures can be well resolved, with the recovered amplitudes accounting for about two-thirds of the input. Therefore, by comparing the two map views at different altitudes, the resolution at higher altitudes (200 km) was slightly worse than that at lower heights (100 km). This reflects better coverage and more penetration of LoS for lower blocks, as shown in Supplementary Figure S1 for altitude 100 km to 400 km. The figure also suggests that the resolution is higher and poorer above the land area and the ocean, respectively. The checkerboard test generally shows the high performance of the researched 3D tomography in the region of interest. Figure 4 shows that the resolution test results are poor at the left top of the bottom panel because the number of LoS distributions is scarce in that area. This is because only a few stations are used in the area, and the observed satellite objects are only GPS and GLONASS.

A second resolution test was also conducted to strengthen the tomographic modeling results by recovering a pattern consisting of a pair of positive and negative anomalies (0.6 × 10¹¹ el/m³) at low and high altitudes, respectively, on a neutral background (Figures 5A,C, as an input). The result, shown in Figures 5B,D, as an output reproduces properly, the assumed pattern of the positive anomaly is reduced to ∼2/3 the amplitude of the input model due to constraints at altitude of 300 and 400 km. Similarly, positive and negative anomalies in the latitude profile recovered properly, with only weak smears in the surrounding blocks not exceeding a few percent of the assumed anomalies.

The results of the 3D tomography modeling are presented in the time range from when the earthquake occurred (10.03 UT) until 10.20 UT in Figures 6(A-R). In Figures 6(A-F), the ionosphere is in a normal state, as indicated by the model’s color, without an anomaly in the ionosphere. At 10.09 UT (Figure 6G), CID began to appear to the south of the epicenter, marked by a positive red anomaly value until 10.16 Figure 6(n). The range of occurrence of ionospheric disturbances is indeed a curt phenomenon (only a few minutes) shown in Figures 6(G-O). This description is in accordance with the Cahyadi et al. (2022) study, which showed CID occurred ∼13 min after the mainshock. At 10.18 UT, the ionosphere has 226 started to return to normal as shown in Figures 6(P-R).

During CID, 13 min after the earthquake (10.16 UT), positive anomaly values were found at an altitude of 300–500 km near the epicenter, as shown in Figure 7. This is in accordance with the characteristics of variations in the ionosphere layer, where an altitude of 300–500 km is included in the F layer, which has the largest ionization and electron density compared to other altitudes. Fernandez (2004) stated that the electron density is mainly concentrated at a maximum between 200 and 500 km, following the Chapman profile. The CID shown by the tomography model matches the SIP position of the GPS satellites 21 and 25, which have succeeded in observing the TEC anomaly. We also estimated the vertical profile of electron density using the IRI model (https://irimodel.org/) on September 28, 2018 and the time when the CID occurred at 10.16 UT. The modeling results show that the highest electron density is located at an altitude of 320.4 km, as shown in Supplementary Figure S2. This supports 3D tomographic modeling, which is modeled at the height of the 3rd ionospheric layer with an altitude of ∼300 km (Chapman, 1931).

Another anomaly detected by the 3D tomographic model shown in Figure 6(n) (black box) could be caused by acoustic waves generated by several earthquakes that occurred before the mainshock (Mw7.5 at 10.16 UT). The USGS (https://earthquake.usgs.gov/earthquakes/) recorded several fairly large earthquakes occurring in the 3 h to 5 min before the main earthquake. The location of the epicenters of these earthquakes can be seen in the figure (Supplementary Figure S3). Meanwhile, the ionospheric anomaly in the 3D tomography Figure 6(n) (blue box) is likely caused by the acoustic wave of the main earthquake.

We performed a statistical test using residual values from the results of 3D tomography modelling. The test was carried out by comparing the distribution results of the post-fit STEC in five epochs, as shown in Figure 8. The post-fit residual showed much smaller dispersion, and its standard deviation is similar to the assumed STEC observation errors (0.2 TECU). The results before and after modeling showed a decrease in the standard deviation value. This decrease in value indicates that the results of the 3D tomographic disturbance of the ionosphere have been well modeled.

Marchetti et al. (2020) analyzed the presumed precursors that appeared before the 2018 Palu earthquake using the Swarm and CSES satellites, which observed changes in the ionosphere before the earthquake. These observations suggests an increase in electron density at night. We analyzed ionospheric disturbances for 40 days before the earthquake using GNSS-VTEC data to check whether there was a long-term preceismic anomaly. From the processing results, we follow the definition of ionospheric disturbance according to the observations made by Borries et al. (2015) in Equation (5), which compares the deviation value with the median value in the form of a percentage. The deviation value is obtained from the difference between the average TEC value and the median TEC value. The average TEC value with the median TEC can be seen in the table in the Supplementary Table S1. ∆ V T E C = V T E C a v e r a g e − V T E C m e d   V T E C m e d × 100 % . (5)

From the processing results, we analyze the magnitude of the deviation that occurs in the time series as in Supplementary Figure S5. It can be said to be an anomaly if the deviation >±30% (Danilov, 2022). The deviation was found at DOY 231—234 which is Dst Index value was checked from https://wdc.kugi.kyoto-u.ac.jp/dst_provisional/. The absolute value of the Dst Index shows relatively low values < 50 nT, which means that the geomagnetic activity is relatively calm, and the resulting anomaly is likely a precursor.

There is a slight difference between the results of processing anomalies before the earthquake in this study and Marchetti et al. (2020). In the 40 days before the earthquake, the electron density observed by the Swarm and CSES satellite anomaly occurred, which was marked by a positive anomaly, while the observation method using GNSS-TEC contained negative anomalies in the same range. This difference occurs due to changes in the electron density profile of the ionosphere layer from the use of different data. Muafiry and Heki (2020) observed a short-term anomaly indicated by a nearly balanced positive and negative anomaly, indicating that the anomaly was created by electron transport in different ionospheric layers.