Aeromagnetic survey in Kusatsu-Shirane volcano, central Japan, by using an unmanned helicopter

Kusatsu-Shirane volcano is one of the active volcanoes in Japan. Phreatic explosions occurred in Mt. Shirane in 1983 and most recently, in 2018, in Mt. Motoshirane. Information on the subsurface structure is crucial for understanding the activity of volcanoes with well-developed hydrothermal systems where phreatic eruptions occur. Here, we report aeromagnetic surveys conducted at low altitudes using an unmanned helicopter. The survey aimed to obtain magnetic data at a high spatial resolution to map the magnetic anomaly and infer the magnetization intensity distribution in the region immediately after the 2018 Mt. Motoshirane eruption. The helicopter used in the survey was YAMAHA FAZER R G2, an autonomously driven model which can fly along a precisely programmed course. The flight height above the ground and a measurement line spacing were set to ~ 150 m and ~ 100 m, respectively, and the total flight distance was 191 km. The measured geomagnetic total intensity was found to vary by ~ 1000 nT peak-to-peak. The estimated magnetization intensity derived from measured data showed a 100 m thick magnetized surface layer with normal polarity, composed of volcanic deposits of recent activities. Underneath, a reverse-polarity magnetization was found, probably corresponding to the Takai lava flow in the Early Quaternary period (~ 1 Ma) mapped in the region. Our results demonstrate the cost-effectiveness and accuracy of using drone magnetometers for mapping the rugged terrain of volcanoes.


Introduction
Kusatsu-Shirane volcano is one of the most active volcanoes located in the central part of Japan. In the area surrounding Kusatsu-Shirane volcano, some Quaternary volcanoes have formed over the Early Quaternary or Tertiary andesitic basement rocks (e.g., Hayakawa and Yui 1989;Terada 2018). At the first stage of the formation history of Kusatsu-Shirane volcano, Matsuozawa lava flow erupted at ~ 0.5 Ma, followed by the Horaguchi lava and Oshi pyroclastic flows (Kaneko et al. 1991).
Following another eruptive stage, the three main pyroclastic cones, Mt. Shirane, Mt. Ainomine, and Mt. Motoshirane, formed approximately 10,000 years ago and are covered by volcanic deposits with a thickness of a few 100 m (Hayakawa and Yui 1989). Recent activities were almost only observed in Mt. Shirane, and frequent phreatic explosions occurred every couple of 10 years.
Mt. Ainomine and Mt. Motoshirane have been inactive for 1500 years (Nigorikawa et al. 2016). At the east flank of Mt. Motoshirane, however, there is a large amount of strong acidic hot spring water discharged (as much as 10 ton per minute; Kikawada et al. 2008), implying that a notable widespread hydrothermal system exists beneath Kusatsu-Shirane volcano (e.g., Ohba et al. 2000).

Open Access
*Correspondence: tkoyama@eri.u-tokyo.ac.jp 1 Earthquake Research Institute, The University of Tokyo, Yayoi 1-1-1, Bunkyo-ku, Tokyo 113-0032, Japan Full list of author information is available at the end of the article In active volcanoes with an aquifer, phreatic explosions occur because of the interaction of hot volcanic volatiles with groundwater. Mt. Shirane has experienced such eruptions. Phreatic explosions are usually not preceded by apparent precursory events; thus, it is difficult to predict their occurrence, which consequently leads to damages and casualties. The detection of tiny geophysical and geochemical signals preceding phreatic explosions as well as geological knowledge of the volcano is the key to accurate and timely predictions (Terada 2018).
The current grounded survey was carried out mainly at Mt. Shirane as it was recently active. Takahashi and Fujii (2014) revealed, through repeated ground-based geomagnetic surveys, the long-term thermomagnetic variations corresponding to volcanic activities caused by magnetizing or demagnetizing sources beneath the crater of Mt. Shirane. Tokyo Institute of Technology installed permanent stations of the magnetometer. In 2014, microseismic swarms and geomagnetic changes were detected at the same time in Mt. Shirane. The stations nearest to the Yugama crater showed a rapid decrease in geomagnetic total intensity. This implies that thermal demagnetization occurred just east of the Yugama crater of Mt. Shirane (Hashimoto et al. 2019). These geomagnetic changes may indicate that volcanic activity is getting gradually activated. In contrast, due to a lack of activity, Mt. Motoshirane was considered less likely to erupt than Mt. Shirane and hence, neither geomagnetic nor seismic stations were deployed in the region.
On January 23, 2018, however, there was a phreatic explosion at the Kagamiike-kita pyroclastic cone of Mt. Motoshirane (Ogawa et al. 2018). Volcanic activity on Mt. Motoshirane was low until then and the eruption was unexpected. The estimated mass of volcanic products in this event was not much, being less than 40 thousand tons (Ogawa et al. 2018;Himematsu et al. 2020), but some casualties happened because of lack of preparedness. Matsunaga et al. (2020) revealed the electrical conductivity structure of Mt. Motoshirane by magnetotelluric surveys. They detected a surface clay-rich alteration layer, ~ 0.5 km thick, and a deep and wide conductor which is interpreted to be a hydrothermal fluid reservoir. This reservoir can be a root of volcanic earthquakes occurring beneath Mt. Motoshirane and can also provide volcanic fluid mixed with meteoric water to the hot springs on the eastern flank of Mt. Motoshirane. A resistive body was also detected, which corresponds with basement rocks due to the geological results (Hayakawa and Yui 1989). Tseng et al. (2020) also revealed the electrical conductivity beneath the crater lake, Yugama, of Mt. Shirane, and detected a pathway between Yugama and a deep reservoir. These previous results imply that Kusatsu-Shirane volcano has a variable subsurface structure with some craters, and it is very important to understand it.
We conducted aeromagnetic surveys, with high spatial resolution, to elucidate the Kusatsu-Shirane volcano. Aeromagnetic surveys are used to investigate the inside of the volcanoes to detect thermal anomalies due to demagnetization of subsurface rocks and predict potential eruptive events (e.g., Nakatsuka 1994;Okuma et al. 1994;Finn et al. 2001;Okubo et al. 2005). We used an unmanned helicopter instead of a manually operated one to prevent safety hazards (Kaneko et al. 2011;Koyama et al. 2013;Hashimoto et al. 2014;Ohminato et al. 2017).
In this paper, we have presented an aeromagnetic survey conducted using an unmanned helicopter and the data processing procedure, followed by the estimation of the magnetization intensity distribution to discuss the subsurface structure of Kusatsu-Shirane volcano.

Aeromagnetic survey in Kusatsu-Shirane volcano
We conducted aeromagnetic surveys of Kusatsu-Shirane volcano, central Japan, in 2018, immediately after the Mt. Motoshirane eruption in January 2018, by using an unmanned helicopter to acquire magnetic field data of high spatial resolution, mitigating potential safety concerns. The unmanned helicopter, YAMAHA FAZER R G2, is driven by an internal combustion engine, and can fly continuously for over 2 h. A common drawback of the usage of internal combustion engines is the lack of combustion efficiency due to decreasing oxygen concentration at high altitude, and the previous model, YAMAHA RMAX G1, did not have the capacity to fly beyond Kusatsu-Shirane volcano the altitude of which is 2171 m (Morimoto 2017). However, FAZER R G2 is a state-of-the-art vehicle with a four-stroke engine (Kido et al. 2016), and can be operated even at altitudes of 2800 m (Morimoto 2017). The helicopter provides advantages, such as safety, and also has the ability to fly along a preprogrammed course, owing to which it can fly along the same course repeatedly and can be very useful for repeated surveys to detect temporal changes in geomagnetic fields. The time and position of the helicopter were recorded accurately up to a few tens of centimeters in the horizontal direction and < 1 m in the vertical direction in a quiet wind condition.
A cesium optical pumping magnetometer (Geometrics G-858 MagMapper Magnetometer) was used to measure the scalar value of geomagnetic total intensity up to 0.1 nT accuracy (Geometrics 2001). Electric power to the magnetometer was supplied by the helicopter generator. Time calibration was performed by the 1 Hz timestamping tool of a small portable GNSS device, Garmin eTrex H (Garmin 2007) via an RS-232C interface to the magnetometer. Note that positioning data acquired by the portable device are generally not very accurate (more than a few meters accuracy) and the positioning data acquired by the helicopter were used instead in the data processing of this study. The time of measurement was adjusted after the survey flights, according to the time of the flight record, containing the precise position information. The console of the magnetometer was fixed to the helicopter body. A magnetic sensor was hung down with a 5.5 m sensor cable to be installed away from the helicopter body in order to reduce the effect of the vehicle magnetization field. Thus, the position of the magnetic sensor can be assumed to be 5.5 m beneath the recorded position of the helicopter. The expected total error of measurement is ~ 10 nT including the effects of a vehicle magnetization field and positioning error due to the swinging of the sensor.
The total weight of magnetometer was ~ 4 kg, which was much less than the maximum payload of the helicopter (35 kg; Morimoto 2017). Measurements were taken at 10 Hz sampling. The flight speed was ~ 10 m/s, and the typical altitude above the ground was 150 m. A typical survey area measured ~ 3.5 km × 3.5 km. The survey line spacing was ~ 100 m, and the total flight distance was ~ 191 km. Figure 1 shows a plane map of the flight courses, colored by the measurement value of geomagnetic total intensity. The measured geomagnetic total intensity anomaly was as high as 1000 nT. A reference station of geomagnetic total intensity was deployed 4 km away from the survey area to eliminate daily variation in geomagnetic total intensity due to the external geomagnetic field. The same model of the magnetometer Geometrics G-858 MagMapper Magnetometer was used. 10 Hz sampling was done during flights, which is the same rate as measured by a magnetometer on the helicopter, and 1 Hz sampling was done in the nighttime to determine a midnight mean value as a baseline.

Data processing and modeling
The observed geomagnetic total intensity data in the helicopter were decimated to 1 Hz sampling for precise position and time calibration by the GNSS receiver, and outliers, such as spike noises, were removed. Thereafter, the magnetic field of both the Earth's core and external origins were subtracted from these data to extract the geomagnetic field component only due to volcanic rock origin. The geomagnetic main field of the core origin was assumed to be expressed by the CHAOS-7 geomagnetic field model (Finlay et al. 2020). This main field model has a better time resolution than a common IGRF model published every 5 years, and it can provide the main field data more precisely. The time variation of the external origin was assumed to be identical to that at the reference station. The variation was estimated using the reference field data by subtracting the midnight mean value. Finally, the extracted geomagnetic field data were decimated to provide evenly distributed spatial data, typically with 100 m spacing.
Using the calculated geomagnetic total intensity anomaly of volcanic origin, the magnetization intensity of Kusatsu-Shirane volcano was estimated, as the orientation of magnetization was assumed to be identical to the direction of the current geomagnetic main field in any position of our targeted area: inclination and declination are 51.0° downward and 7.9° westward, respectively.
The details of modeling are described in Appendix. Regarding the target volume V uni , a much wider area than the survey area, 26 km × 26 km, was considered to avoid artifacts of the edge effect of the model space. Surface topography was referred to as a digital elevation model with a 10 m horizontal square mesh published by GSI Japan. According to the locations of hypocenters, seismicity can be detected up to 1.5-2.0 km depth from the surface (e.g., Matsunaga et al. 2020), and brittle-ductile transition may occur at the bottom of the seismic region, at a depth of 1.4 km (Tseng et al. 2020). Thermal demagnetization of the andesitic rock of Mt. Shirane at 400 °C may also occur at this depth (Yamazaki et al. 1992); thus, the magnetizing thickness was assumed to be 1500 m from the surface in this study. A target volume V uni of 26 km × 26 km × 1.5 km was assumed to be magnetized, and other areas were assumed not to be magnetized at all. The uniform magnetization intensity m uni was estimated as 0.97 A/m. This is close to the results of laboratory measurements of the andesitic rock sampled in Mt. Shirane (Yamazaki et al. 1992). Standard deviation of data misfit between observed data and synthetic data by this uniformly magnetized model was estimated as 134.4 nT, which is too large compared to the expected data error (~ 10 nT).
Regarding the three-dimensional inverse modeling, the target area was assumed to be 19 km × 19 km, still much wider than the survey area, and the magnetizing thickness was again assumed to be 1500 m. A target volume, V ano , of 19 km × 19 km × 1.5 km, was discretized by blocks with dimensions 250 m × 250 m in horizontal spacing, and was divided vertically into four layers. The thicknesses of the four layers are 100, 200, 400, and 800 m from the top layer to the bottom, respectively. Each block was assumed to be uniformly magnetized. Note that each block is not an exact rectangular block, but has fluctuations at the top and bottom due to the topography with a horizontal resolution of 10 m. Figure 2 shows data misfit between observed geomagnetic total intensity and synthetic data using an optimal model. Standard deviation of data misfit was estimated as 12.7 nT, which is close to an expected data error. Figure 3 shows the magnetization intensity distribution of each layer of the optimal model. Note that it shows only an area of our survey and surroundings, and does not show all the inverted model space. The model is well inverted in the survey area over which the helicopter took flights. However, the model was not recovered in the area apart from the survey area due to the lack of sensitivity and model resolution. Additional file 1: Figure S1 shows the deviation of model parameters, which is derived from the ABIC minimization shown in Appendix. The first layer is well resolved, and a deviation of model is less than 3.5 A/m; the second layer is deviated by 1.5-5.8 A/m. The third is deviated by 6 A/m or more. The fourth layer has no resolution, and the deviation is about 20 A/m. Thus, hereafter, we discuss the magnetization in the first to the third layers.

Discussion
The top layer in Fig. 3a indicates that the entire surface is magnetized with positive polarity. A surface layer with a thickness of a few 100 m is composed of volcanic  Kaneko et al. 1991;Uto et al. 2004). Yamazaki et al. (1992) measured the remanent magnetization of pyroxene andesite, sampled close to Yugama Lake, and it exhibited the magnetization intensity as high as 1-2 A/m, which is consistent with a model in the top layer with weak anomalies. The second and third layers, shown in Fig. 3b, c down to a depth of 700 m from the surface, however, indicate negative or almost zero magnetization over a wide area. Akimoto et al. (2002) found that an orientation of the rocks in Yatsugatake, 75 km south of Kusatsu-Shirane, in about 1 Ma is almost opposite to the orientation of a current main field. Thus, an area with negative magnetization can be interpreted as a reversely magnetized region. Although this might have been caused by overweighting on this layer in the inversion calculation, rocks at these depths may be constituted of basement rock in the early Quaternary, Takai lava. This lava is supposed to be approximately the same age as Omeshi volcano, west of Kusatsu-Shirane volcano (Kaneko et al. 1991). K-Ar dating by Kaneko et al. (1991) shows that the age of Omeshi volcano is 1.10 ± 0.09 Ma. This means that the Takai lava erupted in the early Quaternary which was in Matuyama reversed Chron. Nagai et al. (2015) investigated core data from a drill hole, 3 km west of Mt. Motoshirane, and found a tuff dating to 1.29 ± 0.06 Ma at 171.5 m depth or deeper. Hayakawa and Yui (1989) indicated that there is a long inactive stage, spanning a couple of hundred thousand years between two eruptive stages, and the later eruptive stage may not have used magma, unerupted and solidified in the former eruptive stage. It may have cooled and magnetized lava stayed at this depth after the eruptive stage. Also, it might be highly altered and demagnetized, corresponding to a high conductive region of Takai lava (Matsunaga et al. 2020). Figure 4 shows the North-South vertical profile of the magnetization model through Yugama and Mt. Motoshirane, compared with the electrical resistivity model of Matsunaga et al. (2020). The second and third layers of magnetization model correspond to a conductive layer, which can be altered to clay. Therefore, it is consistent that the second and third layers appear less magnetized and partly reversely magnetized.
The detection of temporal changes in the geomagnetic total intensity by aeromagnetic measurements is very challenging. Figure 5 shows repeatability of data acquisition, that is, deviation of measured data for pairs of sites, which are closer than 5 (in gray color) and 10 (in white color) m to each other. A typical data error or standard deviation was as small as 10 nT or less in this survey. On the contrary, a temporal change detected by a groundbased survey was a few nT in the latest measurements (Hashimoto et al. 2019). It seems difficult to detect and is out of the scope of this research. Aeromagnetic survey, however, can measure a wide spatial pattern with high spatial resolution at precisely the same points, which means that spatial stacking to remove the noise can be performed. It may be possible to detect a minor geomagnetic change of volcanic origin with a wide spatial trend

Conclusion
In this study, highly spatially resolved magnetic field data were acquired from Kusatsu-Shirane volcano, and a geomagnetic field anomaly of approximately 1000 nT was detected using an autonomously driven unmanned helicopter. A three-dimensional magnetization intensity distribution was elucidated by using the measured geomagnetic field data. It shows thin normally magnetized lava at the surface underlain by altered and reversely magnetized lavas. This implies that the main basement lava, Takai lava, can be supposed to have erupted at the time in Matuyama reversed Chron, and have altered later. Aeromagnetic surveys have shown the potential to elucidate the wide and deep distribution of magnetization. Recently, some studies have reported aeromagnetic surveys using a common drone, which can be easily performed at a low cost (e.g., Utsugi and Hashimoto 2016;Malehmir et al. 2017;Koyama et al. 2019). It can be significantly improved in the near future and can be utilized in volcanic surveys for the detection of temporal changes in the geomagnetic field due to volcanic activity by repeated surveys, combined with ground-based surveys.

Appendix: Magnetization intensity inversion
In the first step, the average value of the magnetization intensity, that is, uniform magnetization intensity was estimated in the least squares sense for fitting all the data by removing the trend surface of the geomagnetic field anomaly, minimizing the following objective function, uni : where d n is the nth data parameter, which is the measured geomagnetic total intensity subtracted by the core field by the CHAOS-7 model and the external field variation, x n , y n , z n is a vector of the position of the geomagnetic sensor of the helicopter at the nth data parameter, and N is the total number of data parameters. The volume to be uniformly magnetized is given by V uni , and K x n , y n , z n ; x ′ , y ′ , z ′ is a kernel that expressed the geomagnetic total intensity anomaly at position ( x n , y n , z n ) generated by a unit magnetizing source at position ( x ′ , y ′ , z ′ ). The target uniform magnetization intensity is given by m uni and a 0 + a x x + a y y + a z z is an assumed trend surface to be removed, where a 0 , a x , a y and a z are constant parameters and estimated as well as m uni . The volume integral in Eq.
(1) was derived by following a common formulation for a rectangular prism with magnetization (Bhattacharyya 1964;Kunaratnam 1981).
(1) uni = N n=1 d n x n , y n , z n − m uni V uni K x n , y n , z n ; x ′ , y ′ , z ′ dx ′ dy ′ dz ′ − a 0 + a x x n + a y y n + a z z n 2 , The next step is the estimation of three-dimensional anomaly of the magnetization intensity deviating from the average intensity. Various methods of three-dimensional inversion of magnetization have been proposed (e.g., Li and Oldenburg 1996;Nakatsuka and Okuma 2014;Utsugi 2019), and the two-dimensional inversion method by Koyama et al. (2013) was adapted to the three-dimensional case in this study, which is a commonly used least squares method for fitting the data with model constraints for regularization to minimize the following objective function, ano .
(2) where d ′ n is the geomagnetic total intensity anomaly of the nth data, subtracting components generated by a uniform structure and a trend surface. V ano is the target volume for three-dimensional inversion. The magnetization intensity anomaly from a uniform model and the weight function on it are represented by m x ′ , y ′ , z ′ and w x ′ , y ′ , z ′ , respectively.
The objective function, ano , in Eq.
(2) can be rewritten using a discretized model anomaly parameter, m m , as follows: where M is the total number of model parameters.
KdV m x n , y n , z n is a volume-integrated kernel for the mth model parameter, m m . w m is the weight of the mth model parameter. is a hyperparameter that balances the data misfit term and model constraint (or roughness) term, which are the first and second terms of the righthand side of Eq. (3).
The model constraint or model roughness term is supposed to minimize the L2 norm of the weighted magnetization intensity anomaly model deviating from an average uniform magnetization intensity, m uni . The weight of the model anomaly w m is assumed to be the square root of the kernel of the model block to the position just above the block at a typical flight height above the ground.
where (x m , y m ) is the horizontal center position of the mth model block. Note that it is not identical to the weight function of Li and Oldenburg (2000), but our method avoids overweighting on blocks locating far from the survey area.
The hyperparameter, , to weigh the model constraint was determined by the ABIC minimization concept (e.g., Akaike 1980; Honsho et al. 2012;Koyama et al. 2013). An optimal hyperparameter could be uniquely determined to minimize the ABIC, and the ABIC minimization concept could select an optimal solution properly, even according to the L-curve criterion by Hansen (1992) (Additional file 2: Figure S2).  where C M i,j is an element (i, j) of the matrix C M . δ ij is the Kronecker delta.
Thus, a square root of each diagonal element, C M m,m 1/2 , can be supposed to simply represent a deviation of the mth model parameter, m m .