The GNSS tomography provides the opportunity to estimate atmospheric wet refractivity, which is important for precise positioning and navigation, as well as for constraining weather and climate models. The GNSS-derived water vapor is often estimated by implementing voxel-based inversions, but this technique is numerically unstable due to the high number of unknown parameters especially in networks covering large areas or when the spatial resolution increases. To mitigate this problem, in this study, we introduce functional based 3 and 4 dimensional (3D and 4D) inversion formulations. Here, the horizontal changes are modelled by the spherical cap harmonic functions. For the vertical component, empirical orthogonal functions derived from ERA5 (Empirical Reanalysis Fifth generation) are applied as background model. The time-dependency is accounted for by applying polynomial spline functions. Numerical results are based on a network of about 190 GPS stations in Germany during 15 days in summer and winter of 2018. Observations from 8 radiosonde stations are applied for comparisons. Our results indicate that the functional tomography is effective and retrieves wet refractivity indices with the mean RMSE (Root Mean Square Error) of about 1.9 ppm. These values are found to be up to 22% smaller than those derived by comparing ERA5 with the radiosonde data.