Global and Regional High-Resolution VTEC Modelling Using a Two-Step B-Spline Approach

: The ionosphere is one of the largest error sources in GNSS (Global Navigation Satellite Systems) applications and can cause up to several meters of error in positioning. Especially for single-frequency users, who cannot correct the ionospheric delay, the information about the state of the ionosphere is mandatory. Dual- and multi-frequency GNSS users, on the other hand, can correct the ionospheric effect on their observations by linear combination. However, real-time applications such as autonomous driving or precision farming, require external high accuracy corrections for fast convergence. Mostly, this external information is given in terms of grids or coefﬁcients of the vertical total electron content (VTEC). Globally distributed GNSS stations of different networks, such as the network of the International GNSS Services (IGS), provide a large number of multi-frequency observations which can be used to determine the state of the ionosphere. These data are used to generate Global Ionosphere Maps (GIM). Due to the inhomogeneous global distribution of GNSS real-time stations and especially due to the large data gaps over oceanic areas, the global VTEC models are usually limited in their spatial and spectral resolution. Most of the GIMs are mathematically based on globally deﬁned radial basis functions, i.e., spherical harmonics (SH), with a maximum degree of 15 and provided with a spatial resolution of 2.5 ◦ × 5 ◦ in latitude and longitude, respectively. Regional GNSS networks, however, offer dense clusters of observations, which can be used to generate regional VTEC solutions with a higher spectral resolution. In this study, we introduce a two-step model (TSM) comprising a global model as the ﬁrst step and a regional model as the second step. We apply polynomial and trigonometric B-spline functions to represent the global VTEC. Polynomial B-splines are used for modelling the ﬁner structures of VTEC within selected regions, i.e., the densiﬁcation areas. The TSM provides both, a global and a regional VTEC map at the same time. In order to study the performance, we apply the developed approach to hourly data of the global IGS network as well as the EUREF network of the European region for St. Patrick storm in March 2015. For the assessment of the generated maps, we use the dSTEC analysis and compare both maps with different global and regional products from the IGS Ionosphere Associated Analysis Centers, e.g., the global product from CODE (Berne, Switzerland) and from UPC (Barcelona, Spain), as well as the regional maps from ROB (Brussels, Belgium). The assessment shows a signiﬁcant improvement of the regional VTEC representation in the form of the generated TSM maps. Among all other products used for comparison, the developed regional one is of the highest accuracy within the selected time span. Since the numerical tests are performed using hourly data with a latency of one to two hours, the presented procedure is seen as an intermediate step for the generation of high precision regional real-time corrections for modern applications.


Introduction
In the last decades, the use of GNSS for a broad range of commercial and scientific fields has significantly increased. Especially GNSS applications which rely on high precision positioning and navigation became more and more popular. With higher demands on accuracy, the requirements for detecting and correcting errors in the measurements increase. The main error sources in precise GNSS are the effects caused by the atmosphere. Whilst the tropospheric effects can yield approximately 1m error in positioning (in along range, i.e., the height), the ionosphere effects can vary between 1 m and even 100 m. The ionosphere is defined as the part of the upper atmosphere where the number of charged particles (electrons and ions) is high enough to induce significant changes in the propagation of radio waves. The effect mainly depends on the signal frequency; signals with frequencies lower then 30 MHz can be blocked and reflected, signals with a shorter wavelength, however, are affected in their travel times (delay) and their signal paths (bending). The delay of the signal between the satellite S and the receiver R defined as is only an approximation as higher-order effects are neglected. The sign on the right-hand side changes whether it is applied for carrier phase observations ('-') or for pseudorange measurements ('+') [1]. In geodesy the delay can be interpreted in two ways, namely as a 1. disturbing quantity which has to be corrected for precise applications or as a 2. target quantity.
In the latter case, d ion can be used as observation to determine the electron density N e as the ionospheric key parameter. Among many others, e.g., ionosondes and radio occultations, GNSS is the most prominent observation technique used to gain information about N e . However, GNSS is not sensitive to the electron density directly, but on the integrated effect along the ray path. It is denoted as Slant Total Electron Content (STEC) and defined as STEC(x S , x R , t) = R S N e (x, t) ds . ( STEC depends on the electron density N e (x, t) along the ray path between the position x S of the satellite S and the position x R of the receiver R for time t. The position vector x reads x = r cos ϕ cos λ, cos ϕ sin λ, sin ϕ T , where ϕ and λ are latitude and longitude, r is the radial distance within a geocentric coordinate system Σ E .Many ionosphere models rely on the assumption that all electrons in the ionosphere are concentrated on a spherical shell Ω H of infinitesimal thickness with fixed height H = r − R e above the spherical Earth of radius R e . It is denoted as single-layer model (SLM) and allows for the transformation of STEC into the vertical direction as by applying a mapping function M(z) depending on the zenith angle z [2]. It yields vertical total electron content (VTEC) referred to the Ionospheric Pierce Point (IPP), the intersection point of the ray-path with the shell Ω H at the position x IPP . Several approaches exist for the realization of the mapping function; recently UPC developed the Barcelona Ionospheric Mapping Function (BIMF) [3]. Besides the electron density VTEC means another ionospheric key parameter which can be represented as the integration, of the electron density along the height between the lower boundary h l and the upper boundary h u of the ionosphere. There are several approaches for representing VTEC globally by using input data from GNSS. The Ionospheric Associated Analysis Centers (IAAC) of the International GNSS Service (IGS), namely (1) the Jet Propulsion Laboratory (JPL), (2) the Center for Orbit Determination in Europe (CODE), (3) the European Space Operations Center of the European Space Agency (ESOC), (4) the Universitat Politècnica de Catalunya (UPC), (5) the Canadian Geodetic Survey of Natural Resources Canada (NRCan), (6) the Wuhan University (WHU) and (7) the Chinese Academy of Sciences (CAS) routinely provide Global Ionosphere Maps (GIM) representing VTEC based on different mathematical approaches, namely the series expansion or the discretization technique. Whereas in the latter case, usually pixels or voxels are chosen [4], in the case of a series expansion spherical harmonics (SH) are mostly used. The GIM products of the IAACs are also used to generate an IGS GIM as a combined solution [5].
As a matter of fact, the representation of VTEC on a global scale suffers from the inhomogeneous distribution of the input data which does not allow a high-resolution representation on a global scale [6]. Capturing this inhomogeneity, the maximum degree of SHs is generally set to n = 15. The GIMs are usually provided in IONosphere Map EXchange (IONEX) format with a temporal resolution of 2 hours and a spatial sampling of 2.5 • × 5 • with respect to latitude and longitude and with a latency of more than one week [7,8].
Goss et al. [9] developed an approach based on B-spline basis functions and a Multi-Scale Representation (MSR) which allows to model VTEC globally with a higher spectral resolution. B-splines have already proven to be suitable basis function for representing global VTEC [10][11][12]. However, as a matter of fact, finer ionospheric structures can only be monitored and modeled where high-resolution input data are available, i.e., in continental regions such as Europe and North America, where a large number of stations of the global network of the IGS are located. Some of the continental regions are also covered with additional stations from regional networks, e.g., EUREF for Europe, AFREF for South Africa, ARGN and SPRGN for the Australian and Pacific region. Figure 1 shows all IPPs collected at 12:00 UT from the receiver stations of IGS and EUREF. Due to the coverage of observation stations over Europe, a dense cluster of observations is available. The Royal Observatory of Belgium (ROB) [13], for instance, provides Regional Ionospheric Maps (RIM) for the European region. Besides this, there are several approaches for the regional representation of VTEC, e.g., regional B-spline expansions [10,[14][15][16], the thin plate spline (TPS) interpolation [17], the Taylor series [18] and the SH expansion using a spherical cap coordinate system [19].
In this paper, we describe a procedure to represent VTEC for regions with dense data coverage with high spatial and spectral resolution based on a two-step strategy using B-spline basis functions. In the first step, we apply a global B-spline model with polynomial B-splines for latitude and trigonometric B-splines for longitude in order to represent the long-period VTEC part globally. In the second step, we perform a densification into regions with dense data coverage by using the global model as background and estimating B-spline coefficients of tensor products of polynomial B-splines representing the finer structures of VTEC. In combination, the model allows for a high-resolution VTEC representation within the densification regions. The developed approach is denoted as two-step model (TSM).
The paper is outlined as follows: in Section 2 a general description of the modelling of VTEC is given. Additionally in Section 2.1, the two-step model, based on B-spline basis functions and the connection between the two modelling steps is derived. Section 2.2 comprises a detailed description of the estimation of the unknown model parameters using a Kalman filter. Furthermore, in Section 2.3 the products of the TSM are defined. In Section 3 the TSM is applied to the European region during the St. Patrick storm days in March 2015 with a detailed validation in Section 4.

Modelling Strategies for the Regional VTEC Representation
, as introduced in Equation (5), can be represented by means of a series expansion as wherein φ k (x) are given space-dependent basis functions and c k (t) are the corresponding time-dependent series coefficients. Assuming that I s observations y(x i s , t s ) of VTEC are given at IPP positions P i s ∈ Ω H with i s = 1, 2, . . . , I s and discrete time moments t s = t 0 + s · ∆t with s ∈ N 0 and ∆t being a given temporal sampling interval, the observation equation can be written as with e(x i s , t s ) being the measurement errors of the observations. Note, in the sequel of this paper, if not stated otherwise, we omit additional unknown parameters such as the satellite and receiver biases, i.e., Differential Code Biases (DCB) for the geometry free observations as on the right-hand side of Equation (7).

Two-Step B-spline Expansion
At DGFI-TUM we rely on a two-step VTEC model (TSM) in terms of B-splines. The basic equation of polynomial B-spline functions φ k 1 ,J 1 (ϕ) and trigonometric B-spline functions φ k 2 ,J 2 (λ) representing VTEC w.r.t. latitude ϕ and longitude λ; for more information see [9,[20][21][22][23]. The shift parameters k 1 ∈ {0, 1, . . . , K J 1 − 1} and k 2 ∈ {0, 1, . . . , K J 2 − 1} for polynomial and trigonometric B-splines, define the position of the B-spline functions on the sphere Ω H . The total numbers K J 1 = 2 J 1 + 2 and K J 2 = 3 · 2 J 2 of B-spline functions depend on the B-spline levels J 1 and J 2 . In order to select appropriate numerical values for J 1 and J 2 , the two inequality chains have been considered by Goss et al. [9] and relate J 1 and J 2 to the average global sampling intervals ∆ϕ and ∆λ of the input data, i.e., the GNSS measurements, and the maximum degree n max of a SH expansion, which means the measure for the spectral content of a signal given on a sphere. For the numerical values 1 to 6 for the B-spline levels J 1 and J 2 , Table 1 presents the corresponding largest numerical values for n max as well as the corresponding sampling intervals ∆ϕ and ∆λ as already published in [9]. Additionally, the minimum wavelength L ϕ = 2 · ∆ϕ and L λ = 2 · ∆λ are given. Table 1. Values for the B-spline levels J 1 and J 2 , the maximum spherical harmonics (SH) degree n max , the input data sampling intervals ∆ϕ and ∆λ, as well as the minimum wavelengths L ϕ and L λ ; the left part presents the values along the meridian (upper inequalities in Equation (10)), the right part the corresponding values along the equator and its parallels (lower inequalities in Equation (10)). According to [9], the inequality chains (10) comprise two scenarios for determining J 1 and J 2 , depending on 1) the given sampling intervals ∆ϕ and ∆λ and 2) a given maximum degree n max of the spectral content. As mentioned before, most of the GIMs of the IAACs are based on a SH expansion up to degree n max = 15. Consequently, the level values J 1 = 4 and J 2 = 3 can be derived from Equation (10) [9]. Based on these values, average sampling intervals of approximately ∆ϕ = 10 • and ∆λ = 15 • are required.
Regional B-spline Model subject to restrictions: the smallest wavelength represented by the global background must be shorter or equal to the largest one of the regional model. Following this requirement, the inequalities Υ ≥ L ϕ and Γ ≥ L λ (12) can be set up. Similarly to the definition of the global level values, the inequality relations for the regional level values read in dependence of the average sampling intervals ∆ϕ and ∆λ of the observations within the region ∆Ω.
Numerical values from the inequality chains (13) for a specified region can be found in Table 2.
Following this concept, the question remains, how to avoid correlations between the two parts of the TSM. To solve this problem we divide the set D of all available observations into a subset D 1 of global observations y glob and a subset D 2 of observations y reg related to possible regional densification areas ∆Ω. In our current installation we determine the set D 1 by segmenting the Earth's surface into bins of a size related to the global sampling intervals ∆ϕ and ∆λ according to the Equations (10). For all bins in which more than two receiver stations of GNSS networks (e.g., IGS, EUREF, UNAVCO) are located, we select this receiver station which is the nearest one to the center of the bin. Then all observations of the chosen receiver station, related to the corresponding IPPs, are collected in the global set D 1 . The observations of the non-chosen stations within the bins are potential candidates for the regional data set D 2 . It should be noted that the densification data set D 2 may consist of several data sets D 2,a , D 2,b , . . . related to non-overlapping densification areas ∆Ω A , ∆Ω B , . . . such that Note, the realization of the two data sets D 1 and D 2 is currently based on a preliminary procedure and will not be discussed in more detail. There are more sophisticated approaches, i.e., based on Voronoi diagrams considering also the station and receiver qualities [24].

Estimation of B-spline Coefficients
The Kalman-Filter (KF, Kalman, 1960) is a sequential estimator which yields the ionospheric parameters for both modelling steps (global and regional). In the KF the input data from the past have not to be stored and the current state is updated as soon as new observations are available [25]. It is an optimal recursive estimator in terms of minimum variance estimation including a time update (prediction step) and a measurement update (correction step); see e.g., [25][26][27]. The approach consists of a linear formulation of the state equations and of the observation equations In Equation (14) β s is the u × 1 vector of unknown series coefficients at time t s predicted from the state vector β s−1 of the previous time step t s−1 by means of the u × u transition matrix F s and the u × 1 vector w s−1 of the process noise with u as the number of unknown series coefficients. Equation (15) consists of the I s × 1 vector y s = (y(x i s , t s )) of observations, the vector e s = (e(x i s , t s )) of corresponding measurement errors and A s the corresponding I s × u design matrix. The vectors e s and w s are characterized by the expectation vectors E(e s ) = 0 and E(w s ) = 0, as well as the covariance matrices where δ s,l is the delta symbol which equals to 1 for s = l and to 0 for s = l; see [9,28].

Kalman Filter
The solution of the estimation problem as defined by the Equations (14) and (15) can be found by sequential application of the prediction step and the correction step. Hence, the state equation (14) is used to obtain the predicted ionospheric target parameter vector β − s and its covariance matrix Σ − β,s as for details see [28]. Once the prediction of the state vector from time t s−1 to time t s is given, the predicted state vector and its covariance matrix are updated with the new allocated measurements at time t s by and where β and Σ β,s are the updated state vector and its covariance matrix. Thereby the so-called Kalman gain matrix behaves like a weighting matrix between the new measurements and the predicted state vector.

Coordinate System
VTEC exhibits a time-varying phenomenon. Therefore, a proper prediction model is required to take the time variation of the ionospheric parameters into account. Since VTEC glob is represented in the geocentric solar magnetic (GSM) coordinate system which results in much slower variations of the B-spline coefficients in time, a random walk approach could be performed for all the unknown parameters of the global model part within the KF. In this case the transition matrix is defined as identity matrix, i.e., F s = I. Hence, the predicted state vector reads β − s = β s−1 , according to Equation (17). Detailed information about the GSM coordinate system can be found, for instance, in [29,30]. The regional model, on the other hand, is set up in the Earth-fixed geographical coordinate system. Although VTEC varies significantly faster in this coordinates system, since it comprises the diurnal and shorter variations, too, the main part of VTEC is represented by the global background model VTEC glob . Thus, the estimated differences ∆VTEC reg of VTEC reg with respect to the background VTEC values exhibit small values. Thus, the random walk approach can also be used as a prediction model for the unknown regional B-spline coefficients.

Definition of the Global and the Regional Kalman Filter
Since both model parts are defined in different coordinate systems and based on different observation data sets, the estimation of their unknown parameters will be carried out separately. Therefore we rewrite the system of observation equations (15) for the global model as y s,glob + e s,glob = A s,glob β s,glob (22) and for the regional model as y s,reg + e s,reg = A s,reg β s,reg (23) in accordance to the different data sets D 1 for the global model and D 2 for the regional model. The global I D 1 s × 1 observation vector y s,glob consists of y(x i s , t s ) ∈ D 1 , the regional I D 2 s × 1 observation vector y s,reg in Equation (23) consists of the values for all observations y(x i s , t s ) ∈ D 2 . The vectors e s,glob and e s,reg are defined accordingly. The design matrices A s,glob and A s,reg are constructed by the tensor products of the B-spline functions according to the Equations (9) and (11). The state vector β s,glob consists of the unknown B-spline coefficients and receiver DCBs of the receivers that are assigned to D 2 , whereas the satellite DCBs can be omitted, since they are estimated within the global model part.
Since the estimation is performed for both model parts separately, we need to decompose the matrices Σ y and Σ w in the KF Equations (17) to (21) as well. The matrices Σ y,glob and Σ y,reg are set up in accordance to [28] and consist of two diagonal block matrices related to GPS and GLONASS observations of the data sets D 1 and D 2 . The relative weighting between the blocks is performed by adaptively defined variance factors, details can be found by [31]. The block matrices Σ w,glob and Σ w,reg are related to the process noise for the global and the regional model, respectively. They are adapted to the magnitude of the signal, i.e., we have to define different process noise for VTEC glob and ∆VTEC reg . Additionally, the distribution of the observation data is considered for the calculation of the process noise. For more information see [31].

VTEC Model Output
By applying the KF procedure according to the Equations (17) to (20) to the global observations of data set D 1 given in the GSM coordinate system in the first step and to the observations of the data set D 2 given in the Earth-fixed geographical coordinate system for the regional densification step, we obtain the two estimations β s,glob and β s,reg for the unknown series coefficients as well as their estimated covariance matrices Σ β,s,glob and Σ β,s,reg successively. From these estimates we can evaluate the V × 1 vector f s,glob =Ā s β s,glob and (25) as well as the corresponding covariance matrix. The matrixĀ s is set up in a similar way as the matrix A s,glob in the Equation (22). Accordingly, we can evaluate where f s,reg is a W × 1 vector and Σ s,reg the corresponding W × W covariance matrix. Assuming both model parts are uncorrelated, the total VTEC within the region can finally be calculated by applying Equation (8) as for latitude ϕ = ϕ v = ϕ w and longitude λ = λ v = λ w . The same applies for the calculation of the standard deviation where the values σ 2 s,reg (ϕ, λ) and σ 2 s,glob (ϕ, λ) for each point P(ϕ, λ) at time moments t s are given by the diagonal elements of Σ s,glob and Σ s,reg , respectively. A flowchart of the developed TSM approach is shown in Figure 2. Starting from global and regional GNSS station networks, the respective data sets D 1 and D 2 are defined. The left side of Figure 2 depicts the process chain for generating the global model part. The right-hand side of the flowchart shows the process for generating the regional densification. Here the extent of the densification area has to be chosen in dependency of the global levels J 1 and J 2 according to (12). Furthermore, the global model is required as a background for modelling the regional ∆VTEC. Figure 2. Flowchart of the two-step model (TSM), with the process chain for the first step-the global model-on the left side in yellow and the process chain for the second step-the regional densification-on the right side in green.
At the end of the process chain, the regional high-resolution VTEC reg values are calculated by combining the estimated global and regional model parts according to the Equations (29) and (30).

VTEC Products
According to [9], the previously explained procedure allows for the dissemination of two types of products, namely  (25) and (26) as well as b) estimated values ∆VTEC reg (ϕ w , λ w , t s ) and their standard deviations σ ∆VTEC reg (ϕ w , λ w , t s ) from Equation (29) on a regular grid of the regional model at time moments t s taken from the Equations (27) and (28) as well as c) estimated values VTEC reg (ϕ w , λ w , t s ) and their standard deviation σ VTEC,reg (ϕ w , λ w , t s ) on a regular grid of the regional model at time moments t s taken from the Equation (29) and (30) to the user. The two product types reflect the two strategies of dissemination. In terms of product type 2, VTEC is provided to users typically by means of the IONEX format, see [18]. Thereby, the sampling intervals ∆Φ in latitude and ∆Λ in longitude between the grid points P(ϕ v , λ v ) are usually chosen as 1 • , 2.5 • or 5 • . For the calculation of VTEC values at arbitrary points between the grid points an interpolation has to be applied which generally degenerates the quality of the calculated values [9].
In the case of product type 1, the VTEC values and the corresponding standard deviations can be calculated at arbitrary points with higher accuracy using the Equations (25) to (28) directly. However, then an encoder procedure for the estimated B-spline coefficients is necessary. In fact, the presented two-step approach described in Section 2.1 provides two products in terms of Product type 2, namely -TSM-Product 2a the global model VTEC glob (ϕ, λ, t) as well as the corresponding standard deviations σ glob (ϕ, λ, t) for the representation of VTEC on a global scale and -TSM-Product 2c the regional model VTEC reg (ϕ, λ, t) and the corresponding standard deviations σ reg (ϕ, λ, t) for the representation of VTEC with a higher spatial and spectral resolution within the selected region.
Due to the independent estimation of the coefficients d J 1 ,J 2 k 1 ,k 2 for the global and d J 3 ,J 4 k 3 ,k 4 for the regional model part, the global model is not affected by the densification step and can be disseminated as Product type 2a. The second step of the TSM, however, depends on the global background model. Hence, in this case both model steps have to be evaluated in order to provide the TSM-Product type 2c.

Two-Step Model for Europe
In the sequel, the developed approach is applied to real data collected during March 2015. This time span is selected because of the solar storm, the St. Patrick storm, happened on March 17. Thus, stronger spatial variations in VTEC can be observed and modeled. We apply the TSM with a densification for the European region. It is particularly suitable for the application of the TSM due to a dense coverage of receiver stations from both, the global IGS network and the regional EUREF network, which provides a dense observation distribution over Europe, see Figure 1. It shall be noted that the idea of the following investigation is the generation of a VTEC model defined within a region that is characterized by a high spatial and spectral resolution such that a TSM-Product 2c can be established. It is not our objective to create a VTEC model of the most highest spatial and spectral content for a given region.

Global Model
We follow the workflow for the TSM procedure presented in Figure 2 and start with the definition of the different data sets D 1 for the global modelling part. Figure 3 shows in the left-hand panel all available stations from the global IGS network (red dots) and from the regional open network in Europe, EUREF (blue dots). The right-hand panel depicts the selected stations for creating the set D 1 with a more homogeneous distribution over the globe Ω H . The stations are selected by means of a segmentation of the Earth's surface into bins of size 10 • × 10 • w.r.t. latitude and longitude, according to the inequality chains (10). The observations related to the IPPs of the chosen receiver stations define the data set D 1 . Figure 4 shows exemplarily the corresponding IPPs from GPS and GLONASS on 17 March 2015 at 12:00 UT with a less dense data coverage over the continental regions compared to the distribution of IPPs in Figure 1. According to the workflow in Figure 2 (left), the next step is to select the global B-spline levels J 1 and J 2 . As already mentioned in the context of Table 1 the levels J 1 = 4 and J 2 = 3 are chosen here. Given the vectors y s,glob , e s,glob , the matrices Σ y,glob , Σ w,glob , we run the Kalman filter for estimating the unknowns of the global model. We set the Kalman filter step size to ∆t = t s − t s−1 = 10 minutes and estimate the 18 · 24 = 432 global B-spline coefficients as well as the the corresponding standard deviations d 4,3 k 1 ,k 2 (t s ) k 1 =0,...,17, k 2 =0,...,23 , with K J 1 = 2 4 + 2 = 18 and K J 2 = 3 · 2 3 = 24. The top left panel of Figure 5 shows the estimated coefficients and the top right panel the corresponding standard deviations from the Equations (19) and (20), respectively. The estimated parameters are of Product type 1a as defined in Section 2.3. For the calculation of Product type 2a we select the output sampling intervals ∆Φ = 2.5 • and ∆Λ = 5 • w.r.t. the geographical latitude and longitude, resp., as well as the temporal output sampling interval ∆T = 10 minutes. The estimated grid values VTEC glob (ϕ v , λ v , t s ) and the corresponding standard deviations σ VTEC glob (ϕ v , λ v , t s ) are shown in the lower panels of Figure 5. Note, here we applied a bi-linear interpolation according to [9,18] for visualization.

Regional Model
As with the modelling of the global part, the first step for the regional model is to define the data set D 2 . The method for this has already been described at the end of Section 2.1. Figure 6 depicts the receiver stations for both data sets D 1 (red circles) and D 2 (blue stars) within the European region. The observations related to the IPPs of the receiver stations within the regional densification area are collected in the data set D 2 and exemplarily shown in Figure 7. Since the relation D 1 ∩ D 2 = ∅ holds, no observation can be a member of both data sets. Consequently, we assume negligible correlations between the two data sets. UT; the IPPs within the blue box are selected for the regional data set D 2 .
The regional model is set up after the definition of the size ∆Ω of the densification area. According to the inequalities (12) and (13) the regional levels J 3 and J 4 have to be defined independence of the observation distribution over Europe. We proceed as follows: 1.
Selection of the regional extent According to (12)

Selection of the B-spline levels
Although the sampling intervals in the center of the defined region ∆Ω are particularly small, the average sampling intervals ∆ϕ ≤ 2.5 • and ∆λ ≤ 2.5 • are chosen due to the larger gaps at the boundaries. Substituting ∆ϕ, ∆λ, Υ and Γ in the mid part of the chain (13) yields the maximum B-spline levels J 3 = 3 and J 4 = 3 from Table 2 which shows the numbers from the inequality chains (13) analog to Table 1 for the global case. Table 2. Values for the B-spline levels J 3 and J 4 , the maximum SH degree n max , the input data sampling intervals ∆ϕ and ∆λ, as well as the minimum wavelengths L ϕ = 2 · ∆ϕ and L λ = 2 · ∆λ; the left part presents the numbers along the meridian (upper inequalities in Equation 13), the right part the corresponding numbers along the equator and its parallels (lower inequalities in Equation (13)).  Table 1 we identified the global level values J 1 = 4 and J 2 = 3 with the SH degrees value n max = 17 and n max = 12, resp., as cutoff frequency of the spectral signal content of VTEC. From Table 2 we obtain in the same manner for the B-spline levels J 3 = 3 and J 4 = 3 of the regional model. Consequently, our final regional VTEC model contains apectral component up to the cutoff frequencies n max = 64 in latitude direction and n max = 47 in longitude direction.

Latitude
With the chosen level values J 3 = 3, J 4 = 3, the observation vector y s,reg , the measurement error vector e s,reg , the covariance matrix Σ y,reg of the observations and the covariance matrix Σ w,reg of the process noise, we perform the estimation of the regional B-spline coefficients by means of the Kalman filter. Its step size is set to ∆t = 10 minutes, as it was applied for the global model. We obtain the estimated coefficients and their standard deviations as with K J 3 = 2 3 + 2 = 10 and K J 4 = 2 3 + 2 = 10 as the total number of B-spline coefficients. As for the global model part a test of significance was carried out for the estimated coefficients. From the 100 B-spline coefficients around 60 to 70 are significant during night time and 70 to 75 at local noon, with a α-value of 0.01. During the storm time on March 17, the number of significant coefficients increases to 80, because of the stronger signal variations during that time.
For the example in Figure 8, we applied the Equations (27) and (28) to calculate ∆VTEC and their corresponding standard deviations on a grid as Product type 2b introduced in Section 2.3. We chose the output sampling intervals ∆Φ = 1 • in latitude and ∆Λ = 1 • in longitude. For the representation of the map in the upper panels of Figure 8 we use the bi-linear interpolation [9]. For the calculation of the total VTEC values and the corresponding standard deviations within the region ∆Ω we finally apply the Equations (29) and (30) using the estimated values -defined as Product type 2a -VTEC glob (ϕ, λ, t) and σ VTEC glob (ϕ, λ, t) at the same grid points with sampling intervals ∆Φ = 1 • and ∆Λ = 1 • . The lower panels represent finally the maps of Product type 2c with VTEC reg (ϕ, λ, t) as the total VTEC values (left) and σ VTEC reg (ϕ, λ, t), the corresponding standard deviations (right) within the European region for 17 March 2015 at 12:00 UT.
Due to the absolute magnitude of the global signal values with up to 120 TECU, the detailed structures of the regional model with magnitudes of ±4 TECU cannot be visualized appropriately in the bottom left panel of Figure 8. The standard deviation map on the bottom right panel shows larger values at the region boundaries and thus, reflects rather nicely the distribution of the observations, cf. Figures 4 and 7.

Validation
In order to assess the VTEC TSM Products 2a and 2c as defined at the end of Section 2, we use the so-called dSTEC analysis. It is based on real GNSS observations for the comparison with the estimated VTEC maps and features an accuracy of less than 0.1 TECU, [8]. The differences between the observations STEC(x S s , x s,R , t s ) according to Equation (2) at time moments t s and the observation STEC(x S re f , x re f ,R , t re f ) at the reference time moment t re f along the same arc is compared with the differences according to Equation (4). The reference time moment t = t re f is usually referred to the observation with the smallest zenith angle z = z re f , see Figure 9. The quality assessment is finally performed by evaluating the differences with expectation value E(dSTEC(t s )) = 0, for different products, e.g., the TSM-Products 2a and 2c, or the final GIMs of IGS [32], or its IAACs as it was done by Roma-Dollase et al. [7]. Within the dSTEC analysis the difference (35) has to be calculated for a set of a number of selected receiver stations. For the choice of the receiver stations three strategies exist: 1.
The observations of the selected receiver stations are used for the computation of all maps used in the analysis

2.
No station of the selected receivers stations contributed to any of the different products, thus, the stations are independent from all products

3.
A mixture of stations satisfying the items 1. and 2.
Since we define non-overlapping station lists providing observations to D 1 and D 2 , we follow the third strategy. Figure 10 shows the stations used in our dSTEC analysis. Note, we selected also stations that are located close to the boundaries of the region, to use their observations with IPPs located within the region. For the assessment of the product quality we use the following TSM-Products (cf. Goss et al. [9]):

TSM-Product 2a:
-otlg_a: the estimated VTEC We denote the global models as 'otlg' similar to Goss at al. [9]. The first letter 'o' refers to the OPTIMAP processing software, which was developed within a third-party project (see Acknowledgements) and means the basis of the presented study. The second letter indicates the chosen temporal output sampling interval with 't' for ∆T = 10 minutes. The third letter denotes the spectral content of the output signal; it is chosen as 'h' for high and 'l' for low. The last letter indicates the model domain with 'g' for global and 'r' for regional. In the case of the global model we indicate two different grid resolutions with their endings 'a' and 'b', cf. Table 3.
For the validation of our TSM-Products we use other products listed in Table 3. They are characterized by different extents, latencies and output samplings ∆Φ, ∆Λ and ∆T. We compare the global products 'otlg_a', 'otlg_b', 'codg' and 'uqrg', where 'otlg_a' and 'otlg_b' are based on ultra-rapid orbits and hourly GNSS observations with a latency of less than 3 hours [9]. The CODE product 'codg' is a final product and thus, based on post-processed data. UPC's product, denoted as 'uqrg', is a rapid product and based on rapid GNSS data with a latency of approximately one day. Our regional product 'othr' is an ultra-rapid product, too. The only comparable external product for Europe is a final one provided by the Royal Observatory of Belgium (ROB) with a latency of approximately one week.  Figure 11 shows the RMS values of the differences (35) during the time span between 8 March and 23 March 2015 for each product at the 10 receiver stations shown in Figure 10.  Figure 10 for the products listed in the legend; the two global products 'otlg_a' and 'otlg_b' as well as the regional product 'othr', the external regional product 'robg' and the external global products 'codg' and 'uqrg' of the IAACs The RMS values at the different stations vary approximately between 0.5 TECU and 0.9 TECU. Note, during the geomagnetic storm on March 17, the RMS values for all products at all stations increase and can reach values of up to 4 TECU for a short time, but are averaged in Figure 11. In order to rate the RMS values at different station and for the TSM products appropriately, we need to distinguish the stations used in the assessment, which are assigned either to D 1 or to D 2 . Hence, the stations 'BRST', 'HERT', 'GRAZ', 'WTZR' and 'WRST' are stations which are used in the modelling procedure of the global VTEC maps, whereas 'UZHL', 'ORID', 'CEBR', 'ZIMM' and 'TLSE' are used for the regional model part. At a first step we compare the two global products of DGFI-TUM, 'otlg_a' and 'otlg_b' given with different spatial output sampling intervals ∆Φ and ∆Λ. Consequently, different RMS values are obtained for both products at all stations. For the values VTEC(x s,IPP , t s ) and VTEC(x re f ,IPP , t re f ) in Equation (34) a bi-linear interpolation has to be applied between the given grid points. The quality of the interpolated VTEC decreases with an increasing spatial resolution and thus, we achieve a slightly smaller overall RMS value for the product 'otlg_b'. Comparing the bars for 'otlg_a' and 'otlg_b', we find larger differences in the RMS at the stations 'UZHL', 'ORID' and 'CEBR', which are independent stations for both models. The RMS values at all other stations are approximately equal.
In the next step, we compare the global product 'otlg_b' and the regional product 'othr'. Both are given with the same spatial and temporal sampling intervals. From the orange colored and the red colored bars for 'othr' and 'otlg_b', respectively, we can conclude, that the regional model is of higher accuracy at all stations. This means, that we can achieve a significant improvement in the modelling accuracy within the region by applying the second step of the TSM. By comparing the overall RMS of 'othr' and ''otlg_b', given in the parenthesis, we even find an increase in the accuracy of around 18.3% for 'othr'.
In order to rate the TSM-Products 2a against products of the IAACs, we compare the overall RMS values, of 'otlg_a' with 0.71 TECU, otlg_b' with 0.70, 'codg' with 0.72 TECU and 'uqrg' with 0.76 TECU. Consequently, we find comparable results for the TSM products and 'codg' at all stations. However, when comparing the green and red color bars with the dark blue color bar for 'codg', for example at the stations 'CEBR' and 'WTZR', the differences are larger, so that local quality differences can be inferred. The product 'uqrg' performs slightly worse within the European region and for the selected time span, compared to the other global products.
Finally, we compare the regional TSM-Product 2b and the product 'robg'. They are given with different spatial and temporal sampling, cf. approximately the same area. However, 'robg' is a final product, whereas 'othr' is an ultra-rapid product. The product 'othr' shows a slightly higher accuracy than 'robg' of about 0.3 TECU in the overall RMS values. Nevertheless, both products are approximately of the same accuracy.

Summary and Conclusions
The presented approach aims at modelling VTEC with high spectral, spatial and temporal resolution for regions of dense data coverage. Compared to the GIMs provided by the IAACs, which are typically of lower spectral resolution, regional maps can represent finer structures. The global maps are usually provided with spatial output sampling intervals of ∆Φ = 2.5 • and ∆Λ = 5 • and a temporal sampling interval of ∆T = 2h. Regional maps, on the other hand, are usually provided with a finer grid resolution and a temporal sampling of ∆T ≤ 1 hour.
By means of the developed TSM, we generate both, a global model representing the coarser structures and a regional model representing the finer structures of VTEC. Both model parts are adapted optimally to the distribution of the independently defined data sets. The global data set consists of observations collected from selected receiver stations of the IGS network. To be more specific, the selection of these stations has to be adapted to the global situation. This means, that the observation sampling of the global data set has to fulfill the inequalities (10). Depending on the coverage the remaining receiver stations may provide dense clusters of observations within specific regions, cf. Figure 6. We performed numerical tests for Europe for March 2015, including the St-Patrick storm day, the strongest space weather event within the last decade. With the global resolution levels J 1 = 4 and J 2 = 3 we generate as the first step of the TSM the global VTEC model, comparable to the GIMs provided by the IAACs. In the second step, the regional VTEC model with resolution levels J 3 = 3 and J 4 = 3 is created. We apply the developed TSM to a densification area in Europe with cutoff frequencies of n = 64 in latitude direction and n = 47 in a longitude direction. Finally, the applied dSTEC validation shows that our generated regional ultra-rapid model with a latency of less than three hours is of higher accuracy than the global models and comparable to the regional final product of ROB.
Within this paper, we showed, that a higher spectral resolution and a higher accuracy can be achieved using appropriate regional modelling approaches. These high-resolution products can be used to study small scale and fast-changing variations of the ionosphere, e.g., traveling ionospheric disturbances (TID) [35,36], as well as plasma irregularities, such as plasma bubbles. For the use of the regional model in modern GNSS applications, such as autonomous driving or precision agriculture, the developed approach can be seen as an intermediate step towards the following developments: -Since autonomous driving or precise navigation requires real-time information about the state of the ionosphere, the TSM approach must be adapted accordingly. For this purpose, the latency (see Table 3) of the TSM products must be reduced significantly. There are two ways, namely (1) transferring the products to real-time using a suitable forecast algorithm or (2)-which is more preferably-directly applying real-time GNSS observations.
-There is also the question of how ionosphere information should be disseminated to the user in real-time. In general, the dissemination can be undertaken in two different ways: based on estimated series coefficients (Products type 1) or on estimated grid values (Product type 2). In the latter case, an interpolation of the grid values has to be performed to obtain VTEC information at any arbitrary point within the area on the investigation. In case of dissemination of estimated series coefficients, the Radio Technical Commission for Maritime Services (RTCM) message has to be mentioned, which allows us to transfer of the SHs coefficients up to degree 16 to the user. It is in the current form restricted to the sole use of SHs and thus, not suitable for the high-resolution products developed within this paper. Consequently, the dissemination of the products needs to be discussed urgently in order to provide highly precise information to the users.
In a further step, a validation method for ionospheric maps will be developed that is more suitable for navigation applications. A more sophisticated method will allow the determination of positioning accuracy using ionospheric corrections from different products.