The main advantages are from the engineering point of view (as @Alexey mentioned). In the widely used Kriging procedure you can interpret your own "space" by providing a "correlation" (or covariance) model (typically called the variogram ellipsoid) for relations depending on distance and orientation.
There is nothing that prevents other methodologies to have the same features, it just happened that the way kriging was first conceptualized had a friendly approach to people that were not statisticians.
Nowadays with the rise of geostatistics based stochastic methodologies, like Sequential Gaussian Simulation among others, these procedures are getting used in sectors where it is important to define the uncertainty space (which can take thousands to millions of dimensions). Again, from the engineering point of view, geostatistics based algorithms are very easy to include in genetic programming. As so when you have inverse problems you need to able to test multiple scenarios and test their adaptability to your optimization function.
Let's leave the pure argumentation for a moment a state the facts for a modern real example of this use. You can either sample underground samples directly (hard-data) or make a seismic map of the subsurface (soft-data).
In hard data you can measure a property (let's say acoustic impedance) directly without(ish) error. The problem is that this is scarce (and expensive). On the other hand you have the seismic mapping which is literally a volume, pixel-wise, map of the subsurface but does not give you acoustic impedance. For simplicity purposes let's say it gives you the ratio between two values of acoustic impedance (top and bottom). So a ratio of 0.5 could be a division of 1000/2000 or 10 000/ 20 000. It's a multiple solution space and several combinations will do but only one accurately represents reality. How do you solve this?
The way seismic inversion works (the stochastic procedures) is by producing plausible (and this is another story all together) scenarios of acoustic impedance (or other properties), transform those scenarios into a synthetic seismic (like the ratio in the previous example) and compare the synthetic seismic against the real one (correlation). The best scenarios will be used to produce even more scenarios, converging into a solution (this is not as easy as it seems).
Taking this into account and speaking from the point of view of usability I would answer your questions the following way:
1) What makes them popular is usability, flexibility in implementation, a good number of research centers and institutions that keep making newer and more adaptable gaussian based procedures for several different fields (particularly in geosciences, GIS included).
2) The main advantages are, as mentioned before, usability and flexibility from my point of view. If it's easy to manipulate and easy to use you just do it. There are no particular features in gaussian processes that are not reproducible in other methodologies (statistics or otherwise).
3) They are used when you need to include more information into your model than just the data (information such has space wise relations, statistical distributions, and so on...). I can assure that if you have lot's of data with an isotropic behavior using kriging is a waste of time. You can get the same results using any other method that by requiring less information, its faster to run.