How Accurate Geomagnetic Regional Modeling Using Ordinary Kriging For Clustered Distributed Data?: New Insight From Indonesian Archipelago

13 Indonesia relies only on the limited number of repeat station networks due to the archipelago 14 setting with the extensive sea with the clustery distributed pattern. This paper explored 15 geostatistical modeling to overcome that typical data characteristic. The modeling used 16 repeat station data from the 1985 to 2015 epoch. The research used ordinary kriging (OK) 17 compared to the Spherical Cap Harmonic Analysis (SCHA) and Polynomial. The results 18 show that the root means square error (RMSE) of each declination, inclination, and total 19 intensity vary among epochs. OK method for declination component produces smaller 20 average RMSE (7.67 minutes) than SCHA (9.26 minutes) and Polynomial (7.97minutes). 21 For the inclination component, OK has an average RMSE of 9.55 minutes, smaller than 22 SCHA (10.05) but slightly higher than Polynomial (9.36 minutes). For the total intensity 23 component, OK produce an average RMSE of 63.58 nT, smaller than SCHA (82.24 nT) and 24 Polynomial (68.97 nT). The finding shows that the kriging method can be a promising 25 method to model the regional geomagnetic field, especially in the area of limited available 26 data and clustered distributed data.


31
The value and position of the geomagnetic field are dynamic. This is caused by its fluid and 32 material motion (Elsasser, 1956;Brandenburg, 2007). Therefore, some observation and 33 modeling are needed to produce a geomagnetic field map periodically. Geomagnetic field 34 mapping began in 1957 at the International Union of Geodesy and Geophysics (IUGG) 35 meeting in Toronto. The meeting became the birth of the World Magnetic Survey (Heppner,36 1963). In 1960, a mathematical model of spherical harmonic analysis was recommended and 37 the model became the basis for the International Geomagnetic Reference Field/IGRF 38 (Zmuda, 1969;Macmillan and Finlay, 2011). The IGRF is an international agreement 39 generated from the spherical harmonic model of the geomagnetic field and its secular 40 variations (Mandea and Macmillan, 2000). 41 The IGRF model is updated every five years under the auspices of the International 42 Association of Geomagnetism and Aeronomy (IAGA). Improvement to the IGRF's accuracy 43 continues to increase by adding the truncation order of the spherical harmonic analysis 44 equation incorporated by satellite data. However, it still has limitations in terms of accuracy 45 and spatial resolution (Maus et al., 2005). Furthermore, the spherical harmonic analysis 46 equation cannot be used for a geomagnetic field in the form of an orthogonal area (Mandea 47 and Purucker, 2005). This has led to the introduction of new regional geomagnetic field 48 models in the hope to produce a more accurate model. Indonesia repeat stations comply with some special terms determined professionally by the 56 IAGA (Newitt et al., 1996). Moreover, the nature of Indonesia as an archipelago causes the 57 clustered distribution of the repeat stations. The cluster distribution is shown by the nearest 58 neighbor analysis with negative value (Z score = -0.78). Repeat station distribution is 59 difficult to make uniformly. The repeat station locations cannot be placed on the water, the 60 site can't be accessed. No locations were found that met the standards for repeat stations 61 (Newitt et al., 1996). This concern gives more difficult challenges to the regional 62 geomagnetic field modeling framework in Indonesia (Fig. 1).

63
In general, the methods used in regional geomagnetic field modeling are the polynomial 64 modeling method (Alldredge, 1981;Düzgit et al., 1997), the spline and the Taylor expansion 65 (Hall and Meyer, 1976;Wahba, 1990), spherical cap harmonic analysis (Haines, 1985), 66 revised spherical cap harmonic analysis (Thébault et al., 2004), and revised spherical cap 67 harmonic analysis for the Earth's surface (Thébault, 2008). The existing method generally 68 needs dense data for modeling. The boundary effect of the existing method is also affected 69 by the data configuration. Numerical instability due to wide data gaps resulting from sparse 70 data patterns needs to be balanced with statistical or physical regulation to increase model 71 reliability (Torta, 2020).

72
To propose a new regional geomagnetic field modeling with limited and clustered data 73 patterns, the geostatistical method can be an alternative method. Geostatistical methods are 74 very accurate in modeling natural resources (Goovaerts, 1997), such as underground river 75 flows (Huysmans and Dassargues, 2009) and reservoirs of oil and gas (Poon et al., 1993;Xi 76 and Morgan, 2019). Geostatistical methods were developed to predict the distribution of the 77 percentage of ores in mining (Krige, 1951). A geostatistical method is also defined as a 78 method to predict the value of the data property, which is not covered by the sample data 79 distribution or sparse data pattern (Webster and Oliver, 2008).   However, the geostatistical method has never been used for regional geomagnetic field 93 modeling. Searching on the Scopus indexed papers using the keywords "geomagnetic" and 94 "regional," which obtained 3,524 papers. Among the obtained papers, there are only a few 95 which contain the 'geostatistic' keyword ( Fig. 2). There are only two papers that contain the 96 'geostatistic' word. One paper is discussing the history of regression and model-fitting in 97 Earth science (Howarth, 2001), and the second paper discusses the estimation of regional 98 covariance of Total Electron Content (Ghoddousi-fard, 2019). Neither of the papers 99 discusses modeling the regional geomagnetic field. Accordingly, the purpose of this paper 100 is to explore the accuracy of the geostatistical method in regional geomagnetic field 101 modeling and mapping.  104 We used repeat station data from the BMKG in the epoch 1985 to 2015. The BMKG has 105 conducted measurements in repeat stations to cover the Indonesia region starting from 1985.

106
The location for most of the repeat stations is located at the airport. It has the advantage that 107 apart from fulfilling the special requirements for the location of the repetition station, it can 108 also support the calibration of runway azimuth, which is important for aircraft navigation   The Polynomial method is selected as the sample for the interpolation method group because 132 it is more accurate than the SCHA, as is evident in the regional geomagnetic modeling  and equation (2) is used for the polynomial method with θ0 = 2° and ϕ0 = 117.5°. We then 164 analyze the resulted estimation data to determine the mean root square error (MRSE) level. 165 Therefore, we can investigate the accuracy of the geostatistical method in modeling the 166 regional geomagnetic field compared with the existing method.

169
Using the sample data as shown in Table 1

187
The smallest average RMSE value of 6.68 minutes of the inclination component is obtained 188 using the OK method with the stable variogram (Fig. 3b). The OK method, using the 189 exponential variogram and the polynomial methods, produces slightly smaller differences of The RMSE value fluctuation of the Total Intensity Component is the same as that of the 195 Inclination Component (Fig. 3c). This value confirms the contours of the two in the sample

202
The results show that the geostatistical method produces a good estimation or model in the  accurate in straight geomagnetic isoline (Ryan, 1997).

227
The result also shows that the number of repeat stations affects the RMSE, as suggested by Component. Generally, the variogram has intrinsic characteristics that affect to estimate of 242 data in each direction (Webster and Oliver, 2008). Variogram has 3 parameters such as 243 nugget, range, and sill (Armstrong, 1998). Nugget is a constant value for all distances greater

253
The OK method gives a good result in modeling the regional geomagnetic field. This 254 research gives the geostatistical method a chance to play its new role in modeling the 255 regional geomagnetic field in the area of limited and clustered distribution pattern data. Competing interest 285 We declare that we have no significant competing financial, professional or personal 286 interests that might have influenced the performance or presentation of the work described 287 in this manuscript.     Scopus-indexed papers using the keywords "geomagnetic" and "regional". The geostatistic 469 keyword has only two occurrences.