The main objective of this thesis is to develop a procedure for modeling and predicting ionospheric Total Electron Content (TEC) for high precision differential GNSS applications. As the ionosphere is a highly dynamic medium, we believe that to have a reliable procedure it is necessary to transfer the high temporal resolution GNSS network data into the spatial domain. This objective led to the development of a recursive physics-based model for the regular TEC variations and an algorithm for real-time modeling of the medium-scale Traveling Ionospheric Disturbances (MS-TID). The research described in this thesis can roughly be divided into three parts. The main application of these developments can be found in Network RTK. Network-RTK is a technique based on a network of reference receivers to provide cm-level positioning accuracy in real time for users in the field. To get centimeter accuracy after a short (minutes) initialization period the ionospheric delay for the user's receiver needs to be predicted very precisely between the ionospheric pierce points of the reference receivers at the double difference level. Having the cm-level accuracy in the ionospheric interpolation is crucial for the carrier phase ambiguity resolution by the user. To achieve high precision in the ionospheric interpolation, regular and irregular variability of TEC in time and space should be taken into account. The regular TEC variation, which can reach several hundreds TEC units, is mainly a function of solar zenith angle. The irregular (or non-repeatable) variations are mainly wavelike effects associated with Traveling Ionospheric Disturbances (TID). Although TID effects on the TEC are of the order of 0.1 TEC unit, MS-TIDs, with a typical wavelength less than a few hundred kilometers, is one of the main obstacles for accurate spatial interpolation of ionospheric induced delays in a medium-scale reference GPS network. Since most of interpolation methods either use spatial linear (or quadratic) interpolation or fit a lower-order surface, the methods are not capable to model the phase-offset, caused by MS-TIDs, at distinct ionospheric pierce points. There are two major complications. Firstly, interpolation must be done at the double-difference level, which involves taking single differences between ionospheric delays for the same satellite between two different receivers, followed by differencing single differences for different satellites. This means that two different patches of the ionosphere are involved, each related to a different satellite, and each possibly associated with different TIDs. Secondly, for operational network RTK, a real-time strategy for TID detection and modeling is needed. In the first part the performance of several empirical ionosphere models for the regular TEC variation, such as Klobuchar, NeQuick, and the IGS Global Ionosphere Maps (GIM) are studied in the mid-latitude region using GPS data. Our results show that the GIM was able to correct the absolute slant ionospheric delay to better than 80% under different geomagnetic conditions of the ionosphere. The NeQuick model, which performed better than the Klobuchar model, could correct about 60% of the slant ionospheric delay. NeQuick is a real-time ionospheric correction model for the future European Galileo navigation system. A key input parameter for NeQuick is the effective ionization parameter (Az), which will be provided as a second order polynomial in the Galileo broadcast message to single-frequency users. The coefficients of the polynomial will be estimated daily from at least 20 permanent Galileo monitoring stations. As Galileo is under development, we propose an alternative approach for estimating Az using Global Ionospheric Maps (GIM). The main advantages of the alternative approach over the standard approach are: (1) the alternative approach is more reliable, because, each IGS GIM is based on data of up to 300 GNSS stations world-wide and each IGS GIM is the combination of results of up to four analysis centers, (2) the coefficients are more representative for all regions on the world because they are computed from a world-wide grid instead of about 20 distinct locations, (3) with the alternative procedure it is possible to provide Az in a different representation, for instance using a higher order polynomial, grid, or other function types, and (4) the computational effort is much smaller assuming the IGS GIMs have already been computed. In the second part a normal ionosphere is defined using Chapman's ion production theory to approximate the regular variability of the Earth's ionosphere. The normal ionosphere consists of lower and upper region. The lower region is formed in a photochemical equilibrium resulting in a Chapman layer. The upper region is formed in a diffusive equilibrium, whilst ignoring the geomagnetic field, resulting in a new Chapman like ionospheric layer. Integration of the continuity equation of the normal ionosphere over height leads to a Boundary Value Problem (BVP) for the temporal evolution of VTEC. Solution of the BVP results in a novel recursive model for the regular TEC variation as a function of solar zenith angle. The main motivation for developing this model is that the empirical models of the first part were either ill-suited or too complicated to model and predict the regular variation of TEC for high precision differential GNSS applications. The performance of the new model is tested at local and global scales using GIM. In general, despite the geomagnetic field was ignored, the cases analyzed show that the model gives a good overall representation of the regular variation of VTEC in the mid-latitude region under a geomagnetically quiet ionosphere. This is an important result that shows the potential of the model for a number of applications. Since the model has a recursive form it is ideally suited to use as time update equation in a dynamic data processing or Kalman filter. Another application is to use it for removing the geometry-dependent trend from time series of GPS-provided ionospheric delays to provide a pure TID observation, which is carried out in the third part of this thesis. In the third part, a new algorithm for the real-time detection and modeling of MS-TID effects is developed. In order to eliminate effects from large-scale TIDs, the algorithm uses between-receiver single-difference (SD) ionospheric delays in a medium scale GPS network. Although single-differencing also eliminates to some extend the geometry-dependent trend, the remaining part cannot be neglected. In this thesis, we fit the SD data to the recursive model which was developed in the second part of the thesis. Any wavelike fluctuations in the data with respect to the model are assumed to be from MS-TID effects. The detrended SD data are the main input of the algorithm. The algorithm consists of six steps: initialization, detection, scraping, cross-correlation, parameter estimation, and ending. A MS-TID is assumed to be a planar longitudinal traveling wave with spatially independent amplitude that propagates in an ionospheric patch. All characteristic parameters of the MS-TID wave (e.g. period, phase velocity, propagation direction, and amplitude) are considered to be time dependent, while the Doppler-shift caused by the satellite motion is taken into account in the estimation step. The performance of the algorithm is tested with GPS data from a network. Although real TIDs are not perfect waves, the algorithm was able to model (in time and in space) the MS-TID to a large extend. The performance was found to be comparable with the Kriging interpolation method. This is an important first result, in part because these two methods are based on different principles, but also because there is still room for improvement in our algorithm. With our physics based model it is possible to avoid the planar wave approximation and take the phase-offset of the wave into account, something which is not possible with Kriging.