Insight into skywave theory and breakthrough applications in resource exploration

Abstract Skywave refers to the electromagnetic wave reflected or refracted from the ionosphere and propagated in the form of a guided wave between the ionosphere and the Earth's surface. Since the skywave can propagate over large distances, it has been widely used in long-distance communication. This paper explores and demonstrates the feasibility of skywave for deep resource and energy exploration at depths of up to 10 km. Theoretical and technical advancements were accomplished in furthering skywave applications. A new solution method based on Green's function has been developed to study skywave propagation in a fully coupled lithosphere-air-ionosphere full space model. For the first time, the model allows one to study skywave distribution characteristics in the lithosphere containing inhomogeneity such as ore deposits or oil and gas reservoirs. This model also lays a foundation for skywave data processing and interpretation. On a parallel line, we have developed a multi-channel, broadband, low-noise, portable data acquisition system suitable for receiving skywave signals. Using the skywave field excited by a high-power fixed source located in central China, actual field surveys have been carried out in some areas in China including the Biyang depression of Henan Province. The initial results appear encouraging—the interpreted resistivity models prove to be consistent with those of seismic exploration and known geological information, and the exploration cost is only ∼1/4 to 1/10 that of seismic surveys. These initial successful applications of the skywave theory lay a solid foundation for further verification of the new method.


INTRODUCTION
Skywave refers to electromagnetic (EM) waves propagating in the waveguide formed by the ionosphere and the Earth's lithosphere [1][2][3][4][5][6][7]. Due to the large contrast in the electrical property between the ionosphere and air and between the lithosphere and air, the EM waves are reflected back and forth between the ionosphere and the lithosphere as they propagate within the air waveguide. The skywave was first applied in radio communication [8][9][10][11]. Since skywave propagation is not affected by the curvature of the Earth, long-distance communication can be achieved (Fig. 1a). In traditional radio communication, high-frequency (>3 MHz) EM waves are usually used for efficient information transmission. Since the air has very high resistivity (10 13-15 m), even high-frequency EM waves can travel very far in the air, making remote communication with high-frequency radio waves possible. In the middle of the 20th century, with the widespread deployment of submarines, demand emerged for communication between land and underwater vehicles. However, because of the rapid attenuation of highfrequency EM waves in seawater, underwater vehicles could not receive signals transmitted from the land. To tackle the problem, the possibility of using low-frequency skywave for long-distance communication was explored [12][13][14][15][16][17][18][19][20]. For example, the United States established two super-low frequency (SLF, 30-  100 m. The former Soviet Union set up a larger SLF transmitting station (ZEVS) in the high-resistivity area of the Kola Peninsula, and successfully realized the communication between land and underwater submersibles at a distance of 8000 km and a water depth of 100 m. To further increase the depth of underwater communication, scientists from all around the world have also explored the extremely low frequency (ELF, 3-30 Hz) band [21][22][23][24]. Due to the difference in electrical properties of underground resources, the EM method has become the key technology in resource exploration [25][26][27][28][29]. The controlled-source audio-frequency magnetotelluric (CSAMT) method has been one of the most important methods used for resource exploration in the earth [30]. The CSAMT method has the advantages of high resolution, strong antiinterference ability and high operational efficiency, and has gradually become the mainstream method for resource exploration worldwide. However, the CSAMT method has the limitations of low transmitter power, small signal coverage and shallow sounding depth. The actual depth of investigation of the CSAMT method depends on the resistivity structure of the survey area (for the geological environment of China, the depth of investigation of the CSAMT method is usually <1.5 km). With the gradual depletion of shallow resources, demand has become increasingly high for exploration methods with larger depth of investigation [31,32].
Skywave bears great potential in resource exploration in the earth due to its large penetration depth into the lithosphere. Russian scientists have applied ZEVS signals to geophysical studies, earthquake prediction and space physics areas [33][34][35][36]. However, skywave applications for resource exploration have been scarce. The first challenge in applying skywave to resource exploration is the establishment of a proper physical model. In skywave communication, the Earth's lithosphere is assumed to be a nearperfect conductor so that it is largely excluded from the model. However, for resource exploration, this model becomes inadequate as the heterogeneities in the lithosphere now become the main target of interest. Thus, a fully coupled lithosphere, air and ionosphere model must be used (Fig. 1b). The second challenge is the lack of proper data acquisition equipment. The required equipment must be able to acquire multi-channel and wide-band data, and must have extremely low noise level. Meanwhile, the equipment must be sufficiently cost effective for large area exploration. However, no suitable instrument was available on the market to meet both the technical challenges and the economic requirements. Therefore, new equipment had to be developed. Moreover, literature on the data processing of skywave signals was rare. In this paper, we report some of the breakthroughs we have made in recent years in skywave theory and equipment development for resource exploration. We also present a field example to illustrate the application of skywave theory to oil and gas exploration.

Solution method
To study skywave propagation in the lithosphere, a fully coupled lithosphere-air-ionosphere model is first established. This model is different from the traditional skywave model for communication in that the communication model excludes the lithosphere completely, whereas the fully coupled model enables the study of the lithosphere, including any inhomogeneity within. Also, the fully coupled model allows for the study of skywave propagation and attenuation over large distances. The model lays a foundation for skywave survey design and data interpretation.
For skywave propagation over a large distance, a spherical Earth model should be used to accurately account for the Earth's curvature. However, for propagation distances on the order of a few hundred kilometers, a planar model can be adequate. We derive the planar full-model skywave theory by expanding the half-space CSAMT model to include ionosphere and lithosphere. The solution is achieved by the integral equation method based on the Green's function, and derived using the R-function and layer-matrix method [37,38].

Electromagnetic field characteristics
By using the R-function and layer-matrix method, we calculate the Green's function for the horizontal lithosphere-air-ionosphere model over the territory of China. The Green's function is used as the kernel of the integral equation method for three-dimensional (3D) modeling [39]. Figure 2 Natl Sci Rev, 2021, Vol. 8, nwab046

Equipment
For traditional CSAMT or magnetotelluric (MT) applications, different kinds of receivers are available on the market for EM exploration [40]. However, none of the receivers were designed for acquiring skywave signals. Their bandwidths and noise rejection capabilities appear inadequate for skywave applications. In this regard, we have developed a new receiver system for our needs. The specifications of the new broadband and low-noise system are listed in Supplementary Tables 1 and 2. Figure 3 shows the main components of the system. Field experiments proved that the newly developed acquisition system was able to acquire skywave fields for resource exploration at a depth of up to 10 km. It is worth mentioning that the instrument may also be used for traditional CSAMT and MT surveys.

VALIDATION OF THE MEASUREMENTS
Prior to any field survey, we conducted experimental field studies to validate the skywave measurement against the well-known CSAMT and MT methods. This would establish a solid basis for performing massive skywave surveys. It is helpful to first understand how the different methods respond to a layered Earth model. For survey parameters similar to those used in this study, Supplementary  Fig. 3 shows that, over all frequency ranges of interest, both the skywave and the MT responses are characterized by plan wave propagation, whereas the CSAMT method demonstrates the near-field feature at the lower frequencies due to the short transmitter-receiver spacings. This result suggests that the skywave data in the modeled transmitterreceiver distance range can be interpreted as the MT data. Next, the skywave method was validated against the MT method over an area with low levels of EM noise. As Supplementary Fig. 4 shows, over the overlapped frequency band, the skywave sounding curves compare well in general with the MT response curves, with the skywave response being more stable than the MT response. The skywave method was further compared to both the CSAMT and MT methods in a highly noisy area. As Supplementary Figs 5 and 6 illustrate, the skywave data show good consistency over the entire frequency range, whereas the CSAMT or MT data demonstrate uninterpretable variations or outliers due to environmental noise.

CASE STUDIES
Two skywave surveys were carried out in China. The first study area was in Chongqing, southwest China, ∼700 km away from the transmitter. The second test area was in Henan Province, central China, ∼300 km away from the transmitter. In both cases, a total of 30 sets of receiver systems were deployed.
In the Mingyuexia oil and gas area, the skywave exploration result was compared with seismic images. The skywave method was able to outline the largest regional fault, from a resistivity anomaly, as explained by the seismic reflection image. Both electrical and seismic reflection images show the typical chevron anticlines of the Huaying Mountain foldand-thrust belt thin-skinned structures. The interested reader should refer to [41] for more details. This paper shall focus on the Biyang depression case study.
The study area is part of the Biyang depression in the Henan Province. Figure 4a and b show the geological structure of this area as determined by previous studies. The Biyang depression is a fanshaped depression with a southeast-northwest orientation. The major basal structures are normal faults, named F1 through F4. The primary basin fault F1 spans 5 km horizontally and is located at the southern extent of the survey line. Part of the Wangji-Xingzhuang tectonic zone (WXTZ) is generally adjacent to the northern gentle slope zone (located at distances of 25 km to 33 km in Fig. 4b) along the boundary of Fault F3 (marked in Fig. 4b). Most of the hydrocarbons discovered in the Biyang depression are distributed in the layer Eh3 (marked by the pink color in Fig. 4b).
In this survey, a transmitter in central China was used. The transmitter consists of a north-south oriented antenna that is 60 km long and located on the surface with both ends grounded in the earth. A sinusoidal current was injected into the earth through the grounded wire, and the oscillating current created EM waves that propagated outwards. During the survey, 24 frequencies were transmitted in the range of 0.1 Hz to 256 Hz (see Supplementary Table 3).
The length of the sounding profile was ∼34 km, with an azimuth of 174 • , and the spacing between skywave receivers was 80 m (Fig. 4a). The signal was collected in a scalar CSAMT format in this study, where the electric field Ex and the magnetic field Hy were recorded. The data processing was also similar to that in CSAMT [42]. The apparent resistivities were used to perform a laterally constrained inversion [43,44].
The inversion model obtained from the skywave data (Fig. 4c) shows good agreement with the basin geometry. The resistivity model suggests that the study area can be divided into three parts from south to north: an outcrop zone in the south of the basement rock (I), a depression zone in the middle (II) and a gently sloping zone in the north (III). In addition, there are two areas with complex resistivity structures along the survey line: the Anpeng-Yinzhuang tectonic zone (APTZ) and the WXTZ. Combined with the stratigraphy of the study area, a comprehensive interpretation is provided in Fig. 4c.
A comparison of Fig. 4b with Fig. 4c demonstrates that the resistivity structures determined by inversion of the skywave data are consistent with the known fault structures. In particular, the hydrocarbon layer can be distinguished by the resistivity values in the range of 1-30 m, thereby providing important guidance for subsequent drilling. The resistivity model also shows a deep low-resistivity zone toward the northern end of the profile that was not present in the geological mapping.

DISCUSSION AND CONCLUSION
To apply the extremely low-frequency skywave to the exploration of deep Earth resources, the theory on the skywave wave propagation has been extended to a fully coupled lithosphere-air-ionosphere model. The R-function and layer-matrix method are used to solve the problem of skywave based on the 3D integral equation method. The inclusion of displacement currents in the model allows the accurate study of skywave propagation within the lithosphere-ionosphere waveguide on a global scale. For practical application of the skywave theory, we have developed a new multi-channel, broadband, low-noise and portable data acquisition system suitable for field surveys. Field studies have been successfully carried out in the Mingyuexia oil and gas area of Chongqing and the Biyang depression of Henan Province, and the resistivity models were derived from the skywave signal for depths up to 10 km. The obtained lithospheric models are consistent with those of seismic exploration and known geology, which validates the skywave theory and technique developed in this paper.
The skywave method has a few important advantages over the MT and CSAMT methods. By using a much more powerful transmitter than in the CSAMT method, skywaves can be detected at a distance of thousands of kilometers from the transmitter. These waves may be processed as plane waves, as in the MT method, greatly simplifying the data interpretation for the structure of the Earth's crust. Because the skywave method uses human-made sources, the signal-to-noise ratios are more controllable than for the MT method and the data quality is generally better.
Our theoretical research and field case studies show that extremely low-frequency skywave exploration may become one of the practical methods for resource and energy exploration within the upper 10 km of the crust. The cost of this method is significantly less than the seismic exploration method. In addition, the skywave method has the potential of discriminating oil or gas from water based on the resistivity contrasts of the fluids. When combined with seismic methods, the skywave method can be effectively applied in exploration of deep oil and gas.
The field examples show that the skywave signal can be stronger or even much stronger than background noises. Given the broad signal coverage, a large area (e.g. the whole of China) can be surveyed simultaneously as long as a sufficient number of survey stations cover the survey area. This can greatly reduce survey costs. Offshore application can be another area worth exploring. Extremely low-frequency skywave may supplement or supersede the traditional marine EM methods for offshore exploration.
As a resource exploration method, the extremely low-frequency skywave method warrants further study in some aspects. For example, few studies are available to address the response characteristics of the skywave in the near-field, far-field and waveguide areas. The applicability and limitations of the method for mineral exploration in these different areas call for more theoretical and field studies. Another area of study is the intermediate zone effect that describes the influence of lithospheric heterogeneities between the transmitter stations and survey areas. In this study, we have assumed that ground waves are highly attenuative and can be ignored in the acquired skywave data. Nevertheless, a detailed study of the intermediate zone effect may yield new insights into the skywave method for mineral exploration. Also, for use in areas very far from the source, it is necessary to develop the theory of a spherical layered model.

R-function method
For the R-function method, the upper half-space includes two layers: air and the ionosphere. The lower half-space is an N-layer lithosphere, and the N-th layer extends downward infinitely. The integral solution of the model is derived, and the numerical solution of the lithosphere-air-ionosphere model is obtained by using the Fast Hankel transform.
The full space model consists of N + 2 layers, including the ionosphere, air and N lithosphere layers in the lower half-space.
The field of the p-th layer is X P = c p e u p z + d p e −u p z , p that need to be solved are unknown, where p = −1, 0, 1, 2, . . . N. After getting these parameters, we can calculate the EM field. In the solution, the R-function is defined in the lower half-space, p e u p z . The function at the interface of each layer can be obtained by recursions. It can further be proved that Since R 1 (0) and R * 1 (0) are known, the ratio of the field in the first layer of solid rock X 1 (0), its vertical derivativeX 1 (0), V 1 (0), and its vertical derivativeV 1 (0) are determined.
The unknown parameters in the air layer d 0 , c 0 , d * 0 , c * 0 can be obtained from the four boundary conditions at the interface between the ionosphere and the air layer. Then, we can apply the boundary conditions and the left equations in A1 and A2 for calculating the unknown parameters. The EM field components at the Earth surface can be written as where

Layer-matrix method
For the layer-matrix method, the air layer and ionosphere in the upper half-space can have arbitrary layers, and the lithosphere layer in the lower half-space can also have arbitrary layers. The full-space Green's function in the wave number domain is derived by using the layer-matrix method, while the numerical solution of the Green's function in the frequency domain can be obtained by Fourier integration.
In the layer-matrix method, the air layer is set as the 0-th layer, the ionosphere is composed of M-layers and the lithosphere layer is composed of N-layers. The depth of the first layer in the upper halfspace is z − p , where p = 1, 2, . . . , M , the depth of the lower half-space is z p , where p = 1, 2, . . . , N. It can be deduced that ⎡ In the above equations, u x and V z have the same meaning as those in the R-function method X and V , and a R p and a −1 R p are the layer matrix and the inverse matrix for the p-th layer. Similar to the R-function method, after the coefficients of the air layer are obtained by using the boundary conditions, the values of each component of the surface field can be obtained. Based on the boundary conditions, the Fourier integral of the EM field in each layer can also be derived.