Regional geomagnetic core field and secular variation model over the Iberian Peninsula from 2014 to 2020 based on the R-SCHA technique

The Earth’s magnetic field originated in the fluid core, the so-called core field, is the dominant contribution to the geomagnetic field. Since ancient times, the core geomagnetic field has been used primarily for geographical orientation and navigation by means of compasses. Nowadays, thanks to the large amount of geomagnetic data available, core field models can be developed on a global or regional scale. Global models resolve large-scale geomagnetic field features, while regional models can resolve greater detail over a particular region. The spherical harmonic cap analysis is a widely used technique for regional-scale modelling of the geomagnetic field. In this work we have developed a regional model of the core field and its secular variation between 2014.5 and 2020.5 over the Iberian Peninsula, based on data from Swarm satellites, geomagnetic observatories and repeat stations. Its performance has been validated by comparing the fit to the available geomagnetic data using the regional model and the global models IGRF and CHAOS over the whole spatio-temporal range studied. In order to optimise the model, a detailed study of its input parameters has been carried out, showing that not all parameters have an equal influence on the modelling. This new model reproduces the input data with a root mean square error of 2.9 nT, improving the outcome of global models on this region. The results of this work will allow the Spanish Instituto Geográfico Nacional to produce the magnetic cartography of Iberia and the Balearic Islands in 2020.0, which for the first time will be based on a regional core field model, replacing the polynomial variation method used in the past.


Introduction
The first Spanish geomagnetic chart was made in the early decades of the twentieth century for the epoch of 1924.0 by the National Geographic Institute (Instituto Geogràfico Nacional, IGN), at that time Geographic and Cadastral Institute (Instituto Geográfico y Catastral, IGC) (IGC 1927).This cartography consisted of three charts with the isogonic (lines of equal magnetic declination), isoclinic (equal inclination) and horizontal isodynamic curves (equal horizontal component).After that, since 1975 and following the IAGA's (International Association of Geomagnetism and Aeronomy) recommendation, the Geomagnetism Service of the IGN publishes the geomagnetic chart of the Spanish mainland and the Balearic archipelago every 10 years, and the chart of magnetic declination every 5 years, all of them available at IGN website (IGN 2023).These charts are generally used as a reference in navigation and geographic orientation applications and for studies of the core, crustal, and external geomagnetic field characterisations.The last release was published in 2017, and contained the map of declinations, the horizontal and vertical components, and the field strength for the epoch 2015.0 (IGN 2017)).Since these last charts were published, the IGN has not provided a new set of geomagnetic maps.
Up to now, the geomagnetic charts provided by IGN have been developed using ground geomagnetic data.The ground data come from three different sources: 1. Mapping-station data: also called mapping-points, they are absolute observations (sum of the core and crustal geomagnetic fields) taken during the 80s in a dense and regular grid over the Iberian mainland and Balearic archipelago.They provide information about the crustal field.2. Repeat station data: absolute observations every 2-3 years in around 30 locations to provide information about the secular variation.3. Observatory data: data coming from Spanish observatories that provide a good constraint for the core field and its secular variation.
In order to develop these geomagnetic charts, a polynomial fitting was applied to the secular variation (SV, i.e. the first time derivative of the geomagnetic field components) of the core geomagnetic field provided by the repeat stations and observatory data.The temporal fitting translates the absolute geomagnetic observations up to a common epoch, while the spatial fitting provides the lines of equal secular variation, the so-called isoporic lines, at the same epoch.Figure 1 shows an example of this kind of chart, with the 2015-map of declinations obtained by this method (IGN 2017).
In this work, we take advantage of both new sources of geomagnetic data (i.e.satellite data) and novel mathematical approaches.In terms of new data, the satellite era has brought a revolution in geomagnetic studies.Very high-quality and high-resolution magnetic databases are now available, for example, from the European Space Agency's (ESA) Swarm constellation (Friis-Christensen et al. 2006) launched in 2013.These satellite data, together with the ground data detailed above, constitute a robust database through which it is possible to develop highly accurate geomagnetic charts for a target region and a considered epoch.Regarding the mathematical and physical approaches, there are techniques that appropriately fit the data to more realistic models than a simple polynomial fitting and satisfy the Laplace's equation.At global scale, the often used approach is based on spherical harmonic functions in space and cubic splines functions in time.Examples of global geomagnetic models are the CHAOS (Finlay et al. 2020), IGRF (Alken et al. 2021) or CM (Sabaka et al. 2020) families, among others.However, when the distribution of geomagnetic data is concentrated in a small region (as in this work), the optimal way to model the geomagnetic field is to apply regional modelling techniques, which allow a more accurate definition of the geomagnetic field in that region than any existing global model (Thébault et al. 2006;Nahayo et al. 2018;Nahayo and Korte 2022).In this context, the most appropriate technique when the geomagnetic data are located at different altitudes (ground and satellite data) is the Revised Spherical Cap Harmonic Analysis (R-SCHA, Thébault et al. 2006), in which the problem is solved in a spherical cone.This particular modelling technique has recently been applied in other regions of the world, e.g.Talarn et al. (2017), Vervelidou et al. (2018).Here, we provide the regional core field model that will serve as a basis for the geomagnetic cartography of the Iberian Peninsula and Balearic archipelago for the epoch 2020.0 by IGN, replacing the polynomial fitting method previously used by the IGN in earlier times.

Data
The present work used a database of magnetic data from three different sources: from space, via Swarm satellite data of the ESA, and from the ground, via geomagnetic observatories and repeat stations provided by the IGN's Geomagnetism Service.Details are given below.

Swarm satellite data
Provided by the ESA's Swarm constellation.The three satellites Alpha, Bravo and Charlie (hereafter SW-A, SW-B and SW-C) were put in orbit in November 2013 and are still in operation.These three identical satellites have a near-polar orbit, with inclinations of approximately 87 • .In the time interval spanned by this work, SW-A and SW-C flew at approximately 460 km altitude, in orbits separated by 1.4 • , while SW-B flew at approximately 510 km altitude perpendicular to the orbits of SW-A and SW-C.All three satellites are equipped with the same instrumentation, which consists, among others, of a scalar magnetometer (ASM) and a vector magnetometer (VFM) (Olsen et al. 2015); the measurements of the second one are those used in this work (i.e the X, Y and Z geomagnetic components).The data used are those classified by ESA as Level-1b (low resolution at 1 Hz).Before modelling, a battery of filters and corrections have been applied to properly select the satellite data: 1. Data on a spherical cap of 7 This data selection process was carried out in order to retain just the contribution due to the core field.After the application of these filters, the number of data was reduced by more than 80% (see Table S1 of the Supplementary Material for more information).Even so, the number of data points was too large for a computationally efficient modelling.Therefore, a decimation technique was used, which, in practice, consists of retaining the minimum information necessary to generate the regional model preserving a homogeneous spatial and temporal distribution of the data regarding the initial one.The procedure carried out for the decimation of the data consisted of: spatially, the data that were less than 0.5 • away from the points of a uniform grid of 1 • resolution were kept.For the temporal decimation, one datum per month was kept for each point of the grid.Indeed, as can be seen in Figure S1 of the Supplementary Material, the spatio-temporal distribution of the decimated data maintains the homogeneity of the original data.The spatial and temporal distribution of the Swarm satellite dataset used for the model are shown in Fig. 2a and b.Note that the Swarm satellite data are provided in a geocentric coordinate frame.

Geomagnetic observatory data
This ground dataset is composed by 1-min X (North), Y (East), Z (Vertical) data coming from all the geomagnetic observatories situated in the Iberian Peninsula (see locations in Fig. 2c): (1) San Pablo de los Montes Geophysical Observatory (SPT), operated by IGN, it is the reference observatory for the geomagnetic cartography published by IGN and also for the geomagnetic model we have developed.
(3) Ebro Observatory (EBR) operated by the Observatori de l'Ebre, an institution attached to the Ramon Llull University and governed by a foundation supported by Spanish institutions.As in the case of the Swarm satellite data, we have selected data for the interval 2014.5 − 2020.5 and between 00:00 and 2:00 local time (to mitigate the external field disturbances).We have selected only the measurements taken at magnetically quiet times, defined as days with the ap index less than or equal to 5 nT.This selection criterion for magnetically quiet observatory data is less stringent than for satellite data, as the satellite data are much more affected by the external field and the criteria of night-time and low magnetic activity provide sufficiently reliable results.Then, after transforming the observatory coordinates and the geomagnetic field components X, Y and Z from geodetic to geocentric, we have removed the external field from data through the CHAOS − 7.4 model.The crustal contribution was retained, in order to obtain the anomaly biases of the observatories through the modelling process.The temporal distribution of the selected observatory data is given in Fig. 2d.As the Swarm data, the observatory data were also decimated randomly taking 5000 data from each observatory.Geomagnetic observatories make measurements every minute, which makes the temporal distribution of available data dense and homogeneous (Fig. 2d).When the decimation is applied, the data were retained by means of a uniform random distribution.This ensures that the uniformity of the original data distribution is maintained.

Repeat and mapping station data
Repeat stations, also called secular stations, are permanently marked points where precise observations of the Earth's magnetic field are made for a period of a few hours every 2-3 years.They are especially useful to obtain information of the secular variation of the core field.In this work, we analysed the field survey data recorded at geomagnetic repeat stations from official networks located in Spain, Portugal, and France inside the spherical cap (see Fig. 2c) and operated by the IGN, IPMA (Instituto Português do Mar e da Atmosfera, Portugal) and a IPGP (Institut de Physique du Globe de Paris, France), respectively.
We focused on the stations with data strictly included in the temporal window 2015-2020 and with enough number of data for obtaining the SV (at least two in that period).Due to the fact that Portuguese data were outside this time interval, they were not included in our model.We worked with the data from 44 Spanish and 3 French repeat stations (Fig. 2c).These values consist of the D, H, Z, F (declination, horizontal component, vertical component, and intensity) and D, I, F (declination, inclination and intensity), respectively, which were transformed to X, Y, Z for consistency with the rest of data.As in the case of observatory data, we have rotated these data from a geodetic coordinate frame to a geocentric one.
Regarding the mapping stations data, we worked with 753 sets of D, H, Z, recorded at 751 Spanish mapping stations between June 1988 and October 1993, and 348 sets of D, H, Z located at 348 Portuguese mapping stations (grey dots in Fig. 2c).The values from the Portuguese stations correspond to the values translated to the year 1960.0,this is because they were obtained from Peña-Geromini (1963), who collected all the data from the mapping stations of Spain and Portugal that were used in the production of the Geomagnetic Map of the Iberian Peninsula for the epoch 1960.0,available at IGN website IGN (2023).

R-SCHA technique
Just as in other spherical harmonic techniques, the geomagnetic field vector is given by the gradient of a scalar potential that satisfies the Laplace's equation.At a regional scale, the best spatial geometry to solve the Laplace's equation is a truncated cone (see Figure 1 in Thébault et al. 2006) where the geomagnetic data are located inside.This approach is known as spherical cap harmonic analysis technique (SCHA) introduced by Haines (1985).This technique has been widely applied to model geophysical potential fields, i.e. in geomagnetism and gravimetry, at different spatial and temporal scales (see Torta 2020 for more information).For geomagnetic modelling purposes, the technique has been revised (R-SCHA, Thébault et al. 2006) providing a useful tool for modelling the geomagnetic data at different altitudes within the target region.From the fundamental equations of spherical trigonometry (see Kono 2015), the coordinates can be rotated to refer to the centre of the cap.In this new coordinate system, the R-SCHA scalar potential takes the following form (Thébault et al. 2006): where R is the mean radius of the Earth and r is the radial distance, P m n k (m) (cos�) are the Legendre's associated func- tions with real degree n k (m) and order m (k-m is even), cos(m ) and sin(m ) are the Fourier's series being and the colatitude and longitude referred to the centre of the cap.R p (r) and R 0 are the radial functions and K m p (cos�) are the Mehler's functions.Every term of Eq. 1 contains a set of time-dependent coefficients.Similarly to the global spherical harmonic expansion (SHA) coefficients, these are also called Gaussian coefficients.The set of coefficients {g m i,n k , h m i,n k } and {g m e,n k , h m e,n k } are the internal and external coefficients, respectively, while {g m p , h m p , g m 0 , h m 0 } are the Mehler's coefficients both depending on time t.All these coefficients characterise the expansion and may be solved by inverting the geomagnetic data.Unlike in the case of global SHA, the R-SCHA cannot separate the contribution of external origin from that of internal (1) origin.This is why the coefficients have more geometrical interpretation than physical meaning (Schott and Thébault 2011).Potential expansions are truncated up to the maximum spatial degrees K int , K ext , P max and M max .
The cap geometry is determined by the half-angle 0 and the lower and upper limits, a and b, of the truncated spherical cap (see Figure 1 in Thébault et al. 2006).
A geomagnetic field model is obtained from the inversion of magnetic data.The problem of non-uniqueness in geophysics implies that when inverting the data to obtain a model, there is no unique solution but infinite models that reproduce the input data (Backus et al. 1996).In addition, solutions are not stable, which means that large changes in the model parameters can be translated into small (or even no) changes when fitting the data.The problem's solution involves adding information beforehand to choose a single solution among the infinite proposed solutions.This information is contained in the spatial ˆ and temporal ˆ regularisation matrices, which are responsible for stabilising, smoothing and leading the models towards a physically meaningful solution.Similarly to the global modelling using SHA technique, we use analogous regularisation matrices, which minimise the magnetic field norm and the time derivative of the radial component of the field at the core-mantle boundary (CMB) (see Korte and Holme 2003 for further details).

The inversion of the geomagnetic data (d) to obtain the Gaussian coefficients (g) of the model is as follows:
where A is a matrix that depends on the radial distance, latitude and longitude, α and τ are the Lagrange multipli- ers set to provide the best fit trade-off between the data and the model; and and are the spatial and temporal regularisation matrices.In general, when constructing models of the actual field, not all geomagnetic data are of equal quality.It is necessary to reflect this in the inversion process by means of the weights matrix, W, which is the inverse of the variance matrix, i.e. the squared uncertainty of the data (Talarn et al. 2017).In addition, in this kind of approach an initial reference model is used in the modelling process.In this case, we use the dipole term (n = 1) given by CHAOS − 7.4 model.In this way, the neces- sary zero divergence boundary condition is fulfilled (see section 5 of Thébault et al. 2004).
In order to obtain the optimal R-SCHA parameter values, we have performed a stochastic computer simulation, often called Monte Carlo simulation (Rubinstein and Kroese 2016).In the Monte Carlo numerical-statistical method, 100,000 synthetic R-SCHA models were obtained, one for each random combination of parameters.The field values provided by each synthetic model were compared with respect to the input values, i.e. the CHAOS − 7.4 values.The parameter value range was determined based on a starting value of the desired cap size (i.e. 7 • radius and from surface to maximum Swarm height, 590 km) and R-SCHA parameters used in previous studies (e.g.Talarn et al. 2017).

Optimisation of the parameters
Prior to the modelling, a study was performed to optimise the SV trend function and the R-SCHA parametrisation.It was carried out on the basis of synthetic data, i.e. data obtained from a global field model (in our case, the CHAOS − 7.4 model) with a spatial and temporal dis- tribution identical to the real data available in the area and time window under study.It determined that linear or quadratic functions do not reproduce the SV with enough accuracy, especially the vertical component.Therefore, the SV had to be modelled by cubic B-splines (De Boor 2001).
The optimised R-SCHA parameters were: -The limits of the cap, given by the half-angle 0 and the radial limits a and b.With the aim of minimising edge effects and since 0 , together with the limit of the series expansions, determines the maximum harmonic degree which in turn influences the model minimum wavelength to be resolved.-The limits Kmax (inner and outer), Pmax and Mmax of the expansions of Eq. 1. -The number of knot points used in the spline interpolation (N).In a spline interpolation, the knot points represent the points where there is a change in the regression polynomial.-The smoothing parameters ( α and τ ).Optimal val- ues of the spatial and temporal regularisation.These parameters need to be selected to avoid oversmoothing or the creation of fictitious variations in the final model.
According to the previous section, in the Monte Carlo method the optimal parameter combination is the one that minimises the RMS (Root mean square) residuals of the synthetic field values (input data) with respect to the magnetic field values provided by the model (modelled data).It was established that the modelling error is minimised when: the limits of the expansions are set at 19 and 20 for internal and external Kmax, 9 for Pmax and 6 for Mmax; the cubic spline interpolation uses 3 knot points; the cap size is 11 • and radial boundaries are 6365 km and 6924 km.This combination provides, for a spacetime regularisation fixed with α and τ at 1 µT −2 and 1 µT −2 yr 4 , a residual RMS of 0.58 nT.
As spatio-temporal regularisation strongly influences the modelling results, it was analysed separately once the optimal R-SCHA parameters of Eq. 1 were obtained for a fixed spatio-temporal regularisation.We calculate the RMS residual of the synthetic models obtained with 25 parameter combinations of α and τ with logarithmi- cally spaced values between 0.1 and 1000 (in their respective units).We obtained that the values that minimise the RMS of the residuals are 0.1 µT −2 and 0.1 µT −2 yr 4 , respectively; with this regularisation, the synthetic RMS is reduced to 0.23 nT.From these results, we can see that the regularisation parameters have a large impact on the model results.In particular, spatial regularisation influences the results the most, possibly due to the limited size of the cap.Temporal regularisation has less influence on the RMS, possibly due to the fact that the core field changes smoothly over time.
Besides, it has been found that not all R-SCHA parameters have the same weight in the outcome of the model, i.e. there are values or ranges of values that give significantly better results.When the Monte Carlo results are plotted in a bivariate histogram (RMS residual against the parameter range, Fig. 3 and Figure S2 in the Supplementary Material), it can be easily seen that for the inner and outer Kmax limits, the number of knot points and the amplitude of the cap, the trend of the lower residuals is highly polarised towards a particular parameter value.The distribution is characterised for particular intervals of these parameters where a significantly higher number of combinations give the smallest residuals.
On the other hand, the radial size of the cap, a and b, and the expansion limits, Pmax and Mmax, present a more uniform distribution of the RMS.This means that they do not influence the results of the model with as much weight as the rest of the parameters.It is possible that they also have a preferential value, but are masked by the variation in the parameters with the greatest impact.Figure 3 shows an example of each situation, a bivariate histogram of the Monte Carlo results for one parameter (the knot points number) with high weight (left) and another (lower radial limit of the cap) with low weight (right).The remaining histograms can be found in Figure S2 in the Supplementary Material.

Results and discussion
Core field model for 2020.0 After performing the test with synthetic data to establish the optimal parameters, we applied the R-SCHA technique explained in the Methods section to the set of satellite and ground data (Fig. 2) to get the final regional model.However, because the real data are much noisier than those synthesised by a global model, the limits of the expansions had to be reduced and the angular size of the cap had to be enlarged.The values used for the realisation of the model can be found in Table S2 of the supplementary material.Therefore, the final model was obtained by inverting the data for 63 coefficients (equivalent to a global spherical harmonic degree n = 13) through the inversion shown in Eq. 2. Taking into account the spline knot point number, a total number of 441 coefficients were obtained, which constitute the regional model for the time period between 2014.5 and 2020.5 that reproduces both the core field and its secular variation at any point inside the spherical cap.The model shows a mean squared error with the input data of 2.9 nT and a mean value of the residuals of 0.0 nT (considering the three XYZ components together).Eventually, it was possible to obtain the geomagnetic core field and its secular variation maps for the area of interest for the epoch 2020.0 (Fig. 4).
One of the characteristics of the geomagnetic field that may be highlighted is the crossing of the agonic line (line of zero declination) through the Peninsula.In Fig. 4, it is indicated by a dashed line on the eastern component; in the year 2020 the agonic line crossed the Peninsula in a N-S direction from the Bay of Biscay.
The temporal evolution of the geomagnetic field intensity and its secular variation in this period is shown in Figures S6 and S7 of the Supplementary Material.From them we can highlight that the intensity has been increasing throughout the whole period, as it does in the Eurasian region, unlike the global intensity of the geomagnetic field, which has been decreasing in recent years (Alken et al. 2021).When looking at the secular variation of the eastward component of the R-SCHA model, the graph shows a -shape around 2019.0.This behaviour of the secular variation of the model corresponds to the 2020 geomagnetic jerk proposed by Pavón-Carrasco et al. (2021), which also identifies it as -shaped in the Iberian Peninsula region.During 2017, period also covered by our model, a geomagnetic jerk has also been detected in other regions of the planet (Hammer et al. 2021).If we focus on the behaviour of the SV in 2017, no change seems to be observed in the trend of the eastern component.This reaffirms the regional character of the 2017 jerk, which is not observed equally across the planet, in contrast to the global character of the 2020 jerk (Pavón-Carrasco et al. 2021).
Total field for 2020.0 At the Earth's surface, the measured internal magnetic field is the result of the sum of two contributions: the core field, the one generated in the Earth's outer core, and the crustal field, i.e. the remanent magnetic field recorded by the rocks located in the upper lithosphere above the Curie isotherm.Thanks to the mapping-station data (see Methods Section) and with the R-SCHA model for the core field, it has been possible to represent the total internal field maps at 2020.0 epoch (Fig. 5).The procedure carried out is described in the next statements.Firstly, each mapping-station point was translated from its original location and date to the year 2015.0 using the secular variation given by the COV-OBS.x2global geomagnetic field model (Huder et al. 2020).Thereafter, the data referring to the DHZ elements have been transformed into XYZ.The R-SCHA core field in 2015.0 has been subtracted, obtaining the crustal field for the Iberian Peninsula and Balearic Islands.As the crustal field is time-invariant, it was added to our 2020.0core field, thus obtaining the total field maps (Fig. 5).It can be observed that the field isolines are not parallel as in the case of the main field (Fig. 4), where the spatial variations are of a much larger scale than the region studied.This deformation in the isolines is greatest where the crustal magnetic field is strongest (see Figure S8 of the Supplementary Material).This corresponds to the Central Iberian and Ossa-Morena geological zones, as well as other local anomalies such as the Basque magnetic anomaly (North of Spain).The analysis of their geological origin is beyond the scope of this paper, but several specialised studies have been carried out by other authors (e.g.Ayarza et al. 2021;García Lobón et al. 2017;Calvin et al. 2014).

Comparison with global models
Global models of the geomagnetic field can also be used to obtain the magnetic field over the Iberian Peninsula.For instance, the global models IGRF and CHAOS, which provide a maximum harmonic degree of 13 and 20 for the core field, can resolve wavelengths of 3000 km and 2000 km, respectively.In this section, we compare our regional model with these global ones to see the advantages of using regional models in areas with enough data in contrast to global ones.
The spatial differences obtained between our regional model and the global ones (the IGRF-13 up to harmonic degree 13 and the CHAOS − 7.4 up to harmonic degree 15) in both the geomagnetic components and their secular variation at surface height in 2020.0 are plotted in Figs. 6 and 7.The root mean square of the mean difference between model R-SCHA and CHAOS of X, Y  S3 of the Supplementary Material.We have also compared how well the different models reproduce the magnetic dataset.A table with the mean values and standard deviation of the data-model differences (Table S4) as well as histograms of the distribution of the residuals (Figures S3, S4 and S5) are provided in the Supplementary Material.The result is that the R-SCHA model slightly reduces the standard deviation (2.9 nT) and the mean value (0.1 nT) of the model-dataset residual (for IGRF were 4.0 ± 8.0 nT, while for CHAOS 0.0 ± 3.4 nT).
In addition, the time evolution of the global models has been compared with the data from the geomagnetic observatory of San Pablo de los Montes, located exactly in the centre of the spherical cap of validity of the model.As can be seen in Figure S6 of the Supplementary material, the R-SCHA model fits the time evolution of the field measured at the observatory lightly more accurately than the global models, especially the IGRF, which fits the secular variation with a straight line (due to the linear interpolation between 2015 and 2020).The CHAOS model is much more similar to our R-SCHA, as it also uses a spline interpolation (but a sixth-order base of splines one), to model the secular variation (Finlay et al. 2020).The CHAOS model reproduces the same trend as the R-SCHA but with a little larger standard deviation.Moreover, the regional and CHAOS models present similar behaviour when the SV is analysed by means of the San Pablo observatory data (see Figure S7).However, some differences are observed in the North component and in the rest of the components from 2019 onwards, when the 2019-2020 geomagnetic jerk occurred (Pavón-Carrasco et al. 2021).In light of these results, the regional R-SCHA model improves the reconstruction of the field values and their secular variation for the given study area and time with respect to the global models.

Conclusions
In this work, we have developed a model of the core field and its secular variation of the Iberian Peninsula and Balearic Islands between 2014 and 2020.We have used a spherical harmonic analysis in a spherical cap, the socalled R-SCHA technique.The model was developed by inverting the geomagnetic data measured by the Swarm satellites, four geomagnetic observatories and repeat stations in the target area.The obtained model has a mean residual with the input data of 0.1 nT and a standard deviation of 2.9 nT, improving the results of the global From the optimisation study with the Monte Carlo method based on synthetic data, the parametrisation of the R-SCHA and the spatio-temporal regularisation that minimises the residual has been established.The parametrisation obtained as optimal is as follows: the limits of the expansions set at 19 and 20 for internal and external Kmax, 9 for Pmax and 6 for Mmax; 3 knot points for the spline interpolation; the cap size is 11 • and radial boundaries are 6365 km and 6924 km.
From the distribution of the residuals for the different random parametrisations of the Monte Carlo method, it has been found that not all parameters influence the results of the model equally.There are some parameters that have more weight, forcing most of the residuals closer to zero which correspond to parametrisations with a certain value or range of values of these parameters.These are internal and external: Kmax limits, the number of spline knot points and the amplitude of the cap.On the other hand, the remaining parameters (radial size of the cap, a and b, and the expansion limits, Pmax and Mmax) provided, in the range studied, a uniform distribution of the residuals.This indicates that their weight in the model is minor, and that their influence is masked by the variability of the rest of the parameters in the random Monte Carlo process.
In summary, for the first time a regional geomagnetic field model satisfying the Laplace's equation has been developed for the Iberian Peninsula and Balearic Islands, using not only ground but also satellite data to obtain the total internal field maps that will constitute the geomagnetic charts of Iberia and the Balearic Islands by 2020.0 provided by IGN.

Fig. 2
Fig. 2 Spatio-temporal distribution of the final magnetic data.a Decimated Swarm spatial and b temporal distribution from 2014.5 − 2020.5.The width of the histogram stacked-bars represents one month.c Spatial distribution of the ground data, observatories (triangles), repeat stations (green circles) and mapping stations (grey dots).d Ground data temporal distribution; repeat stations (above), the width of the histogram bars represents 1 year, and observatories (below) with the width of the histogram stacked-bars representing 1 month

Fig. 3
Fig. 3 Bivariate histograms of the RMS of the residual of the range of parameters, limited to RMS<10.The purple stars indicate the value obtained as optimal in the Monte Carlo method.a Number of nodes (N) and b Lower limit of the cap (a)

Fig. 4
Fig. 4 R-SCHA core magnetic field components: a X, b Y, c Z and d intensity (colormaps) and its secular variation (blue lines) in nT/yr in 2020.0 calculated at surface height.The dashed white line of b indicates the null east component (Y) line, i.e. the agonic line

Fig. 5
Fig. 5 Total geomagnetic field components in 2020.0.a X b Y c Z d F in nT calculated at surface height

Fig. 6 Fig. 7
Fig.6Difference of the geomagnetic (from top to bottom: X, Y, Z and F) given by our regional model and the global models CHAOS − 7.4 (left column) up to degree 15, and IGRF-13 (right column) up to degree 13 at surface height in 2020.0