Multiyear day-by-day counts of nighttime temperatures were used to compute the daily temperature variations near the vicinity of the earthquake epicenter. Data near the ground surface were obtained from Tribhuvan International Airport, Nepal, through NOAA Surface Data Hourly Global Database. We analyzed surface air temperature and nighttime data for 2011–2015 close to the terminator time 0,500–0,600 LT to define the normal and abnormal state of the air temperature. The pre-terminator times have been found as one of the most sensitive indicators for buildup of thermal anomalies, because of the limited solar radiation impact during these hours (Ouzounov et al., 2006). The time series for January 1 to May 31, 2015, is shown in Figure 2C. We computed the residual values to distinguish between the current value and the multi-year year mean of the air temperature variability. The maximum offsets from the mean value reached near +5°C on April 20 and +4°C on May 5 (Figure 2C) with a confidence level of more than +2 sigma for all the observations from 2011 to 2015. This transient rapid increase in the surface air temperature is a little more significant than the remote satellite observations shown in Figure 2A, which agrees with the thermodynamic processes explained by the LAIC concept. (Pulinets and Ouzounov, 2011). To understand the day-by-day variability of the daytime and nighttime temperature during the earthquake events, we analyzed the hourly temperature near the epicenter in a new way. We used 3-h global surface air temperature data from the GEOS-5 assimilation (Figure 2B). The presence of ions in the atmosphere creates a possibility for water vapor molecules to join with these ions through the hydration process, which is different from condensation. The evaporation/condensation process, the phase transition of the first order, always occurs during the potential chemical equality. However, newly formed ions have different chemical potentials; what is to introduce should consider in the one-component approximation we introduce the correction to the chemical potential. The increase of the water molecules chemical potential indicates the strength of the nucleation process and can be used as an indicator of an imminent earthquake according to Boyarchuk et al., 2010 and demonstrated already within multiparameter analysis (Ouzounov et al., 2018a; Pulinets et al., 2020).

Long-term array observations of radon anomalies have been used for pre-earthquake studies (Fu et al., 2011, 2015, Zoran et al., 2012; Giuliani et al., 2013, Karastathis et al., 2019). Anomalous radon variations occurred a few days to a few months before the earthquake and are likely to be associated with earthquake development’s dilatancy and micro fracturing stages. Radon plays a crucial role in developing the lithosphere-atmosphere-ionosphere coupling (LAIC) concept (Scholz et al., 1973; Pulinets and Ouzounov, 2011; Pulinets et al., 2018) associated with the pre-earthquake process. During the 2015 earthquake sequences in Nepal, simultaneous measurements of soil radon-222 were recorded at the main campus of Jadavpur University, Kolkata, India. Deb et al. (2016) studied the precursory seismogenic radon. To observe coherent responses, the soil Rn222 concentration profiles were measured simultaneously at two nearby (~200 m apart) locations (Figure 3A), A and B, having uniformly clayey soil, within the Jadavpur University premises (22.5667^◦N, 88.3667^◦E). The soil moisture content at location A was lower than location B. It was expected that simultaneous recording of radon time series at these two locations would add confidence in identifying pre-seismic responses (Figure 3b). The 4-months time of their observations fortunately overlapped with the Nepal seismo-active period in 2015 (Figure 3). Their sensor, a solid-state nuclear track detector method, was used to detect the alpha-radiation in the radioactive radon gas. Two simultaneous 4-months of observations have been analyzed. During the observation period, four statistically significant anomalies (above the ±2σ level) were obtained on April 20, April 29, May 19, and May 29, 2015, simultaneously at both locations A and B. April 20 anomaly preceded the April 25 M7.8, and the April 29 was identified as precursors to the May 12. The May 29 radon anomaly is likely to be identified as a pre-seismic event to the M5.6 earthquake at Kokrajhar, Assam, on June 28, 2015 (Day et al. l, 2016). Despite that radon fluctuation on May 19 registered at locations A and B with >2σ significance, no earthquake occurred within a 1,000 km radius from Kolkata. Probably this anomaly is not associated with seismic origin but with geodynamics transition associated with the new Moon event on May 18, 2015, 1 day before the Radon anomaly occurrence. The connection between Lunar phases, geodynamics, and seismicity has been statistically established (Kolvankar et al., 2010). The overall interpolation of radon anomalies related to the Nepal 2015 earthquakes demonstrates that the reliable detection of radon anomaly due to seismicity requires simultaneous measurement of soil radon concentration by a broad network of radon/gamma sensors (Fu et al., 2011, 2015; Karastathis et al., 2019).

All thermodynamic models of the atmosphere, considering the phase transitions of water and latent heat fluxes, operate with latent heat (per mole or molecule) as a constant at a given temperature. It is equal to 0.422 eV per one water molecule. However, if one carefully deals with accurate data, one can often find that violations of the gas equation are observed (inexplicable additional variations in temperature, humidity, and pressure). We approached this phenomenon from a different approach different when studying the processes of ionization of the surface manner layer of the atmosphere by radon, the release of which from the Earth’s crust sharply increases at the final stage of preparation for a strong earthquake. The detailed calculations could be found in Pulinets et al., (2015). Here we discuss the final findings. As it turned out, atmospheric ions formed during the ionization of atmospheric gases by energetic alpha particles emitted by radon during decay become instantly hydrated. Hydration is not equivalent to condensation because the process takes place at any relative humidity level and does not require saturated steam. However, all the same, a phase transition of water molecules from a free to a bound state takes place, and, just as during condensation, latent heat is released. It turned out that the released amount of latent heat is more significant at a high rate of ion production and high concentration of ions than for ordinary condensation. How much larger? For earthquakes, the value ranges from 0.01 to 0.1 eV, i.e., in extreme cases, it can reach about 25% of the latent heat constant. (it’s the difference between the released latent heat during hydration and normal condensation). This difference we call the correction of the chemical potential of the water vapor in the atmosphere. We call it the atmospheric chemical potential (ACP) because, in the act of evaporation/condensation, the energy of the water molecule is equal to its chemical potential. To understand the day-by-day variability of ACP during the earthquake events, we analyzed 3 h global ACP data computed from the GEOS-5 assimilation (Figure 2B).

One of the main parameters we used to characterize the Earth’s radiation environment is outgoing long-wave-earth radiation (OLR). OLR has been associated with the top of the atmosphere (TOA) integrating the emissions from the ground, lower atmosphere, and clouds (Ohring G. and Gruber, 1982), and primarily was used to study the Earth’s radiative budget and climate (Gruber and Krueger, 1984; Mehta and Susskind, 1999).

Daily OLR data were used to study the OLR variability in the zone of earthquake activity (Liu, 2000; Ouzounov et al., 2007, 2018b; Xiong et al., 2010). An augmentation in radiation and a transient change in OLR was proposed to be related to thermodynamic processes in the atmosphere over seismically active regions. We can determine the atmospheric anomaly in Euler’s reference frame as a first approximation by subtracting the mean. The average can be defined at the same day each year, local time, and location over 11 years (i.e., more than one solar cycle). The advantage of this approach is its efficiency in the presence of long-term satellite observations. The OLR anomalous variations were defined as an E_index (Ouzounov et al., 2007) as the statically defined maximum change in the rate of OLR for a specific spatial location and predefined times. They have constructed analogously to the anomalous thermal field proposed by (Tramutoli et al., 2005, 2013).

The E_index was defined as statically estimated variability in OLR values for specific locations and periods: E _ Index i , j ( t ) = ( S ∗ ( x i , j  , y i , j  , t ) − S ∗ ¯ ( x i , j  , y i , j  , t ) ( / σ i , j       (1) Where: t = 1, K days, ( S ∗ ( x i , j  , y i , j  , t ) ) is the current OLR and S ∗ ¯ ( x i , j  , y i , j  , t ) is the computed mean of the OLR field, defined for multiple years of observations over the same location and same local time, σ i , j  is the standard deviation of   S * , and K is the total number of analyzed days. We use the Thermal radiation anomaly (TRA) index, a modified version of E_Index (Eq. (1)), TRA   Anomal y i,j ( t ) = ( A ∗ E _ Index i , j ( t ) ) / B (2)

A and B are regional coefficients, A—a mask, mainly defined by the regional seismo-tectonic patterns and TRA appearance frequency (A = 0.1–0.9). B normalizes each of NOAA 15 and 18-time series of OLR data to the same time coverage over 11 years averaging period (B = 1–3.5). The TRA indexes have been processed with additional preprocessing to avoid aliasing short wavelengths and spatial filtering based on a “minimum curvature” algorithm (Ouzounov et al., 2018a). We used OLR data from NOAA’s Advanced Very High-Resolution Radiometer (AVHRR) for satellite thermal analysis over Nepal. The results of the 2015 Nepal earthquake showed that there was a rapid increase in transient infrared radiation in satellite data in mid-March 2015. An anomaly is observed near the epicenter; it peaked around April 4–7, about 2–3 weeks before the M7.8 earthquake (Figure 4, Figure 6). Further analysis revealed another temporary OLR anomaly on May 2–3. The second M7.3 event occurred on May 12, 2015. (Ouzounov et al., 2015). Our results show that long-wave radiation signals associated with earthquake processes were observed near the epicentral regions several days before the corresponding earthquakes. TRA hotspots appeared quickly, remained in the same regions for several hours, and then quickly dissipated. During March-May 2015, about TRA 15 anomalies were detected around the Nepali M7.3 epicentral areas (Figures 5,6). With red dots showing the anomaly locations, the date is given in the text, and a circle shows the confidence area. The centers TRA are computed based on Eqs. (1), (2). The energetic threshold for determining the anomalous patterns was >2.5 sigma STD. TRAs on April 5, 23, and 29 are within the maximum level for the entire period inside the entire region. The possible reason for several TRA’s is the activating of gas releases over Nepal and Central Asia. The triggered ionization inside the ABL generates zones of “thermal-bursts.” The cross-section of several TRA areas probably indicates the spatial clustering of the degassing along joined tectonic lineaments. The appearance of TRA anomalies almost disappears after the May 12 Earthquakes. The spatial distribution indicates a large appearance area, but the distance from the epicenter is allowed inside the Dobrovolsky (Dobrovolsky et al., 1979) and Bowman et al., 1998 area (R = 10^0.43M/R = 10^0.44M), which is about 2,250 km according to Dobrovolsky. This rapid enhancement of radiation could be explained by an anomalous flux of the latent heat over the area of increased tectonic activity. (Pulinets and Ouzounov, 2011, 2018; De Santis et al., 2017). Analogous observations were observed within a few days before the most recent major earthquakes in Japan (M9, 2011), China (M7.9, 2008), Italy (M6.3, 2009), Samoa (M7, 2009), Haiti (M7, 2010), and Chile (M8.8, 2010) (Ouzounov et al., 2011b, c; Pulinets and Davidenko, 2014).

For our analysis of ionospheric data, we used two sources of information: global maps of the total electron content in IONEX format provided by JPL and a time series of calculated vertical TEC of two GPS receivers in the region (lhaz and lck3). We also controlled the solar-geophysical conditions to purify the data from the possible solar-geomagnetic activity effects. Considering that we deal with the equatorial ionosphere, the primary source of geomagnetic activity was the equatorial geomagnetic index Dst provided by Kyoto World Data Center for Geomagnetism (http://wdc.kugi.kyoto- u.ac.jp/dstdir/index.html).

The epicenter of the M7.8 earthquake was at the outer slope of the northern crest of the Equatorial Ionization Anomaly (EIA). To detect the ionospheric precursors, it is necessary to carefully analyze the geophysical conditions around the time of the Nepal earthquakes. For this purpose, the Solar-geophysical conditions during April- May 2015 is shown in Figure 7, where the Solar radio flux F10.7 were analyzed. To put all parameters together (which are expressed in different values), all parameters were normalized. Where F10.7—it is Solar electromagnetic radiation on the Wavelength 10.7 cm: GEC—it is the Global ionospheric content which is the sum of all values in the IONEX table; REC (R5) - the Regional Ionospheric Content—is the sum of all values from the IONEX index within the 500 km radius circle; and REC (R10)—t is the Regional Ionospheric Content - is the sum of all values from the IONEX index within the 1,000 km radius circle.

In Figure 7B we plotted the Dst (equatorial geomagnetic) index for the same period. One can see that we are dealing with a moderately disturbed period. From the Dst we can conclude that we have three moderate (with Dst ~ −80 nT) geomagnetic storms: before the first earthquake and one started during the second earthquake. The other disturbances we can classify as small, not exceeding −30 nT (continuous line, bottom panel of Figure 8). The period is also characterized by the essential variations of F10.7 (min 100, max 155). Normalized F10.7 is shown as the green line in the upper panel of Figure 7. The F10.7 variations should be eliminated in the dTEC variations (defined in Eq.(3)) because they contribute to the running mean calculations. The first attempt to reveal the pre-earthquake variations having disturbed conditions is considering the difference between the Global TEC and Regional TEC (blue and red lines in Figure 7). As we see from the upper panel of Figure 7, the Global TEC follows the F10.7 with a delay of nearly 2 days (Afraimovich et al., 2008; Marchitelli, V et al., 2020). We exclude the difference in the vicinity of day 101 (April 10), which is a magnetic storm day. The most suspicious days are 111 (April 21) and 114 (April 24), when we may expect negative and positive anomalies. The differential GIM maps were computed for these days to check this supposition, which is presented in Figure 8. We see the strong negative (in crests) anomalies on April 21 and a powerful positive anomaly on April 24. The strength of the equatorial anomaly should change. Considering that the epicenter is inside the northern crest of EIA, we will express the strength of the TEC anomaly at the Northern crest to the TEC in the trough of EIA (Figure 9). We also calculate S - the strength of EIA on undisturbed days as a relation between the TEC value in the crest of anomaly and in trough between crests. On day 111 (April 21), the EIA completely disappears which gives its strength less than 1: S = 0.98. S- It is the strength of equatorial anomaly (relation between the TEC value in crest of anomaly and in trough between crests). On day 114 (April 24), the anomaly strength is 1.69, and on disturbed days before the earthquake—day 98 (April 07) and after the earthquake - days 117 (April 27) and 120 (April 30), the strength varied from 1.22 to 1.39.

From the historical data of the ionosphere, monitoring has elaborated a conception of the precursor mask, which generalizes the pattern of ionosphere parameter variations (foF2 or GPS TEC) a few days before the earthquakes (Pulinets et al., 2007, 2018). Δ T E C = 100 ⋅ ( T E C − T E C m ) / T E C m (3) Where TEC—is the Total Electron Content in TECU, where TECU = 10¹⁶ el/m², TEC_m 15-days running average value of TEC. The precursor mask was developed based on the analysis of strong (M ≥ 6) earthquakes in Greece and Italy (Pulinets and Ouzounov, 2018; Davidenko and Pulinets, 2019). We used this approach to analyze the dTEC variations for the lhaz GPS receiver around the two Nepal earthquakes presented in Figure 10. The horizontal axis is the Day of Year 2015 of the mainshock designated by zero. Vertical axis—is the local time (LT), and VTEC deviation is color-coded in percentages relative to the 15-days mean. We can see that the positive anomaly appears near 8 PM and exists continuously up to 4 AM the next day. In the present figure, the anomalies appear 3 days before the mainshock. However, further analysis for other earthquakes revealed that such anomalies might appear up to 5 before the mainshock. We can see similar night-time positive deviations before both strong earthquakes lasting more than 1 day and repeating every day at the same local time. The only difference with the pattern presented in Figure 10 is that the favorable variations start at 6 PM and finish at 6 AM. It can relate to the different terminator times for Greece and Nepal (Zolotov, 2015).

One of the possible experimental techniques that can monitor the ionization’s perturbations within the lower ionosphere uses Very Low Frequency/Low Frequency (VLF/LF) probing. Waves in the VLF frequency range are trapped between the lower ionosphere and the Earth and are reflected by the D region at an altitude of ~65 km during the daytime and ~85 km during nighttime. The received signals contain information about the reflection height’s region and its variability (Barr et al., 2000). The propagation of sub ionospheric VLF/LF (15–50 kHz) signals from navigation or time service transmitters over distances of thousands of kilometers (with low attenuation ~2–3 dB per Mm) enables remote sensing over large regions of the upper atmosphere in which ionospheric modifications lead to changes in the received amplitude and phase of the signals. The regular monitoring of many years at the Far East network, which operates in a highly seismic zones such as Japan and Kurile Islands, has established a statistical correlation between anomalies of the VLF/LF signal parameters in the nighttime and earthquakes with М ≥ 5.5. When observed, the possible time intervals of seismic-related phase and amplitude anomalies are about 1 week before an earthquake and 1 week after the event (Rozhnoi et al., 2004,2009,2015; Maekawa et al., 2006; Hayakawa et al., 2010, Biagi et al., 2011; Hayakawa 2015). Anomalies for such earthquakes were found in 20–25% of all cases. However, for strong earthquakes (M > 6.8), anomalies VLF/LF were observed in 60–70% of the earthquakes (Rozhnoi et al., 2013).

The VLF/LF receiving stations deployed both in Europe (Eastern Europe, Sheffield, Graz), the Far East (Kamchatka, Sakhalin, Kuril Islands), North America (Orange, CA) and Asia (Bishkek, Varanasi) are equipped with the UltraMSK receivers (http://ultramsk.com/). All the stations simultaneously receive the amplitude and phase of MSK (Minimum Shift Keying) modulated signals with fixed frequencies in a narrow band 50–100 Hz around the central frequency and adequate phase stability. The receivers can record signals with time resolutions ranging from 50 msec to 60 s. For our purpose, we use a sampling frequency of 20 s. The VLF/LF observation network is shown in Figure 11A, and the epicenters of earthquakes with M > 5 for the last 3 years. The analysis reported in this paper about earthquakes in Nepal in April-May 2015 is based on these data recorded by the VLF/LF stations in Bishkek (KGZ) and Varanasi (VAR).

Figure 11B, shows the relative positions of our observing stations and transmitters VTX (17.0 kHz) in India, NWC (19.8 kHz) in Australia, and JJY (40 kHz) in Japan, together with the positions of the epicenters of earthquakes in April-May 2015 in the region under analysis. The areas of earthquake preparation where precursors can be found (Dobrovolsky et al., 1979) are shown by the pink circle for the first strong earthquake on April 25 and the yellow circle for the second earthquake on May 12, 2015. The station in Varanasi (VAR) began a regular reception in April; therefore, analysis for both stations were made from the beginning of April. Two wave paths—NWC-KGZ and JJY-VAR pass directly above the epicenter of the Nepal earthquake. The signal from the JJI transmitter received at the Varanasi (VAR) station also passes above the epicenter, but it was not used for the analysis because it is very noisy. Two sub-ionospheric paths (VTX-KGZ and VTX-VAR) are outside the epicenter but inside the earthquake preparation area. We used a residual signal of amplitude calculated as the difference between the actual signal and the monthly averaged signal for our analysis.

The last was calculated using the data from undisturbed days. Since VLF/LF signals are very stable during daytime and are unaffected by any force except X-rays emitted during solar flares, the analysis was made only for nighttime. The results of these analyses are shown in Figure 12A. For validation of the results, we analyzed the signals from the same NWC and JJY transmitters propagating far away from the seismic zone, the significant negative nighttime amplitude anomalies for the four paths crossing the area where the possible precursors of the earthquake can be found 4–5 days before the first earthquake while the signals in the “aseismic” paths vary insignificantly (Figure 12B). Then after several days, when the signals were quiet, the second series of anomalies can be seen. These anomalies continue after the earthquake during the period of aftershock activity.